浅谈离散傅里叶变换简化算法--Goertzel算法

契机是之前对Matlab与Mathematica实现FFT算法(Fast Fourier Transformation快速傅里叶变换算法)感到好奇, 加之学习《数字信号处理》这门基础课, 于是就好好研究了下这部分内容. 边看边写, 目测写不到什么高深的内容里去. 说了半天Goertzel算法又是什么东西呢?在讨论FFT之前, 这次先来谈谈这个Goertzel算法.

在Matlab中实现离散傅里叶变换用fft函数:

a0 = fft(f1)/length(f1)*2;

在Mathematica中实现离散傅里叶变换用Fourier函数:

a1 = 2 Fourier[f1, FourierParameters -> {-1, -1}];

或

a2 = 2 Fourier[f1, FourierParameters -> {1, -1}]/Length[f1];

Goertzel算法简单说就是计算离散时间信号傅里叶变换的一种方法, 由Gerald Goertzel 在1958年首次提出. Goertzel算法是直接计算离散信号傅里叶变换的简化, 从效率上说, 这种算法所需要的计算实数乘法次数大约是直接计算的一半, 而加法次数稍逊于直接计算. 尽管如此, 这在当时主要靠人工手算的时代, 还是十分欢迎. 当N值很大, 正比于$N^{2}$次的浮点运算量十分庞大, 直到1965年J. W. Cooley和J. W. Tukey合作发表了快速傅里叶变换算法后, 计算量下降了几个数量级, 才开启了数字系统处理信号的时代.

虽然Goertzel算法的地位与FFT算法比起来像是个配角, 但当我们不需要计算所有点的傅里叶变换时, Goertzel算法是可取的, 理解实现起来也相对容易, 因此比较适合小型处理器中实现.

计算离散傅里叶变换还得从源头说起, 也就是直接手算开始, 并以此为参照, 容易体会到计算的简化.

$$X[k]= \sum\limits_{n=0}^{N-1} x[n] W_{N}^{kn}\quad k=0,1,...,N-1$$

这是大家熟悉的长度为N有限长序列的离散傅里叶变换计算方法（以下简称DFT）,用这个式子我们可以讨论下直接计算X[k]的计算量.

考虑最大计算量时视x[n]与$W_N^{kn}$皆为复数.

$$x[n]W_{N}^{kn}=(Re\{x[n]\}Re\{W_N^{kn}\}-Im\{x[n]\}Im\{W_N^{kn}\}) \\ \quad \quad \quad \quad +j(Re\{x[n]\}Re\{W_N^{kn}\}+Im\{x[n]\}Im\{W_N^{kn}\})$$

可以看到两个复数相乘需要4次实数乘法与2次实数加法. 计算每一个k的X[k]需要$4N$次实数乘法以及$2N+2(N-1)=4N-2$次实数加法, 因此如果计算所有N个k值再分别乘以N就行了, 也就是$N^{2}$级的计算总量.

后面改善DFT计算效率的算法大多利用了$W_N^{kn} (=e^{-j(2 \pi /N)kn})$这个因子的对称性与周期性, 也就是:

$$W_N^{k(N-n)}=W_N^{-kn} \quad$$ $$W_N^{kn}=W_N^{k(N+n)}=W_N^{(k+N)n}$$

以上两式都容易看出来, 第一个式子说明了这个因子的对称性, 第二个式子说明这个因子对于$k$和$n$有周期性, 周期为$N$.

---

Goertzel算法便是利用了周期性改善计算的栗子.

首先根据周期性:

$$W_N^{-kN}=e^{j(2 \pi /N)Nk}=e^{j2\pi k}=1$$

在计算的的目标$X[k]$中插入这个因子

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

然后我们需要构造这么一个序列, $u[n]$为单位阶跃序列

$$y_k [n]=\sum\limits_{r=-\infty}^{\infty}x[r] W_N^{-k(n-r)} u[n-r]$$

发现当$n=N$时

$$y_k [N]=\sum\limits_{r=-\infty}^\infty x[r] W_N^{-k(N-r)} u[N-r]$$

由于$x[n]$是有限长序列, 对于上式, r只在[0, N-1]取的值$x[r]$才不为0, 因此, 对比之后发现, $X[k]=y_k [N]$. 这样就把序列$y_k [n]$与$X[k]$联系起来了.

回过头再来看我们构造的离散序列$y_k [n]$, 可以看成有限长序列$x[n]$与$W_N^{-kn} u[n]$序列的离散卷积和是不是. 联想到离散卷积可以求得系统响应, 故$y_k [n]$也可以看作单位脉冲响应为$W_N^{-kn} u[n]$的系统对输入$x[n]$的响应, 这点对于之后的推导有重要意义.

利用上述这一点, 能够表述出系统的差分方程, 这个我们能从系统单位脉冲响应的Z变换中得出.

$$H(z)=Z\{W_N^{-kn} u[n]\}=\frac{1}{ 1-W_N^{-k} z^{-1} }$$

我们知道系统函数$H(z)$能够描述离散LTI系统, 并从中得出系统差分方程的系数, 关系是这样的(可以从对系统差分方程两边同时Z变换得到):

$$\sum\limits_{k=0}^{N}a_k y[n-k]=\sum\limits_{k=0}^{M}b_k x[n-k]$$ $$H(z)=\frac{Y(z)}{X(z)}=\frac{\sum\limits_{k=0}^M b_k z^{-k} }{\sum\limits_{k=0}^N a_k z^{-k} }$$

对比就可得到描述这个LTI系统差分方程的系数, 即

