Реальный выход FFT

Question

Я реализовал fft в ucontroller at32ucb, используя библиотеку поцелуя fft и в настоящее время борюсь с выходом fft. Мое намерение состоит в том, чтобы анализировать звук, исходящий от пьезодинамика. В настоящее время частота эхолота составляет 420 Гц, которую я успешно получил от выходного сигнала fft (перекрестная проверка с помощью осциллографа). Тем не менее, выходная частота составляет всего половину ожидаемого, если я поставлю в генераторе сигнал функции генератора. Я подозреваю, что его формула вычисления бинарной частоты, которую я получил неправильно; в настоящее время используется fft_peak_magnitude_index * частота дискретизации/fft_size. Мои данные реальны и делают реальные fft. (выходные отсчеты = N/2) А также делает фильтрацию и оконтирование перед fft. Любое предложение было бы большой помощью!

  // 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

Основная проблема, вероятно, заключается в том, как вы заявили fft_input. На основе your previous question вы выделяете fft_input в виде массива kiss_fft_cpx. С другой стороны, функция kiss_fftr ожидает массив скаляров. Бросая массив ввода в kiss_fft_scalar с:

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

KissFFT по существу видит массив вещественных данных, которая содержит нули каждый второй образец (то, что вы заполнили как мнимые части). Это эффективная версия с улучшенной дискретизацией (хотя и без интерполяции) исходного сигнала, т. Е. Сигнал с эффективностью в два раза большей частоты дискретизации (которая не учитывается в вашем преобразовании-frequency). Чтобы исправить это, я предлагаю вам упаковать свои данные в kiss_fft_scalar массив:

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); 

Заметим также, что при поиске пика величины, вы, вероятно, заинтересованы только в конечном величине пика, вместо бега максимума. Таким образом, можно ограничить цикл только вычисления пики (с использованием freq_bin вместо ctr в качестве индекса массива в следующих sprintf отчетности, если это необходимо):

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" 

Наконец, при вычислении частоты, связанную с бункером с наибольшая величина, вам нужно обеспечить, чтобы вычисления выполнялись с использованием арифметики с плавающей запятой. Если, как я подозреваю, n является целым числом, ваша формула будет выполнять коэффициент 10989/n, используя целочисленную арифметику, приводящую к усечению. Это можно просто исправить:

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