当我们想计算某个频段能量时，FFT一般是首选项。但是当我们并不关注所有的频段，仅仅关注某一个频段，FFT较大的计算量使得其并不是一个较好的选择，那么有没有不用FFT就可以计算频段能量的方法呢？答案当然是肯定的，这就是我们今天要介绍的戈泽尔算法（Goertzel algorithm）。

I、戈泽尔算法

戈泽尔算法把离散傅立叶转换看成是一组滤波器，将输入的讯号与滤波器中的脉冲响应做卷积运算，求得滤波器的输出，即可得到频率域其中一点的频率。此算法利用旋转因子的周期性，将离散傅立叶转换转换为线性的滤波运算。首先回顾DFT的定义

​其中

​我们知道W是周期为N的函数，那么公式（1）可以改写成如下形式

​如果用变量n替换掉常量N，可以得到

​其中 Fk(n)在戈泽尔算法中定义为

​写成卷积的形式

​公式(5)可以看成是f(n)与复指数exp(-j2kπ/N)的卷积，其中k表示时间序列的中心频率对应的索引，这意味着DFT的输出可以视为时间序列经过窄带滤波器的输出。

II、计算步骤

Step1: 当我们确定了采样率和目标频率之后，我们需要预先计算一些常量：

​Step2: 然后我们定义Q0，Q1，Q2三个变量

Step3：我们安装如下的方式进行计算

​Step4: 最后可以得出目标频率的幅值

​从计算步骤可以看出，戈泽尔算法实质就是一个两极点的IIR滤波器。

III、应用

戈泽尔算法有什么用呢，电话系统中电话机与交换机之间的一种用户信令，通常用于发送被叫号码，使用的是双音多频信号（dual-tone multifrequency，DTMF）。双音多频的拨号键盘是4×4的矩阵，每一行代表一个低频，每一列代表一个高频。每按一个键就发送一个高频和低频的正弦信号组合，比如'1'相当于697和1209赫兹(Hz)。交换机可以解码这些频率组合并确定所对应的按键，解码采用的正是基于戈泽尔算法。

​我们首先生成号码1234对应的DTMF信号，得到的结果如下图所示

​然后我们采用戈泽尔算法进行解码，简单来说就是计算DTMF对应频率的能量，其中最大值就对应着DTMF的号码

goertzel_filter(CONSTANTS[0], CONSTANTS[1], internalArray, &T[0], &T[1], SAMPLES); 
goertzel_filter(CONSTANTS[2], CONSTANTS[3], internalArray, &T[2], &T[3], SAMPLES); 
goertzel_filter(CONSTANTS[4], CONSTANTS[5], internalArray, &T[4], &T[5], SAMPLES); 
goertzel_filter(CONSTANTS[6], CONSTANTS[7], internalArray, &T[6], &T[7], SAMPLES); 
goertzel_filter(CONSTANTS[8], CONSTANTS[9], internalArray, &T[8], &T[9], SAMPLES); 
goertzel_filter(CONSTANTS[10], CONSTANTS[11], internalArray, &T[10], &T[11], SAMPLES); 
goertzel_filter(CONSTANTS[12], CONSTANTS[13], internalArray, &T[12], &T[13], SAMPLES); 
goertzel_filter(CONSTANTS[14], CONSTANTS[15], internalArray, &T[14], &T[15], SAMPLES); 
goertzel_filter(CONSTANTS[16], CONSTANTS[17], internalArray, &T[16], &T[17], SAMPLES);

本文代码在公众号语音算法组菜单栏点击Code获取。

参考文献：

[1]. Handbook of Digital Processing

[2]. https://www.embedded.com/the-goertzel-algorithm/

[3]. https://github.com/cleversonahum/dtmf-generator/blob/main/dtmf-generator.py