深入浅出Goertzel算法

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

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

一、从DFT说起

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

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

X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi k n}{N}} = \sum_{n=0}^{N-1} x[n] \cdot W_N^{kn}

其中 $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 引入一个"恒等式"

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

W_N^{-kN} = e^{j \frac{2\pi k N}{N}} = e^{j 2\pi k} = 1

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

X[k] = W_N^{-kN} \cdot \sum_{n=0}^{N-1} x[n] \cdot W_N^{kn} = \sum_{n=0}^{N-1} x[n] \cdot W_N^{-k(N-n)}

2.2 构造卷积/滤波的形式

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

y[n] = \sum_{m=0}^{N-1} x[m] \cdot W_N^{-k(n-m)}

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

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

X[k] = y[n] \Big|_{n=N} = y[N]

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

2.3 推导递推关系

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

H(z) = \frac{1}{1 - W_N^{-k} z^{-1}}

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

H(z) = \frac{1 - W_N^{k} z^{-1}}{(1 - W_N^{-k} z^{-1})(1 - W_N^{k} z^{-1})} = \frac{1 - W_N^{k} z^{-1}}{1 - 2\cos\left(\frac{2\pi k}{N}\right) z^{-1} + z^{-2}}

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

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

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

s[n] = x[n] + 2\cos\left(\frac{2\pi k}{N}\right) \cdot s[n-1] - s[n-2]

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

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

y[n] = s[n] - W_N^{k} \cdot s[n-1]

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

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

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$ ）。这时可以进一步简化。

利用结算公式：

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

计算其模的平方：

|X[k]|^2 = s[N]^2 + s[N-1]^2 - 2\cos\left(\frac{2\pi k}{N}\right) \cdot s[N] \cdot s[N-1]

发现了吗？以上全程零复数运算！ 递推阶段本身就是纯实数的，结算阶段也变成了纯实数运算。这意味着：
计算量减少：不需要维护实部和虚部，乘法和加法次数都更少
易于定点化：全实数的运算改写为定点数版本非常直接，这对于某些没有FPU的嵌入式处理器来说十分友好

四、计算复杂度对比

方法	计算目标	复杂度
直接DFT（单个 $k$ ）	$X[k]$	$O(N)$ 次复数乘法
FFT（全部 $N$ 个）	$X[0] \sim X[N-1]$	$O(N \log N)$ 次复数运算
Goertzel（单个 $k$ ）	$X[k]$	$O(N)$ 次实数乘法 + 1次复数乘法

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

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

M < \frac{N \log_2 N}{N} = \log_2 N

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

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

5.1 $k$ 值的确定

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

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

k = \frac{f_{\text{target}}}{f_s} \cdot N

5.2 $N$ 的选择

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

\Delta f = \frac{f_s}{N}

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

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

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

8.1 窗函数

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

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

8.2 滑动Goertzel

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