真正的FFT输出

Question

我已经使用kiss fft库实现了fft到at32ucb系列ucontroller，并且目前正在与fft的输出挣扎。 我的意图是分析来自压电扬声器的声音。 目前，测深仪的频率是420Hz，我从fft输出成功得到（用示波器交叉检查）。但是，如果将函数发生器波形放入系统，输出频率只是预期的一半。 我怀疑它的频率箱计算公式，我错了;目前使用的是fft_peak_magnitude_index *采样频率/ fft_size。 我的输入是真实的，并做到真正的fft。 （输出采样= N/2） 并且还在fft之前进行iir滤波和窗口化。 任何建议将是一个很大的帮助！

  // 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_fftr另一方面期望一个标量数组。通过铸造输入数组成kiss_fft_scalar与：

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

KissFFT基本上看到实值数据的阵列，其包含零每第二个样品（你在作为虚数部分填充的）。这实际上是您原始信号的上采样版本（尽管没有内插），即信号采样率实际上是两倍（这在您的freq_bin到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); 

还要注意的是，虽然寻找峰值幅度，你可能只是在最后的最大峰感兴趣，而不是运行最大值。因此，可以将循环限制（如果需要的话，如以下sprintf语句数组索引使用freq_bin的ctr代替）只计算峰：

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