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

Question

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

1. 获取 N 个（变换大小）样本
2. 应用window 函数
3. 使用样本进行快速傅里叶变换
4. 标准化输出
5. 生成频谱图

在下图中，您可以看到两个使用hanning 窗口函数的 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

这并不完全是错误的答案，而是一个逐步调试的过程。

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

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 或类似工具中的颜色选择器）。

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

Answer 2

您生成的代码几乎是正确的。所以，你并没有给我留下太多要纠正的地方：

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 点中指出了第二个错误。但在他的第 3 点中。他添加了 -Inf -bug as log(0) 未定义，如果您的波形文件包含一系列精确的 0 值，则可能发生这种情况。为避免这种情况，常量 1e-6 就足够了。

没有问，但有一些优化：

1. 将p = fftw_plan_dft_r2c_1d(transform_size, in, out, FFTW_ESTIMATE);放在主循环之外，

2. 在主循环外预计算窗口函数，

3. 放弃数组processed，只使用一个临时变量一次保存一条谱线，

4. out[i][0] 和 out[i][1] 的两个乘法可以被放弃，取而代之的是下一行中的一个常数乘法。我把这个（和其他东西）留给你改进

5. 感谢@Maxime Coorevits，另外还可以避免内存泄漏：“每次调用fftw_plan_dft_rc2_1d() 时，FFTW3 都会分配内存。在您的代码中，您只能在外循环外调用fftw_destroy_plan()。但在事实上，您每次请求计划时都需要调用它。”

Answer 3

Audacity 通常不会将一个频率箱映射到一条水平线，也不会将一个采样周期映射到一条垂直线。 Audacity 中的视觉效果可能是由于重新采样频谱图以适应绘图区域。