Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

I'm currently designing a receiver which has to determine whether a signal contains specific frequencies or not. The frequencies are at constant 215Hz difference:

{16,351  16,566  16,781  16,996  17,211   17,426  17,641  17,856  18,071}

I sample at 44100kHz and I have 1024 samples each time. Therefore my resolution is 44100/1024=43.06Hz, and the bins centers are multiples of 43.06:

Bin[1]=43.06Hz,... BIN[385]=16578.10    BIN[386]= 16621.16   
BIN[387]=16664.22    BIN[388]=16707.28     IN[389]=16750.34
BIN[390]=16793.40    BIN[391]=16836.46    BIN[392]=16879.52    
BIN[393]=16922.58    BIN[394]=16965.64    BIN[395]=17008.70   
BIN[396]=17051.76    BIN[397]=17094.82,....

My algorithm is based on comparing the values of the desired bins (or the nearest ones) to the values of the bins in between. for example, if I wish to know if 16,566 exists or not, I compare BIN[385] to BIN[386],BIN[387],BIN[384],BIN[383] and if it's value is significantly larger (say by a hundred) then the average value of the side ones then I conclude that this frequency does exist.

My question is as follows: considering the fact that I work in high noise environment, and the algorithm I use is the one above, what would be the best window for me to apply before the FFT, other than the obvious rectangular window I currently use.

Thanks

share|improve this question
    
Do you have to use an FFT? If not I suggest using a correlation receiver which would yield optimal performance for Gaussian noise. –  Deve Jan 16 '13 at 9:33
add comment

3 Answers

up vote 3 down vote accepted

I'm presuming that you have some or all 9 of the input frequencies available at one time, and you just want to find out which of them is present (you didn't specify in your problem statement).

Two methods come to mind: 1) compute DFTs at the frequencies of interest, or 2) use the Gaussian ratio technique, as described in reference_1 and reference _2.

Regarding 1) above: When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer. In other words, you can compute a DFT bin centered at ANY frequency, not just those that correspond to the usual k = integer of your FFT. For example, given your sample rate and N, a frequency of F = 16351 would correspond to a DFT bin computed at k = (N*F)/sample_rate = 1024*16351/44100 = 379.669 (ie: between bins 379 and 380 of your FFT).

So for the DFT, X(F) = (sum over n){x[n]*e[-j*twopi*n*k/N]}, where X[F] is the complex result, N = 1024, F = 16351, k = 379.669, and n = 0 to N-1.

Some of the other values for your problem are:

Frequency and value of k

16,566 ==> 384.662

16,781 ==> 389.659

16,996 ==> 394.646

I'll leave it to you to figure out the rest. And, as you can see, you've got about 5 bins of separation between the input frequencies. That may, or may not be enough based on your SNR and method used. And you might want to use N = 1025 instead – it makes your input frequencies much more 'bin-centered.'

Method 2) above uses an FFT, which, given the right circumstances, can be extraordinarily accurate when computing frequencies. Note that in the first reference above, by using a Gaussian window, and then taking the log ratio of amplitudes from adjacent bins, you are actually solving for frequency, subject to the requirements of the technique. Some code (DevC++ compiler) for your specific problem, followed by the programs' output is shown below (explanation to follow):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft_recur(long, double* r, double* i);
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O
const long N = 1024;   double sample_rate = 44100., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ;
long n;      double  t, twopi = 2*3.141592653589793238;
double amp_sq[N], F, fbin, amp, power, c = 2.3;

for (n = 0; n < N; n++)  {  // generate test data
    t = twopi*n/sample_rate;
    r[n] = 1.*cos(16351.*t)+1.*cos(16996.*t)+1.*cos(17211.*t)+1.*cos(17426.*t) ;
} // end for

//Gaussian window as computed in [MCEACH94] program
double z;
for (n = 0; n < N/2; n++) {
    z = (.5+(N/2)-1.-n)*(twopi/2.)/N;
    w[n] = exp(-c*z*z);
    w[N-1-n] = w[n];
}

//multiply input data by window points
for (n = 0; n < N; n++)  {
    r[n] = w[n]*r[n];    // we're overwriting r[n] with the windowed r[n]
} //end for

fft_recur(N, r, i);       // note: fft is not scaled

// compute amplitude squared
for (n = 0; n < N/2+1; n++)  {
    amp_sq[n] = r[n]*r[n]+i[n]*i[n];
} //end for

