FFT - 计算频率间的精确值

Question

我正在使用一个我在网上找到的很好的FFT库来尝试编写一个音高检测程序。目前为止，我已经成功地让库对包含几个正弦波（其中一个频率为440Hz）的测试音频信号进行FFT计算（我将采样大小设置为16384，采样率设置为44100Hz）。

FFT输出如下：

433.356Hz - Real: 590.644 - Imag: -27.9856 - MAG: 16529.5
436.047Hz - Real: 683.921 - Imag: 51.2798 - MAG: 35071.4
438.739Hz - Real: 4615.24 - Imag: 1170.8 - MAG: 5.40352e+006
441.431Hz - Real: -3861.97 - Imag: 2111.13 - MAG: 8.15315e+006
444.122Hz - Real: -653.75 - Imag: 341.107 - MAG: 222999
446.814Hz - Real: -564.629 - Imag: 186.592 - MAG: 105355

正如您所见，441.431Hz和438.739Hz的频率均显示出同等高的幅度输出（“MAG:”后面的最右侧数字），因此很明显目标频率440Hz介于两者之间。增加分辨率可能是一种逼近的方法，但这会增加计算时间。

我该如何计算介于两个频率区间之间的确切频率？

更新：

我尝试了DSPGuru网站上讨论的{{link1：Barry Quinn的“第二估计器”}}，并取得了很好的结果。以下是440Hz方波的结果 - 现在只有0.003Hz的偏差！

这是我使用的代码。我只是改编了我找到的这个例子，它是为Swift准备的。谢谢大家提供宝贵的意见，这是一次很棒的学习之旅 :)

- morezoom

---

你所寻找的方法被称为DESA（离散能量分离算法）。 - DaBler

谢谢您的提示，我会研究一下。我还发现了这篇文章，其中展示了几种解决方案... - morezoom

我预计在计算上最具节俭性和可靠性的方法是过采样FFT（即，对输入进行零填充，特别是到某个2的幂）。但这里有一些其他技术可以参考：https://gist.github.com/fasiha/957035272009eb1c9eb370936a6af2eb - Ahmed Fasih

我曾想过是否可以利用谐波内容来帮助估算。我一定会更深入地探索这些谐波产品谱技术。谢谢！ - morezoom

Answer 1

Sinc插值可以准确地插值（或重建）FFT结果间的频谱。 零填充FFT将产生类似的插值频谱。 您可以使用高质量的插值器（例如窗口化的Sinc核）并采用逐步逼近法来估计实际的频谱峰值，以达到信噪比所允许的任何分辨率。 除非在插值核中包括频谱的共轭图像的影响，否则该重建可能无法在DC或Fs/2 FFT结果间正常工作。

有关时域重建的详细信息，请参见https://ccrma.stanford.edu/~jos/Interpolation/Ideal_Bandlimited_Sinc_Interpolation.html和https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula，但是相同的插值方法也适用于频率或时间上的带限或时间限制信号。

如果你需要一个精度较低、计算负担较小的估算器，那么抛物线插值（以及其他类似的曲线拟合估算器）可能会起作用。详见：https://www.dsprelated.com/freebooks/sasp/Quadratic_Interpolation_Spectral_Peaks.html 和 https://mgasior.web.cern.ch/mgasior/pap/FFT_resol_note.pdf 关于抛物线的细节，以及 http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.555.2873&rep=rep1&type=pdf 关于其他曲线拟合峰值估算器。

- hotpaw2

---

感谢您为我解释了所有这些不同的插值技术。我确实尝试了二次和其他更简单的插值方法，遵循此页面，并取得了改进，但仍然存在最多0.5Hz的误差。所以我想我需要更加认真地阅读有关Sinc插值和窗口化Sinc滤波器的内容 :) - morezoom

Answer 2

要计算“真实”频率，我曾经使用过抛物线拟合算法。对于我的用例来说，它非常有效。

以下是我用来找到基本频率的方法：

• 计算DFT（WOLA）。
• 找到DFT bin中的峰值。
• 找到谐波产品谱。虽然不是最可靠和精确的方法，但这是一种非常简单的找到基本频率候选者的方法。
• 基于峰值和HPS，使用抛物线拟合算法找到基本音高频率（如果需要，还可以找到振幅）。

例如，HPS表示基本（最强）音高集中在您的DFT的bin x中；如果bin x属于peak y，则抛物线拟合频率将从peak y中取出，并且那就是您正在寻找的音高。

如果您不是在寻找基本音高，而是要在任何一个频率区间内精确测量频率，请对该区间应用二次曲线拟合。

以下是一些代码供您参考：

struct Peak
{
    float         freq     ; // Peak frequency calculated by parabola fit algorithm. 
    float         amplitude; // True amplitude.   
    float         strength ; // Peak strength when compared to neighbouring bins.         
    uint16_t      startPos ; // Peak starting position (DFT bin).
    uint16_t      maxPos   ; // Peak location (DFT bin).
    uint16_t      stopPos  ; // Peak stop position (DFT bin).
}; 

void calculateTrueFrequency( Peak & peak, float const bins, uint32_t const fs, DFT_Magnitudes mags )
{
    // Parabola fit:
    float a = mags[ peak.maxPos - 1 ];
    float b = mags[ peak.maxPos     ];
    float c = mags[ peak.maxPos + 1 ];

    float p   = 0.5f * ( a - c ) / ( a - 2.0f * b + c );
    float bin = convert<float>( peak.maxPos ) + p;

    peak.freq      = convert<float>( fs ) * bin / bins / 2;
    peak.amplitude = b - 0.25f + ( a - c ) * p;
}

- Danijel

---

谢谢您详细的回答！我猜这也被称为二次方法？振幅公式也会很有用。对于音频信号，最大的幅度往往似乎是基频，但我会看看使用HPS是否可以更好地帮助。 - morezoom

关于“最大振幅”，它远不止于此。我肯定会实现HPS，因为它并不复杂。祝好运。 - Danijel