$$H(z)=\frac{1}{ 1-W_N^{-k} z^{-1} }=\frac{b_{0}z^{0}}{a_{0}z^{0}+a_{1}z^{-1}}$$

于是, 可写出差分方程:

$$y_k [n]=W_N^{-k} y_k [n-1]+x[n]$$

它是一个递推式, 这意味着每一个$y_k [n]$的计算都可由前一个计算值计算. 若画出递推计算的信号流图, 大概是这样的:

这是Mathematica做的信号流图, 画了一下午看来还是没能画出想要的效果, 我只能说用GraphPlot画比之前用Graph函数画的效果好一点´_>`.

有了简单明了的信号流图, 再来看看求$X[k]$的计算量. 计算每一个新的$y_k [n]$需要4次实数乘法与4次实数加法, 为了递推计算到$y_k [N]$需要计算之前的N-1个值, 也就是大约4N次实数乘法, 4N次实数加法, 而到这只是计算了X[k]的其中一个k的值. 如此看来乘法次数差不多而加法次数比直接计算稍逊一筹.

然而这种递推的计算方法也带来了改善效率的可能, 通过变换离散信号结构使用直接II型来实现系统, 就能得到令人满意的结果.

之前我们得到了$H(z)$(因为是某个k值下的, 其实更应该写做$H_k (z)$), 进行如下变形:

$$H(z)=\frac{1-W_N^{-k} z^{-1}}{(1-W_N^{-k} z^{-1})(1-W_N^k z^{-1})}=\frac{1}{1-2Cos(2\pi k/N) z^{-1}+z^{-2} } (1-W_N^{-k }z^{-1})$$

$H(z)$重新看作两个系统级联:

$$H(z)=H_{1}(z)H_{2}(z)$$

系统函数分别为:

$$H_1 (z)=\frac{1}{1-2Cos(2\pi k/N) z^{-1}+z^{-2} }$$ $$H_2 (z)=1-W_N^{-k} z^{-1}$$

用之前的方法可以求得两个系统的差分方程, 引入中间信号$v_k[n]$:

$$H_{1}(z)=\frac{V_{k}(z)}{X(z)}$$ $$H_{2}(z)=\frac{Y_{k}(z)}{V_{k}(z)}$$ $$v_k [n]=2Cos(2\pi k/N) v_k [n-1]- v_k [n-2]+x[n]$$ $$y_k [n]=v_k [n]-W_N^{k} v_k [n-1]$$

此时系统的信号流图可表示为:

图中可以发现为了计算$y_k [N]$不用再计算之前的所有n从0到N-1的值了, 从图上看只需左边回路循环计算N次, 右边回路计算一次就行了. 也就是:

$$y_k [N]=v_k [N]-W_N^{k} v_k [N-1]$$

对于计算每一个$v_k [n]$, 由于计算式中的系数都是实数, 需要两次实数乘法与4次实数加法, 迭代N次需要2N次实数乘法与4N次实数加法, 接着计算一次右边回路需要4次实数乘法与4次实数加法. 总结一下, 计算某个k值下$y_k [N]$一共需要大约$2N+4$次实数乘法与$4N+4$次实数加法, 因此如果计算所有k值再乘以N, 大约需要$2N(N+2)$次实数乘法与$4N(N+1)$次实数加法, 说了半天可见实数乘法计算量约是直接计算的一半.

Goertzel算法在N值较小或只是计算部分N值时能体现出这种方法的优势. 另一个不得不提的优势来自于迭代的计算方法, 这使得计算在输入第一个样值点后就可以开始.

---

下面是几个简单实现的例子, 先来最近折腾的Mathematica

(*Mathematica*)
Clear[Gortzel]
Gortzel[signal_List]:= 
  Block[{n = Length[signal], ini = signal // First, signallist = signal~Append~0}, 
   Table[#2 - E^(-I (2 Pi/n) k) #1&@@ Block[{i = 2}, 
Nest[{#[[2]], N[2 Cos[(2 Pi k)/n]#[[2]]] - #[[1]] + signallist[[i++]]} &, {0, ini}, n]], 
{k, 0, n - 1}]];

算出结果没问题, 但不够高效, 计算几百个点的序列就得等好一会儿了.

%Matlab
function[xk]=Goertzel(xn)
clear xk;
signallist = [xn, 0];
N=length(xn);
xk=ones(1,N);
for k = 0:N-1
    box=[0, signallist(1)];
    for ii=2:N+1
        box=[box(2), 2*cos(2*pi*k/N)*box(2)-box(1)+signallist(ii)];
    end
    xk(k+1)= box(2)-exp(-1j*2*pi*k/N)*box(1);
end

看了下原来Matlab自带这个算法, 使用goertzel函数, 试了下还挺高效.

//C++
#include <iostream>
#include <cmath>
#include "mycomplex.h"
#include <vector>

const double E=2.718281828459;
const double PI=3.1415926535898;

using namespace std;

vector<Complex> Goertzel(vector<Complex> xn)
{
    vector<Complex> xk;
    int n = xn.size();
    xn.push_back(0);
    for (int k = 0; k < n; k++)
    {
        Complex box[2] = { 0, xn[0]};
        for (auto it = xn.begin() + 1; it != xn.end(); it++)
        {
            Complex temp = box[1];
            box[1] = box[1] * 2 * (cos(2 * PI*k / n)) - box[0] + *it;
            box[0] = temp;
        }
        xk.push_back(box[1] - Complex(cos(2 * PI*k / n), -sin(2 * PI*k / n)) * box[0]);
    }
    return xk;
}

参考: 《离散时间信号处理（第三版）》Alan V. Oppenheim, Ronald W. Schafer 著.
.
.
.
.
.
.
.
.
.
終わり