使用 fftw 和窗函数生成正确的频谱图

对于一个项目，我需要能够从 .WAV 文件生成频谱图。我读过以下应该做的事情：

1. 获取N（变换大小）个样本
2. Apply a window http://en.wikipedia.org/wiki/Window_function功能
3. 使用样本进行快速傅里叶变换
4. 标准化输出
5. 生成频谱图

在下图中，您可以看到两个 10000 Hz 正弦波的频谱图，均使用hanning http://en.wikipedia.org/wiki/Window_function#Hann_.28Hanning.29_window窗函数。在左侧您可以看到由以下命令生成的频谱图audacity http://audacity.sourceforge.net/右边是我的版本。正如你所看到的，我的版本有更多的线条/噪音。是否泄漏在不同的垃圾箱中？我怎样才能获得像大胆产生的清晰图像。我应该做一些后期处理吗？我还没有进行任何标准化，因为不完全了解如何做到这一点。

update

I found this http://www.gemmantics.com/spectrograms-in-c-with-libav-ffmpeg-libpng-and-sphinx/教程解释如何用 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);
}

这并不完全是问题所在的答案，而是调试此问题的逐步过程。

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] and out[i][1]，其中第一个是该 bin 结果的实部，第二个是虚部。如果数组在内存中是连续的（确实如此），那么out[i][12]可能与out[i+6][0]等等。其中一些will越过数组末尾，添加随机值。

2. 你的窗口函数正确吗？打印出每个 i 的 windowHanning(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 或 gimp 或类似工具中的颜色选择器）。

此时，您将从头到尾验证所有内容，并可能发现错误。