定义
戈泽尔算法是一种计算离散傅里叶变换的递归方法,它给出 DFT 中单个频点系数的估算过程[3]。对长度为 N 的离散信号 x(n),其第 k 个频点的 DFT 系数可写为对输入序列与旋转因子乘积的求和,戈泽尔算法把这一求和改写为滤波运算,从而逐点独立地求出 X(k)[4]。
与直接按定义计算 DFT 相比,戈泽尔算法所需的实数乘法次数约为后者的一半,实数加法次数大致相当,而三角函数求值次数则降为大约 N 分之一,这是它最突出的收益[5]。该算法可以按输入样点的到达顺序逐个处理,不必等到全部 N 个样点采集完毕才开始运算,这一点与需要整块数据的快速变换不同[2]。
原理
算法的出发点利用了旋转因子的周期性。由于在 k 为整数时旋转因子的 N 次幂等于 1,可以在 DFT 求和式中乘入该项而不改变结果,于是原式被整理成输入序列与一个无穷长指数序列的卷积形式,该卷积对应一个冲激响应为旋转因子幂次的系统[6]。这样,x(n) 可以看作激励,X(k) 则等于系统在 n = N 时刻的输出[6]。
对应的一阶系统函数为 $H_k(z)=\dfrac{1}{1-W_N^{-k}z^{-1}}$,它是一个系数为复数的一阶递归滤波器,每计算一个新输出需要 4 次实数乘法和 4 次实数加法[6]。为减少复数乘法,将系统函数的分子分母同乘一个因子,得到二阶形式 $H_k(z)=\dfrac{1-W_N^{k}z^{-1}}{1-2\cos(2\pi k/N)z^{-1}+z^{-2}}$,其极点为一对共轭极点、另有一个零点,反馈支路只含实数系数[5]。
二阶结构只需对反馈部分进行 N 次迭代,前馈部分在迭代结束后计算一次即可。中间变量满足递推式 $v_k[n]=x[n]+2\cos(2\pi k/N)\,v_k[n-1]-v_k[n-2]$,初值取零,最后由 $X(k)=v_k[N]-W_N^{k}\,v_k[N-1]$ 得到复数结果;若只需要幅度谱,则可由 v_k[N] 与 v_k[N-1] 直接算出模平方[2]。二阶形式只需存储一个余弦参数,总运算量约为 2N 次实数乘法和 2N 次实数加减法[1]。
flowchart LR
x["x(n) 输入样点"] --> A["反馈支路:2cos(2πk/N) 与延时单元"]
A --> B["中间变量 v_k(n) 迭代 N 次"]
B --> C["前馈支路:一次复数乘加"]
C --> D["X(k) 第 k 个 DFT 系数"]
发展历程
该算法在 1958 年由杰拉德·格策尔提出,他在当年发表的论文中给出了有限三角级数的求值方法,这一方法后来被引入数字信号处理领域并以其姓氏命名[3]。此后,数字信号处理的经典教材将其收入频域分析章节,并系统讨论了从一阶到二阶形式的推导与实现[6]。
随着电话网数字化的推进,双音多频拨号信号的识别成为一个典型需求。拨号时每个按键对应低频与高频各一个频率的叠加,需要同时判断若干个已知频率是否存在,这一场景与戈泽尔算法只针对给定频点计算的特点相契合,使其成为音频信号解码中的常用手段[1]。
在实现层面,研究者后来针对该算法的串行依赖问题提出了改进。由于每一步迭代都依赖前一步的结果,处理过程几乎没有并行空间,有专利文献提出让递推式跨越多个样点来拆分计算以获得并行度,但同时指出这会增加乘法次数,并可能带来数值稳定性问题[7]。
应用
电话系统是该算法最典型的应用场合。在双音多频信号检测中,需要同时评估 8 个不同频率,对每个输入样点分别以不同的 k 值运行算法即可完成判断;在北美的一些系统中,用于区分呼叫进展音调的频率数量可达 16 个[8]。除了拨号音识别,该算法还用于移相键控与频移键控调制解调器的信号处理[4]。
在测量与成像领域也有应用。例如在成像光谱仪的遥感数据处理中,需要对每个像素列求逆离散傅里叶变换的幅度以形成光谱图像,此类计算可以在流水线处理器上用戈泽尔算法完成,因为它能够按序接收样点而不是要求整块数据同时可用[2]。在电力信号处理中,该算法可用于电压或电流采样的谐波提取,因为很多场景只关心少量特定频点,而不需要完整的全谱信息[9]。
局限
该算法的效率优势只在频点数量很少时成立。若要获得完整频谱,其运算量仍与直接按定义计算 DFT 相当,呈 N 的平方量级,此时快速傅里叶变换的复杂度要低得多[5]。因此它并不适合替代快速变换用于全谱分析,而是在只需少数频点时作为补充手段[8]。
算法要求待检测的频率事先已知。目标频点由索引 k 确定,必须先限定在 0 到 N-1 的范围内选取,无法在不预先设定频率的情况下自动扫描任意未知成分,这一限制在把它用于音调检测时尤为明显[8]。
结构上的串行依赖也带来实现上的困难。递推式相邻项之间存在依赖,使算法在并行处理器上难以展开,最少需要 N 个时钟周期才能完成一次长度为 N 的运算,对 40 点信号意味着超过 40 个时钟周期[7]。此外,该滤波结构在定点实现中还存在定标问题,处理较长信号段时可能出现不稳定,因此通常只适合短信号段[7]。
参见
参考资料
- [科普中国]-格兹尔算法 [引用日期2026-10-14]
- \[X[k]=\sum_{n=0}^{N-1}x[n]W_{N}^{kn}\] (1) [引用日期2026-10-14]
- references [引用日期2026-10-14]
- The sliding DFT - Loading web-font TeX/Main/Regular [引用日期2026-10-14]
- F [引用日期2026-10-14]
- DTSP2e_text(PDF) [引用日期2026-10-14]
- US9058209(PDF) [引用日期2026-10-14]
- google.com 上的网页 [引用日期2026-10-14]
- [0050] 图1为本发明实施例提供的一种格策尔算法下的频谱确定方法的流程图,如图1所示,该方法包括: [引用日期2026-10-14]