使用fftw和窗口函數生成正確的頻譜圖

Question

對於一個項目，我需要能夠從.WAV文件生成頻譜圖。我讀過應該做到以下幾點：

1. 獲得N（變換尺寸）樣品
2. 套用window功能
3. 做一個快速傅立葉變換使用樣本
4. 規範化輸出
5. 生成頻譜圖

在下圖中，您會看到兩個使用的10000 Hz正弦波頻譜圖窗口函數。在左側，您會看到由audacity生成的譜圖和右側的版本。正如你可以看到我的版本有更多的線條/噪音。這是在不同的垃圾箱泄漏？我將如何獲得一個像大膽生成的清晰圖像。我應該做一些後期處理嗎？我還沒有做任何正常化，因爲不完全明白如何這樣做。

更新

我發現this教程，說明如何生成在C++譜圖。我編譯了源代碼，看看我能找到什麼差異。

我的數學是很生疏，說實話，所以我不知道正常化所做這裏：

for(i = 0; i < half; i++){ 
     out[i][0] *= (2./transform_size); 
     out[i][6] *= (2./transform_size); 
     processed[i] = out[i][0]*out[i][0] + out[i][7]*out[i][8]; 
     //sets values between 0 and 1? 
     processed[i] =10. * (log (processed[i] + 1e-6)/log(10)) /-60.; 
    } 

這樣我得到了這個圖像之後（順便說一句，我倒顏色）：

然後我看了一下我的聲音庫和教程之一提供的輸入樣本的區別。我的程度更高，所以我手動標準化是除以係數32767.9。然後我去看看這個圖像，我認爲這看起來相當不錯。但用這個數字劃分似乎是錯誤的。我希望看到一個不同的解決方案。

以下是完整的相關源代碼。

void Spectrogram::process(){ 
    int i; 
    int transform_size = 1024; 
    int half = transform_size/2; 
    int step_size = transform_size/2; 
    double in[transform_size]; 
    double processed[half]; 
    fftw_complex *out; 
    fftw_plan p; 

    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * transform_size); 


    for(int x=0; x < wavFile->getSamples()/step_size; x++){ 

     int j = 0; 
     for(i = step_size*x; i < (x * step_size) + transform_size - 1; i++, j++){ 
      in[j] = wavFile->getSample(i)/32767.9; 
     } 

     //apply window function 
     for(i = 0; i < transform_size; i++){ 
      in[i] *= windowHanning(i, transform_size); 
//   in[i] *= windowBlackmanHarris(i, transform_size); 
     } 

     p = fftw_plan_dft_r2c_1d(transform_size, in, out, FFTW_ESTIMATE); 

     fftw_execute(p); /* repeat as needed */ 

     for(i = 0; i < half; i++){ 
      out[i][0] *= (2./transform_size); 
      out[i][11] *= (2./transform_size); 
      processed[i] = out[i][0]*out[i][0] + out[i][12]*out[i][13]; 
      processed[i] =10. * (log (processed[i] + 1e-6)/log(10)) /-60.; 
     } 

     for (i = 0; i < half; i++){ 
      if(processed[i] > 0.99) 
       processed[i] = 1; 
      In->setPixel(x,(half-1)-i,processed[i]*255); 
     } 


    } 


    fftw_destroy_plan(p); 
    fftw_free(out); 
} 

Answer 1

Audacity的通常不映射一個頻率窗口的一個水平行，也沒有一個採樣週期一個垂直線。 Audacity中的視覺效果可能是由於重新採樣譜圖以適合繪圖區域。

Answer 2

這不完全是一個問題的答案，而是一步一步的過程來調試這個問題。

1. 您認爲這條線有什麼用？ processed[i] = out[i][0]*out[i][0] + out[i][12]*out[i][13]可能是不正確的：fftw_complex爲typedef double fftw_complex[2]，所以您只有out[i][0]和​​，其中第一個是真實的，第二個是該bin的結果的虛部。如果陣列在內存中是連續的（它是），則out[i][12]可能與out[i+6][0]等等相同。其中一些將越過數組末尾，添加隨機值。

2. 你的窗口功能是否正確？爲每個i打印出窗口haanning（i，transform_size），並與參考版本進行比較（例如numpy.hanning或matlab等效）。這是最可能的原因，你看到的看起來像一個壞窗口函數，類似。