// compute natural log of FFT outputs
for (n = 0; n < N/2+1; n++)  {
    r[n] = .5*log( amp_sq[n] ) ; // r[n] is now the natural log of the amplitude of FFT bin 'n'
} // end for                     // .5*log(A squared) = log(A)

cout<<"\n\nfrequency and amplitude estimation results\n\n n        magnitude      frequency       amplitude\n\n";
for (n = 370; n < 415; n++)  {
    fbin = n+.5+.5*c*(r[n+1]-r[n]);   // r[n+1] and r[n] are ln(b) and ln(a)
    F = (sample_rate/N)*fbin;
    amp = sqrt(c)*exp(r[n]+((fbin-n)*(fbin-n))/c)*(2.0*sqrt(twopi/2.)/N) ;
    printf("%2d\t%9.5f\t%9.5f\t%9.5f\n",n,amp_sq[n],F,amp);
} //end for

system ("PAUSE");
return 0;
} // end main
//******************** fft_recur ***********************
void fft_recur(long N, double *r, double *i)  {
long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1;
double c = -1.0, s = 0.0, t1, t2, u1, u2;

for (h = 0; h < N-2; h++) {    // ***** bit reverse starts here ***
    if (h < j) {
       t1 = r[h]; r[h] = r[j]; r[j] = t1; 
       t2 = i[h]; i[h] = i[j]; i[j] = t2;
    }
    k = i2;
    while (k <= j) {
          j = j - k;      k = k/2;
    }
    j = j + k;
}   //****** bit reverse done ******

for (k = 1; k < N; k = k*2) {
    l1 = l2;     l2 = l2*2;
    u1 = 1.0;    u2 = 0.0;
    for (j = 0; j < l1; j++) {
        for (h = j; h < N; h = h + l2) {
            i1 = h + l1;
            t2 = (r[i1] - i[i1])*u2 ;
            t1 = t2 + r[i1]*(u1 - u2) ;
            t2 = t2 + i[i1]*(u1 + u2) ;
            r[i1] = r[h] - t1;
            i[i1] = i[h] - t2;
            r[h]  = r[h] + t1;
            i[h]  = i[h] + t2;              
        } // end for over h
        t1 = u1 * c - u2 * s;
        u2 = u1 * s + u2 * c;
        u1 = t1;  //x = u1 - u2;    y = u1 + u2;
    } // end for over j
    s = - sqrt((1.0 - c) / 2.0);
    c =   sqrt((1.0 + c) / 2.0);
} // end for over k

} // end fft


frequency and amplitude estimation results

 n        magnitude      frequency       amplitude

370       0.01156       15959.53142       0.00065
371       0.01328       16002.86947       0.00070
372       0.01542       16046.20441       0.00076
373       0.01810       16089.42371       0.00082
374       0.02138       16136.92108       0.00095
375       0.03020       16209.62311       0.00211
376       0.14116       16370.55303       3.20360
377      77.00512       16349.81449       0.95808
378     3193.91932      16351.21160       1.00396
379     24622.71009     16350.90700       0.99982
380     32939.04014     16351.10559       0.99850
381     7803.22610      16350.67285       1.01020
382     319.11953       16353.81879       0.87280
383       2.60322       16324.16039       7.57917
384       0.00113       16605.10228       0.00051
385       0.00724       16606.09917       0.00052
386       0.00851       16650.89175       0.00058
387       0.01072       16694.51462       0.00065
388       0.01381       16737.66957       0.00074
389       0.01787       16785.61197       0.00090
390       0.02814       16862.78585       0.00251
391       0.17568       17013.74935       2.83597
392      85.58946       16994.84074       0.95935
393     3413.62331      16996.22071       1.00408
394     25288.20406     16995.88948       0.99974
395     32472.56144     16995.79357       1.00043
396     7296.89833      16957.17339       4.09340
397      60.55471       17218.13259       1.24775
398     3330.74781      17212.33274       1.01810
399     25463.82887     17210.95691       0.99978
400     32350.54927     17210.69761       1.00169
401     7144.91966      17170.91626       4.35710
402      55.60789       17436.25628       1.40301
403     3423.65695      17427.15709       1.01474
404     25642.98140     17426.03425       1.00002
405     32244.80235     17425.93874       1.00073
406     7095.46971      17426.35995       0.98918
407     278.99401       17421.20777       1.27242
408       1.56515       17494.25360       0.02622
409       0.02946       17619.31213       0.00091
410       0.01521       17677.62446       0.00071
411       0.01452       17719.83432       0.00069
412       0.01340       17762.64014       0.00066
413       0.01224       17805.61294       0.00063
414       0.01113       17848.67454       0.00060

