c++ - 使用 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. 针对相同的 i 但针对所有 x(水平线)打印出 power_in_db[i] 的值。它们大致相同吗？

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 或类似软件中的颜色选择器)。

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