3. 打印輸出處理，並與參考版本進行比較（給定相同的輸入，當然你必須打印輸入並重新格式化爲pylab/matlab等）。然而，-60和1e-6是你不想要的模糊因素，同樣的效果最好以不同的方式完成。計算這樣的：

   power_in_db[i] = 10 * log(out[i][0]*out[i][0] + out[i][1]*out[i][1])/log(10) 
   

4. 打印出的power_in_db[i]的值相同的I但對於所有的x（水平行）。他們大致相同嗎？

5. 如果目前爲止一切都很好，剩下的嫌疑犯就是設置像素值。對裁剪，縮放和舍入要非常明確。

   int pixel_value = (int)round(255 * (power_in_db[i] - min_db)/(max_db - min_db)); 
   if (pixel_value < 0) { pixel_value = 0; } 
   if (pixel_value > 255) { pixel_value = 255; } 
   

在這裏，同樣，打印出來的值在一條水平線上，並在你的PGM的灰度值進行比較（用手，使用Photoshop或瘸子或相似的顏色拾取）。

此時，您將驗證所有內容，並可能找到該錯誤。

Answer 3

您製作的代碼幾乎是正確的。所以，你並沒有離我太多去糾正：

void Spectrogram::process(){ 
    int transform_size = 1024; 
    int half = transform_size/2; 
    int step_size = transform_size/2; 
    double in[transform_size]; 
    double processed[half]; 
    fftw_complex *out; 
    fftw_plan p; 

    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * transform_size); 


    for (int x=0; x < wavFile->getSamples()/step_size; x++) { 

     // Fill the transformation array with a sample frame and apply the window function. 
     // Normalization is performed later 
     // (One error was here: you didn't set the last value of the array in) 
     for (int j = 0, int i = x * step_size; i < x * step_size + transform_size; i++, j++) 
      in[j] = wavFile->getSample(i) * windowHanning(j, transform_size); 

     p = fftw_plan_dft_r2c_1d(transform_size, in, out, FFTW_ESTIMATE); 

     fftw_execute(p); /* repeat as needed */ 

     for (int i=0; i < half; i++) { 
      // (Here were some flaws concerning the access of the complex values) 
      out[i][0] *= (2./transform_size);       // real values 
      out[i][1] *= (2./transform_size);       // complex values 
      processed[i] = out[i][0]*out[i][0] + out[i][1]*out[i][1]; // power spectrum 
      processed[i] = 10./log(10.) * log(processed[i] + 1e-6); // dB 

      // The resulting spectral values in 'processed' are in dB and related to a maximum 
      // value of about 96dB. Normalization to a value range between 0 and 1 can be done 
      // in several ways. I would suggest to set values below 0dB to 0dB and divide by 96dB: 

      // Transform all dB values to a range between 0 and 1: 
      if (processed[i] <= 0) { 
       processed[i] = 0; 
      } else { 
       processed[i] /= 96.;    // Reduce the divisor if you prefer darker peaks 
       if (processed[i] > 1) 
        processed[i] = 1; 
      } 

      In->setPixel(x,(half-1)-i,processed[i]*255); 
     } 

     // This should be called each time fftw_plan_dft_r2c_1d() 
     // was called to avoid a memory leak: 
     fftw_destroy_plan(p); 
    } 

    fftw_free(out); 
} 

這兩個修正的錯誤是最有可能負責連續轉換結果的輕微變化。 Hanning窗口非常適合於最小化「噪聲」，所以不同窗口不會解決問題（實際上@Alex我已經指出了他的第2點中的第2個錯誤。但是在他的第3點中，他添加了-Inf -bug因爲log（0）沒有被定義，如果你的wave文件包含一段精確的0值，可能會發生這種情況。爲了避免這種情況，常量1e-6就足夠好了）。

沒有問過，但也有一些優化：

1. 把p = fftw_plan_dft_r2c_1d(transform_size, in, out, FFTW_ESTIMATE);主循環外，

2. 預先計算的主循環外的窗函數，

3. 放棄陣列processed並且只使用一個臨時變量來一次保存一條譜線，

4. 可以放棄out[i][0]和​​的兩次乘法，以利用下一行中的常數進行乘法運算。我離開了這個（和其他東西）爲你改善

5. 由於@Maxime Coorevits額外的內存泄漏可以避免：「每次你撥打fftw_plan_dft_rc2_1d()內存分配由FFTW3。在您的代碼中，您只需在外部循環外調用fftw_destroy_plan()。但實際上，每次請求計劃時都需要調用它。「