KISS FFT output with or without windowing

Question

I am currently trying to implement fft into avr32 micro controllers for signal processing purposes using kiss fft. And having a strange problem with my output. Basically, I am passing ADC samples (testing with function generator) into fft (real input, 256 n size) and retrieved output makes sense to me. However, if I apply Hamming window to ADC samples and then pass these to FFT, the frequency bin of the peak magnitude is wrong (and different from previous result without windowing). ADC samples have DC offset, thus I eliminated the offset but still it doesnt work with windowed samples.

The below is the first few output values through rs485. First column is the fft output without window whereas second column is the output with window. From Column 1 the peak is at row 6 (6 x fs (10.5kHz) / 0.5N) gave me the correct input freq result where column 2 has a peak magnitude at row 2 (except dc bin) which does not make sense to me. Any suggestion would be helpful. Thanks in advance.

    488         260   //dc bin
    5            97
    5            41
    5            29  
    4            26
    10           35
    133          76
    33           28
    21           6
    17           3

kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_input[n];
kiss_fft_cpx fft_output[n];

for(ctr=0; ctr<n; ctr++)
{
    fft_input[ctr].r = zero;
    fft_input[ctr].i = zero;
    fft_output[ctr].r =zero;
    fft_output[ctr].i = zero;
}

// IIR filter calculation

for (ctr=0; ctr<n; ctr++)
{       
    // filter calculation
    y[ctr] = num_coef[0]*x[ctr];

    y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
    y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
    //y1[ctr] += y[ctr] - 500;
    // hamming window
    hamming[ctr] = (0.54-((0.46) * cos(2*PI*ctr/256)));
    window[ctr] = hamming[ctr]*y[ctr];

    fft_input[ctr].r = window[ctr];
    fft_input[ctr].i = 0;
    fft_output[ctr].r = 0;
    fft_output[ctr].i = 0;

}

kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);


for (ctr=0; ctr<n; ctr++)
{   
    fft_mag[ctr] = (sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);

    //Usart write
    char filtResult[10];
    sprintf(filtResult, "%04d %04d\n", (int)x[ctr], (int)fft_mag[ctr]);
    //sprintf(filtResult, "%04d %04d\n", (int)x[ctr], (int)window[ctr]);
    char c;
    char *ptr = &filtResult[0];
    do
    {   
        c = *ptr;
        ptr++;
        usart_bw_write_char(&AVR32_USART2, (int)c);
        // sendByte(c);

    } while (c != '\n');
}

kiss_fft_cleanup();
free(fftConfig);        

Stannary

Answer 1

Frequency domain output clarifications

In the frequency domain, the rectangular and Hamming windows look like:

As you may be aware, multiplication in the time domain by a window corresponds to a convolution in the frequency domain, which essentially spreads the energy of the signal over multiple frequency bins in what is often called spectral leakage. For the specific windows that you have chosen (depicted above in the frequency domain), you may notice that the Hamming window spreads much less energy outside the main lobe, but that main lobe is a little bit wider than that of the rectangular window.

As a result, the spike of DC energy ends up spreading past bin 0, and into bin 1 when using the Hamming window. It isn't so much that you have a strong peak in bin 1. In fact, if you plot the data you provided, you should see that the spike you were seeing at index 6 is in fact still there at the same frequency with and without the use of the Hamming window:

If you want to get rid of the DC spike and leakage in surrounding bins, either remove the bias in your data (essentially applying a notch filter), or you'll have to ignore a few more low-frequency bins when looking for your "first strong spike".

Filtering issue

Finally, note that there is also an issue with the way the IIR filter is implemented whereby the arrays x and y will be indexed out of bounds when ctr==0 and ctr==1 (unless you've made some special provision and x & y are in fact pointers with an offset from the start of the allocated buffers). This could affect the results for both with and without window. If you are only filtering a single block of data, a common assumption is that earlier samples were zeros. In that case, you can avoid the out-of-bound indexing with:

// filter calculation
y[ctr] = num_coef[0]*x[ctr];

if (ctr>=1)
{
  y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
}
if (ctr>=2)
{
  y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
}

If on the other hand you want to filter multiple blocks of n samples, you'll have to remember the last few samples from the previous block. This could be done by allocating buffers that are slightly larger than the block size:

x = malloc((n+2)*sizeof(kiss_fft_scalar));
y = malloc((n+2)*sizeof(kiss_fft_scalar));
// initialize "past samples" for the first block, assuming data was zero
x[0] = x[1] = 0;
y[0] = y[1] = 0;

Then you can use an offset in those buffer. Index 0 & 1 represent past samples, whereas the rest of the buffer from index 2 are filled with the current input data block. This leads to the following slightly modified filtering code:

// filter calculation
y[ctr+2] = num_coef[0]*x[ctr+2];

y[ctr+2] += (num_coef[1]*x[ctr+1]) - (den_coef[1]*y[ctr+1]);
y[ctr+2] += (num_coef[2]*x[ctr]) - (den_coef[2]*y[ctr]);

// hamming window
hamming[ctr] = (0.54-((0.46) * cos(2*PI*ctr/256)));
window[ctr] = hamming[ctr]*y[ctr+2];

Finally, at the end of each block you then have to update the "past samples" at index 0 & 1, with the last samples of the current block to be ready to process the next input block:

// remember last 2 samples of block
x[0] = x[n-2];
x[1] = x[n-1];
y[0] = y[n-2];
y[1] = y[n-1];

Subcortex

Wow, thank you for really nice clarification. Since I am already segmenting real time input from ADC, does this imply that rectangular windowing has been done to compensate Fourier Transform assumption (periodic signal)? And dont even have to bother about other windowing? And for your last comment about the filter, could you elaborate bit more? I am a newbie to programming world. thanks! – Stannary

Rectangular windowing is implicit when performing FFT on a block of data (from a possibly much longer sequence of data). However you may still want to use another window to adjust the tradeoff between frequency resolution (the width of the window's main lobe in the frequency domain) and sidelobe attenuation. – Subcortex