真正的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