Goertzel 算法及其应用

5.9 Goertsel 算法及其应用
5.9.1 Goertsel 算法
如前所述, N 点时域序列)(n x 的离散付里叶变换式为
∑=-=-1
0)()(N n kn
N W n x k X ， 1,....,2,1,0-=N k （5.9.1）
这N 点频域序列是同时被算出的，不可能只计算其中某一个或几个指定点。

Goetzel 算法是为了解决
这个问题而提出的。

这个算法把离散付里叶变换看作一组滤波器，将输入端的时域序列与其中一个滤波器的冲激响应序列进行卷积运算，求滤波器的输出序列，即得)(k X 序列的一点。

这种算法利用旋
转因子k
N W 的周期性，使DFT 运算化为线性滤波运算。

由于
1))(2(
==-π
--kN N
j kN N
e
W
故式（5.9.1）可化为 )
(1
10
)()()(m N k N N m km N N m kN
N
W m x W
m x W
k X ---=-=-∑=∑=， 1,....,2,1,0-=N k （5.9.2） 定义序列)(n y k 为
)
(1
0)()(m n k N
N m k W m x n y ---=∑= （5.9.3A ）
可见)(n y k 是由两个序列卷积而得到的序列。

)()()(n h n x n y k k ⊗= （5.9.3B ） 其中，)(n x 是输入的N 点序列，另一个序列被看作滤波器的冲激响应序列
)()(n u W n h kn
N k -=。

（5.9.4）
对比式(5.9.2）和式(5.9.3），可知：按式（5.9.3A ）进行卷积运算，当N n =时，滤波器的输出)(N y k 就是)(k X ：
N n k n y k X ==|)()( （5.9.5）
对式(5.9.4)进行Z 变换，可得滤波器的系统函数
1
11)(---=
z
W
z H k k （5.9.6）
这是一个一阶系统。

图5.9.1示出这个系统的信号流图，相应的差分方程为
)()1()(n x n y W n y k k
N k +-=-， 0)1(=-y (5.9.7)
按照此式进行递推运算，到了 N n = 时刻，即可依据式（5.9.5）得到)(k X 。

按照式(5.9.6)进行运
算时，可先算好旋转因子k N W -，储存起来。

每次递推包含一次复数乘法。

按式(5.9.1)直接计算N 点离
散付里叶变换,需要2
4N 次实数乘法和)N-N(24次实数加法。

按照上述Goertzel 算法，所需的实数乘法和实数加法都是2
4N 次。

所以当N 不大时，上述算法的效率稍差。

下面介绍改进的Goertzel 算法，这种算法所需的实数乘法次数约为直接方法的一半。

图5.9. 1 用一阶系统实
现Goertze
l 算法 图 5.9.2 用二阶系统实现Goertzel 算法
把式（5.9.6）的分子和分母都乘上因子)1(1
--z W k N ，就得到第k 个滤波器的系统函数为
)
1)(1(1)(1
11
-------=z W z W z
W z H k N k N k N k 2
11
))/2cos((211---+π--=z
z k N z
W k N （5.9.8） 与此相应的信号流图示于图5.9.2。

由式(5.9.8)可见，滤波器是一个二阶系统，有一对复数共轭极点和一个复数零点。

为了便于运算，在图5.9.2所示的流图中，设立状态变量 u 和v 。

按照图5.9.2计算)(k X
时，步
骤有二，即：
1． 实现一对复数极点
输入点依次取 )1(),...,2(),1(),0(-n x x x x ，
进行递推运算。

每次运算中，更新状态变量u 和v 。

作N 次迭代所需的计算量是N 2次实 数乘法和N 4次实数加法。

2． 实现复数零点。

)(n x 是一个N 点序列， 1,...,2,1,0-=N n 。

在0)(=N x 点上。

计算状态变量u 和v 。

这时，按照图5.9.2算出滤波器的输出)(N y k ，此即)(k X 。

所需的计算量是4次实数乘法和N 4次实数加法。

综上所述，计算一点)(k X 需要进行)2(2+N 次实数乘法和)1(4+N 次实数加法。

这种算法要求的乘法次数约为直接算法的一半。

在这种较为有效的方案中，仍具有这样的优点，即必须计算和存
储的系数只有))/2cos((k N π和k N W 。

最后，还要说明图5.9.2所示的算法的另一个优点。

当输入序列为实序列时，离散付里叶变换序列
)(k X 是对称的，即)()(*k N X k X -=。

容易证明，图5.9.2的网络形式在计算)(k N X -时和计算)(k X 时具有完全相同的极点，但前者的零点系数与后者的零点系数成复共轭关系。

由于零点仅在最后的迭代中实现，所以诸极点要求的N 2次乘法和 N 4 次加法可以用来计算离散付里叶变换的两个
值。

因此，若用Goertzel 算法计算离散付里叶变换的所有N 个点的值，需要的乘法次数近似为2
N ，加法次数近似为2
2N 。

然而，它同直接计算离散付里叶变换一样，计算量仍然正比于2
N 。

5.9.2 双音多频（DTMF ）信号
在数字电话中，每个拨号键产生的信号由两个音频组成，其一是低音频，另一是高音频。

4个低音频（697Hz 、770Hz 、852Hz 和941Hz ）和4个高音频（1209Hz 、1336Hz 、1477Hz 和1633Hz ）分别组成低音组FL 和高音组FH 。

这种信号称为双音多频（dual-tone multifrequency, DTMF ）信号。

可视化程序B0416演示DTMF 信号的产生与检测。

其界面示于图5.9.3。

图中，左下角是DTMF 键盘，由10个数字键和6个功能键组成。

每个数字键产生高、低音频信号各一个。

例如，按下“7”键时，产生频率为Hz 852FL =和Hz 1209FH =的正弦序列。

其合成序列为)(n x 。

用 Goertzel 算法对序列
)(n x 进行解码，得8点频谱|)(|k X ，其中以22=k （相应于Hz 852FL =）和31=k （Hz 1209FH =）
这两条谱线的幅度最大。

由此得知被按下的键为“7”。

下一节介绍用Goertzel 算法实现DTMF 检测。

图5.9.3 DTMF 键盘及其检测
5.9.3用Goertzel 算法实现DTMF检测
在实践中，由于通信线路中还包含语音信号。

而语音信号的特点是二次谐波比较大，所以在线路终端进行解码时，除了检查DTMF的8个频率点外，还要检查它们的二次谐频。

如果二次谐波小，就可以确信这是DTMF信号，进一步检测哪两点的频谱幅度最大，从而确定键值。

由于我们只对DTMF
的8个频率点及其二次谐频感兴趣，所以最好用Goertzel 算法代替DFT进行分析。

图5.9.4 基于Goertzel 滤波器的DTMF 解码原理
图5.9.4示出基于Goertzel 算法的DTMF检测方案的框图。

整个方案由8对Goertzel 滤波器组成。

每对滤波器检出一个DTMF 信号及其二次谐波。

二次谐波用来甄别语音和DTMF 信号。

语音含有很
大的偶次谐波，而DTMF 信号则无。

将每个输出求平方，以测出每个DTMF 频率及其二次谐波的信
号强度。

问题是怎样选择Goertzel 滤波器的参数。

首先，数字电话系统中的采样频率已确定为
Hz 8000S =f 。

这个选择是合适的，因为DTMF 信号的8个频率及其二次谐波频率都落在奈奎斯特
区间中。

其次，主要问题是怎样选择数据点数N 以得到可靠的频率分辨率。

N 越大，频率分辨率就越高。

但是，对于按键的最短时间要求是40ms ，所以，可用的数据点数可能少至800004.0⨯=N 。

在8个基频中，相邻基频差最小值为Hz 73697770=-。

若要保证两个相邻基频序号之差至少为2，则频率分辨率应为5.362
73
=>
D ，故N 应取200以上。

另一方面，在频谱图上，8个基频及其二次谐频由谱点序号k 标识。

由此算出的频率与实际的基频和二次谐频相比，是有误差的。

为了使误差最小，应在300~200=N 范围内再做细致选择。

人们研究的结论是：对于已选定的8个基频，205=N 是最好的选择。

表5.9.1列出详细的数据(引自Ref.12)。

程序
B0416和M050901.M 是基于以上考虑而编制的。

可视化VB 程序B0416的界面已见于图
5.9.3 。

程序M050901.M 运行时，命令窗要求用户输入8 位电话号码code ，之 后，程序检测该码，并予以显示。

例如，若键入 code = “67392477”，则程序显示检测结果，如图5.9.5所示。

图5.9.5 8位电话号码的DTMF 检测结果。