Signal detection algorithm in the frequency domain

Question

I would like to detect signals in the frequency domain. That is, having a spectrogram, automatically determine the signal in it (its frequency and bandwidth). I tried to apply the algorithm described here Peak signal detection in realtime timeseries data, but I get incorrect data (noise and incorrect signal values ​​are detected). An example of a spectrogram with which I work is in the picture. The detection of broadband signals and signals with frequency modulation is especially interesting. Perhaps someone will suggest similar implementations in python. Thank you in advance. This is how I process the spectrum:

    self.L = 15  # L-point filter
    self.b = (np.ones(self.L)) / self.L  # numerator co-effs of filter transfer function
    self.a = np.ones(1)  # denominator co-effs of filter transfer function
    fft_data = np.fft.fft(balanced_signal, n=self.fft_size) / self.fft_size
    _fft_log = (np.abs(np.fft.fftshift(fft_data)))
    _fft_log = ss.lfilter(self.b, self.a, _fft_log)
    spectrum = savgol_filter(_fft_log, 100, 5,mode='nearest')

Then I pass this data to the signal detection function expecting to get 1 where there is a signal and 0 where there is a noise value spectrum of mine signal:

    lag = 30
    threshold = 5
    influence = 0
 def thresholding_algo(self,y, lag, threshold, influence):
    signals = np.zeros(len(y))
    filteredY = np.array(y)
    avgFilter = [0] * len(y)
    stdFilter = [0] * len(y)
    avgFilter[lag - 1] = np.mean(y[0:lag])
    stdFilter[lag - 1] = np.std(y[0:lag])
    for i in range(lag, len(y)):
        if abs(y[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1]:
            if y[i] > avgFilter[i - 1]:
                signals[i] = 1
            else:
                signals[i] = -1

            filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i - 1]
            avgFilter[i] = np.mean(filteredY[(i - lag + 1):i + 1])
            stdFilter[i] = np.std(filteredY[(i - lag + 1):i + 1])
        else:
            signals[i] = 0
            filteredY[i] = y[i]
            avgFilter[i] = np.mean(filteredY[(i - lag + 1):i + 1])
            stdFilter[i] = np.std(filteredY[(i - lag + 1):i + 1])

    return dict(signals=np.asarray(signals),
                avgFilter=np.asarray(avgFilter),
                stdFilter=np.asarray(stdFilter))

But I get detected noises. As soon as I raise the threshold to 6, I have no detection at all.

Answer 1

For noisy signal like this is best to use sliding average and threshold that or diference between it and original signal. Which method to use depends on what you want to achieve and what you consider to be useful signal so either:

if ( abs(signal - sliding_average(signal)) >= positive_threshold )  this_is_signal;

or

if ( sliding_average(signal) <= negative_threshold )  this_is_signal;

Note that your values are negative that is why I used negative threshold and <= in the second example.

Also I hope the spectrum is power spectrum (You know FFT produce complex domain results so you have to convert to real domain first or use different frequency domain conversion like band pass filters DCT or DST in order to use real domain approaches)

In your case (signal with "horizontal base line") I would use the second option... Not a python coder but I tried it in C++/VCL with sliding average window size of 11 samples and threshold -117.0 the result:

In yellow is your "original" input (extracted from image) , the sliding average is in aqua the threshold value is in dark green and threshold result is in green.

Here C++/VCL code I did this with:

//$$---- Form CPP ----
//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "win_main.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
int xs,ys;
Graphics::TBitmap *bmp;
float *data_in=NULL,*data_out=NULL;
int n;
//---------------------------------------------------------------------------
float* load_graph(int &size,float y0,float y1,AnsiString filename)
    {
    int x,y,xs,ys;
    DWORD **pyx,c;
    float a,my,*data=NULL;
    Graphics::TBitmap *bmp;
    bmp=new Graphics::TBitmap;
    bmp->LoadFromFile(filename);
    bmp->HandleType=bmDIB;
    bmp->PixelFormat=pf32bit;
    xs=bmp->Width;
    ys=bmp->Height;
    size=xs;

    data=new float[size];
    pyx=new DWORD*[ys];
    for (y=0;y<ys;y++) pyx[y]=(DWORD*)bmp->ScanLine[y];
    my=(y1-y0)/float(ys-1);

    for (x=0;x<xs;x++)
     for (data[x]=y0,y=0;y<ys;y++)
      if ((pyx[y][x]&0x00FFFFFF)==0x00FFFF00)
        {
        a=y0+(float(ys-1-y)*my);
        data[x]=a;
        break;
        }

    delete[] pyx;
    delete bmp;
    return data;
    }
//---------------------------------------------------------------------------
void draw()
    {
    int i;
    float x,y,dx,y0,y1,dy;
    bmp->Canvas->Brush->Color=clBlack;
    bmp->Canvas->FillRect(TRect(0,0,xs,ys));

    y0=-123.0;
    y1= -70.0;
    dx=float(xs)/float(n);
    dy=float(ys)/float(y1-y0);

    // render input data
    bmp->Canvas->Pen->Color=TColor(0x00008080); // yellow
    for (x=0.0,i=0;i<n;i++,x+=dx)
        {
        y=(y1-data_in[i])*dy;
        if (!i) bmp->Canvas->MoveTo(x,y);
        else    bmp->Canvas->LineTo(x,y);
        }

    // render output data
    bmp->Canvas->Pen->Color=TColor(0x00F0F000); // aqua
    for (x=0.0,i=0;i<n;i++,x+=dx)
        {
        y=(y1-data_out[i])*dy;
        if (!i) bmp->Canvas->MoveTo(x,y);
        else    bmp->Canvas->LineTo(x,y);
        }

    // render threshold
    float thr=-117.0;
    bmp->Canvas->Pen->Color=TColor(0x00008000); // dark green
    y=(y1-thr)*dy;
    bmp->Canvas->MoveTo(0,y);
    bmp->Canvas->LineTo(xs-1,y);
    bmp->Canvas->Pen->Color=TColor(0x0000F000); // green
    for (x=0.0,i=0;i<n;i++,x+=dx)
        {
        y=ys>>3;
//      if (fabs(data_in[i]-data_out[i])>=5.0) y-=ys>>4;
        if (fabs(data_out[i])<=-thr) y-=ys>>4;
        if (!i) bmp->Canvas->MoveTo(x,y);
        else    bmp->Canvas->LineTo(x,y);
        }

    Form1->Canvas->Draw(0,0,bmp);
    bmp->SaveToFile("out.bmp");
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    bmp=new Graphics::TBitmap;
    bmp->HandleType=bmDIB;
    bmp->PixelFormat=pf32bit;

    data_in=load_graph(n,-123.2,-70.0,"noise.bmp");
    data_out=new float[n];

    int i,j,k,m=5; // sliding average
    float a;
    for (i=0;i<n;i++)
        {
        for (a=0.0,j=i-m;j<=i+m;j++)
            {
            k=j;
            if (k< 0) k=0;
            if (k>=n) k=n-1;
            a+=data_in[k];
            }
        data_out[i]=a/float(m+m+1);
        }

    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    delete[] data_in;
    delete[] data_out;
    delete bmp;
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    xs=ClientWidth;
    ys=ClientHeight;
    bmp->Width=xs;
    bmp->Height=ys;
    draw();
    }
//-------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    draw();
    }
//---------------------------------------------------------------------------

And input Image I used for data extraction:

Most of the code is not important just look at the stuff marked with these comments:

// render threshold
// sliding average

the data_in[] is your input data, data_out[] is sliding average (with m as its half strength) and thr is the threshold.