深入浅出Goertzel算法

在数字信号处理（DSP）中，当需要分析信号的频谱的时候，最常用也最为大众熟悉的应该就是快速傅里叶变换（FFT）。但是在某些特定场景，我们只关心某个或某几个特定频率的时候，计算完整的FFT很划不来，尤其是在性能受限的嵌入式设备上或要求高实时性的场景中。

这时候，就该Goertzel算法登场了，它不需要扫描整个频谱，它只关心其中的一个（或几个）频率分量，能够高效地计算出其能量。

一、从DFT说起

要理解Goertzel算法，我们首先回顾一下离散傅里叶变换（DFT）的定义。

对于长度为$N$的离散信号$x[n]$，其第$k$个频率分量的DFT定义为：

其中$W_N = e^{-j\frac{2\pi}{N}}$是旋转因子（twiddle factor）。

如果我们只需要计算某一个$k$对应的$x[k]$，直接按定义计算就需要$N$次复数乘法和$N-1$次复数加法。而完整的FFT虽然高效（$O(N \log N)$），但它还是一口气算出了所有$N$个频率分量。如果我只关心某一个频率，FFT就做了大量无“无用功”，有没有办法只算我要的那个呢？

有的有的，答案就是Goertzel算法。

二、Goertzel算法的推导

Goertzel算法的核心思想非常巧妙：把DFT的直接计算转化为一个二阶IIR滤波器（递归结构）。

2.1 引入一个"恒等式"

注意到旋转因子有一个周期性质：

利用这个性质，我们可以在DFT表达式上乘以一个值为1的因子，不改变结果：

2.2 构造卷积/滤波的形式

我们定义一个中间序列$y[n]$，将其看作信号$x[n]$通过一个滤波器后的输出：

更准确地说，我们定义滤波器的冲激响应为$h[n] = W_N^{-kn} \cdot u[n]$（其中$u[n]$是单位阶跃函数），那么$y[n]$就是$x[n]$和$h[n]$的卷积。

对比上面的式子，不难发现：

也就是说，DFT的第$k$个分量，等于信号通过特定滤波器后在$n=N$时刻的输出值。

2.3 推导递推关系

$h[n] = W_N^{-kn} \cdot u[n]$对应的系统函数为：

这是一个一阶IIR滤波器，但它涉及复数系数，计算起来并不方便。Goertzel的精妙之处在于，将分子分母同时乘以$(1 - W_N^{k} z^{-1})$，得到：

注意分母展开利用了$W_N^{-k} + W_N^{k} = 2\cos\left(\frac{2\pi k}{N}\right)$以及$W_N^{-k} \cdot W_N^{k} = 1$。

这样我们就把系统拆成了两个部分：

第一部分（二阶递推，全实数运算）：

初始条件：$s[-1] = 0$，$s[-2] = 0$

第二部分（最终输出，只需计算一次）：

而我们需要的结果是$X[k] = y[N]$，即：

2.4 小结：算法流程

整个Goertzel算法可以归纳为两个阶段：

1. 递推阶段（执行$N$次迭代）：

令$\text{coeff} = 2\cos\left(\frac{2\pi k}{N}\right)$

对$n = 0, 1, 2, \dots, N-1$：

$s[n] = x[n] + \text{coeff} \cdot s[n-1] - s[n-2]$

每次迭代只需要 1次实数乘法 + 2次实数加法

2. 结算阶段（只执行 1 次）：

$X[k] = s[N] - W_N^{k} \cdot s[N-1]$

这一步需要 1 次复数乘法

三、只需要功率谱？更简单！

在很多实际应用中（比如DTMF检测），我们并不关心$X[k]$的相位，只关心功率（即$|X[k]|^2$）。这时可以进一步简化。

利用结算公式：

计算其模的平方：

发现了吗？以上全程零复数运算！ 递推阶段本身就是纯实数的，结算阶段也变成了纯实数运算。这意味着：

计算量直接减半：不需要维护实部和虚部，乘法和加法次数都更少

易于定点化：全实数的运算改写为定点数版本非常直接，这对于没有FPU的嵌入式处理器来说十分友好

代码实现更简单：不需要复数数据结构，用几个普通变量就能搞定

四、计算复杂度对比

单独看一个频率，Goertzel和直接DFT都是$O(N)$，但Goertzel的常数因子更小（实数运算 vs 复数运算）。

什么时候Goertzel比FFT更划算？ 粗略估计，当你需要检测的频率数量$M$满足：

当然，这只是一个粗略的估计。实际中还需要考虑缓存命中率、指令流水线、是否有硬件FFT加速器等因素。

五、频率分辨率与参数选择

5.1 $k$值的确定

Goertzel算法中的$k$并不要求是整数。当$k$取非整数时，我们实际上是在计算“非标准频率点”上的DFT，这给了我们更大的灵活性。

给定目标频率$f_{\text{target}}$和采样率$f_s$，对应的$k$值为：

5.2 $N$的选择

$N$（采样点数）直接决定了频率分辨率：

$N$越大，频率分辨率越高，但计算延迟也越大。在实际工程中，需要在分辨率和实时性之间做权衡。

与FFT不同的是，Goertzel算法不要求$N$是2的幂次。这在实际应用中非常方便——你可以自由选择$N$ 来精确匹配你想要的频率分辨率。

六、实际使用中的注意事项

8.1 窗函数

和FFT一样，Goertzel也会受到频谱泄漏的影响。如果目标频率不恰好落在DFT的频率格点上（即$k$不是整数），能量会“泄漏”到相邻频率。

解决方案同样是在输入信号上施加窗函数（如Hanning窗、Hamming窗等）。

8.2 滑动Goertzel

在需要连续检测的场景中（如实时音频流），可以使用带滑动窗口的Goertzel算法，避免对每个新数据块从头开始计算。基本思路是：当窗口滑动一个样本时，增量更新$s[n]$的值，而不是重新计算整个递推过程。这在实时系统中可以进一步降低计算量。