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

Question

对于项目,我需要能够从.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);
}

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

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