11

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

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

在下图中，您可以看到两个使用汉宁窗函数的 10000 Hz 正弦波的频谱图。在左边你会看到一个由audacity生成的频谱图，在右边是我的版本。如您所见，我的版本有更多的线条/噪音。这是在不同的垃圾箱中泄漏吗？我将如何获得像大胆生成的清晰图像。我应该做一些后期处理吗？我还没有做任何规范化，因为不完全理解如何去做。

在此处输入图像描述

更新

我发现本教程解释了如何在 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);
}
4

3 回答 3

5

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

  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]等等。其中一些将超出数组的末尾，添加随机值。

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

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

于 2014-01-31T14:39:42.737 回答
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 点中指出了第二个错误。但在他的第 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()在外循环之外调用。但实际上，每次请求时都需要调用它计划。”

于 2014-02-01T12:14:43.043 回答
2

Audacity typically doesn't map one frequency bin to one horizontal line, nor one sample period to one vertical line. The visual effect in Audacity may be due to resampling of the spectrogram picture in order to fit the drawing area.

于 2014-01-27T21:11:30.260 回答