Real FFT output

Question

I have implemented fft into at32ucb series ucontroller using kiss fft library and currently struggling with the output of the fft. My intention is to analyse sound coming from piezo speaker. Currently, the frequency of the sounder is 420Hz which I successfully got from the fft output (cross checked with an oscilloscope). However, the output frequency is just half of expected if I put function generator waveform into the system. I suspect its the frequency bin calculation formula which I got wrong; currently using, fft_peak_magnitude_index*sampling frequency / fft_size. My input is real and doing real fft. (output samples = N/2) And also doing iir filtering and windowing before fft. Any suggestion would be a great help!

            // IIR filter calculation, n = 256 fft points       
        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] - 510; //eliminate dc offset

            // hamming window
            hamming[ctr] = (0.54-((0.46) * cos(2*M_PI*ctr/n)));
            window[ctr] = hamming[ctr]*y1[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);

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

                if(fft_mag[ctr] > peak)
                {
                    peak = fft_mag[ctr];
                    freq_bin = ctr;
                }

            frequency = (freq_bin*(10989/n)); // 10989 is the sampling freq
                //************************************
                //Usart write
                char filtResult[10];
                //sprintf(filtResult, "%04d %04d %04d\n", (int)peak, (int)freq_bin, (int)frequency);
                sprintf(filtResult, "%04d  %04d  %04d\n", (int)x[ctr], (int)fft_mag[ctr], (int)frequency);
                char c;
                char *ptr = &filtResult[0];
                do
                {
                    c = *ptr;
                    ptr++;
                    usart_bw_write_char(&AVR32_USART2, (int)c);
                    // sendByte(c);

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

Answer 1

The main problem is likely to be how you declared fft_input. Based on your previous question, you are allocating fft_input as an array of kiss_fft_cpx. The function kiss_fftr on the other hand expect an array of scalar. By casting the input array into a kiss_fft_scalar with:

kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);

KissFFT essentially sees an array of real-valued data which contains zeros every second sample (what you filled in as imaginary parts). This is effectively an upsampled version (although without interpolation) of your original signal, i.e. a signal with effectively twice the sampling rate (which is not accounted for in your freq_bin to frequency conversion). To fix this, I suggest you pack your data into a kiss_fft_scalar array:

kiss_fft_scalar fft_input[n];
...
for (ctr=0; ctr<n; ctr++)
{       
    ...
    fft_input[ctr] = window[ctr];
    ...
}
kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, fft_input, fft_output);

Note also that while looking for the peak magnitude, you probably are only interested in the final largest peak, instead of the running maximum. As such, you could limit the loop to only computing the peak (using freq_bin instead of ctr as an array index in the following sprintf statements if needed):

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

  if(fft_mag[ctr] > peak)
  {
    peak = fft_mag[ctr];
    freq_bin = ctr;
  }
} // close the loop here before computing "frequency"

Finally, when computing the frequency associated with the bin with the largest magnitude, you need the ensure the computation is done using floating point arithmetic. If as I suspect n is an integer, your formula would be performing the 10989/n factor using integer arithmetic resulting in truncation. This can be simply remedied with:

frequency = (freq_bin*(10989.0/n)); // 10989 is the sampling freq