KISS FFT exit with or without window
I am currently trying to implement fft in avr32 microcontrollers to handle signals using kiss fft. And a strange problem with my exit. Basically, I am skipping ADC samples (testing with a function generator) in fft (real input, size 256 n) and the resulting result makes sense to me. However, if I applied a Hamming window to the ADC samples and then fed them to the FFT, the peak frequency bit is in error (and different from the previous windowless result). The ADC samples are DC offset, so I removed the offset, but still it doesn't work with windows.
Below are the first few output values via rs485. The first column is fft output without window and second column is output with window. From column 1, the peak is in row 6 (6 x fs (10.5 kHz) / 0.5N) gave me the correct input frequency result where column 2 has the maximum amplitude in row 2 (except dc bin), which for me does not have meaning, 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);
source to share
Frequency domain spreading profiles
In the frequency domain, rectangular and chrome windows look like this:
As you know, multiplying in the time domain by the window corresponds to convolution in the frequency domain, which significantly expands the signal energy over several frequency bins in the so-called spectral leakage . For certain windows you have selected (pictured above in the frequency domain), you may notice that the Hamming window has significantly less energy outside the main lobe, but this main blade is slightly wider than the rectangular window.
As a result, a burst of DC energy ends up spreading past bin 0 and into bin 1 when using the Hamming window. It is not that you have a strong peak in box 1. In fact, if you generate the data you provided, you should see that the surge you saw at index 6 actually still exists at the same frequency , and without using the Hamming window:
If you want to get rid of DC spike and leakage in surrounding containers, either remove bias in your data (mainly by applying a notch filter), or you will have to ignore a few lower frequency boxes when looking for your "first big burst".
Filtering problem
Finally, note that there is also a problem with how the IIR filter is implemented, causing the arrays x
and y
will be indexed out of bounds, if ctr==0
and ctr==1
(unless you have made some special conditions and x
and y
are actually pointers with an offset from the beginning of the allocated buffers ). This can affect results with or without a window. If you are only filtering one block of data, the general assumption is that the previous patterns were zeros. In this case, you can avoid non-anchor indexing:
// 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 samples n
, you will need to remember the last few samples from the previous block. This can 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 the offset in that buffer. Indexes 0 and 1 represent passed samples, while the rest of the buffer from index 2 is filled with the current block of input. This results in 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 need to update the "past samples" at index 0 and 1, with the last samples of the current block 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];
source to share