First, 4 tones of the input are generated (I use twopi when computing time because it's easier than writing it 4 times into the cosines). Then, the Gaussian window is generated, after which it is applied to the input (by the way, the window shown here, which is from reference_1, is slightly wrong – I don't want to use the 'corrected' version of it because I don't want to write a page describing how I derived it). Then we take a 1024 point FFT, followed by the (non-normalized) calculation of amplitude squared. Then comes the calculation of the natural log of each FFT bin. Finally, I compute the frequencies and amplitudes via the equations in reference 1 and reference 2above.

As you can see from the results following the code, we have an isolated peak between bin 379 and 380 from the first input frequency, which is estimated as 16350.90700 and an amplitude of .99982. We can also see that the other 3 inputs exhibit a little 'peak merging' from bins 392 to 407, but the peaks and associated estimated frequencies and amplitudes are quite accurate.

The results could be even more accurate if I had: 1) adjusted parameters, 2) used the 'correct' window, 2) used a non-recursive FFT, or 3) used one of the variations of the technique (which I am not going to show here).

You might also notice that, as your value of the parameter 'c' gets smaller ('c' is related to the 'beamwidth' of the Gaussian window), you can actually improve some of your estimates, but it can cause some results to blow up due to numerical sensitivities. The effective 'Q' or 'gain' of the Gaussian ratio technique is on the order of 5,000 – 10,000, and it can be much higher if you know how to adjust the parameters properly.

Of course, if your SNR is too low, you may have to increase N to get more integration time, subject to the limits of signal stability.

share|improve this answer
    
Why are you using a Gaussian window? As is correctly stated in other answers, using no window will provide the best ability to resolve nearby frequency components. Regarding method 1, the algorithm typically used for calculating what amounts to a DFT bin at an arbitrary location (i.e. a fractional bin index) is the Goertzel algorithm. For a small number of frequencies (as in this example, 9), using the Goertzel algorithm to estimate the frequency content at each frequency of interest can be easier than doing a full DFT. –  Jason R Jan 24 '13 at 14:06
    
The key is the ability to estimate frequency after you take the natural log. Please see mgasior.web.cern.ch/mgasior/pap/biw2004.pdf for the max error and gain characteristics of different windows under Gaussian interpolation. After taking the log, a Gaussian window is vastly superior to other windows for frequency estimation (in another paper, the gain of the rectangular after the log is shown to be about 3, compared to thousands for the Gaussian ones). –  Kevin McGee Jan 26 '13 at 6:45
    
As for the Goertzel, a DFT computed for a single bin is likely to be more efficient than Goertzel when computed on a modern processor which is optimized for FIR-like operations, as noted in : groups.google.com/group/comp.dsp/browse_thread/thread/… A full N-point DFT would make no sense, since that would be less efficient than an FFT (or nine single bin DFTs). –  Kevin McGee Jan 26 '13 at 6:52
add comment

As is well-described in this answer, window functions are characterized by their mainlobe width and maximum sidelobe height. There is an inherent tradeoff between the two (decreasing one will increase the other). For the case where you are trying to discriminate between two nearby frequencies at similar power levels, you want a window with as small of a main lobe as possible. That window is the rectangular window (i.e. none at all). However, it does have rather poor sidelobe performance, which can hurt you if you're trying to look for signals that are separated significantly in frequency but are at different power levels.

If you need even better frequency resolution than what the unwindowed approach will give you, the correct answer is to observe the signal longer and then analyze it using a longer DFT. Note that there is no free lunch; you can't just zero-pad the existing signal and then take the DFT (which will just interpolate between the bins you already have). You'll need to actually collect a longer block of samples to start with.

share|improve this answer
add comment

A rectangular window (no window) will give you the deepest notch in the immediately adjacent bins, however it is not very robust against "spectral leakage" from noise. Any other window will cause significant spreading of any carrier into the immediately adjacent bins, so you won't get a 100X notch in those adjacent bins. For looking 2 or 3 bins away, you could specify a window with a low stop-band in that region and beyond.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.