第5章〓快速傅里叶变换

快速傅里叶变换 (Fast Fourier Transform)是指在经典离散傅里叶变换（DFT)的基础上，以降低运算量和提高运算效率为目的，对DFT进行的各种改进，各种改进算法均称为傅里叶变换的快速算法，简称FFT。采用快速算法可以使离散傅里叶变换所需要的乘法次数大为减少，而且序列的点数越大，FFT算法对计算量的减少就越显著。



视频讲解



5．1引言
◆

快速傅里叶变换并不是一种新的变换，而是对DFT的一种算法改进。因此，FFT算法的基础是DFT，核心是对DFT算法进行高效改进。

有限长序列的离散傅里叶变换的时域和频域均为离散形式，可以采用计算机进行计算。DFT在数字信号处理中应用非常广泛。例如，FIR滤波器设计需要由h(n)求H(k)以及由H(k)求h(n)，这些都需要进行DFT计算。信息通信、图像传输、雷达、声呐等领域需要进行信号频谱分析，这也需要进行DFT计算。在数字系统分析与设计中，也经常要进行DFT计算。然而，在20世纪60年代中期以前，由于DFT的计算量非常大，且计算机的运算效率低，故很少采用计算机进行DFT运算。1965年，库利（J.W.Cooley）和图基(J.W.Tukey)在《计算数学》(Mathematics of Computation)杂志上发表了著名的“机器计算傅里叶级数的一种算法”的论文，提出了DFT的一种快速算法，傅里叶变换的快速算法才引起广泛重视，从此，快速傅里叶变换的应用情况发生了根本性的变化，FFT经过不断完善和改进，其运算时间一般可缩短至两个数量级以上，运算时间的有效缩短有效推动了FFT的广泛应用。

党的二十大报告提出“实现高水平科技自立自强，进入创新型国家前列”，实现这一目标，需要我们发奋努力，在科学和技术的各个领域不断创新、不断突破。快速傅里叶变换具有非常广泛的应用，尖端医疗器械、先进电力控制设备、5G通信技术、星空探测、航母控制系统、高铁信号处理系统等都离不开FFT。


5.2DFT的运算量及改进的途径
◆


5.2.1DFT算法运算量的估计 

设序列x(n)为N点有限长序列，其DFT算法如下。

正变换（DFT）： 


X（k）=∑N-1n=0x(n)WnkN,k=0,1,…,N-1（51）


逆变换（IDFT）:


x(n)=1N∑N-1n=0X(k)W-nkN,n=0,1,…,N-1（52）


比较DFT和IDFT可以发现，二者的算法结构非常类似，差别很小，它们的差异性主要体现在如下两点。

（1） 正变换和逆变换算法中WN的指数符号相反； 

（2） 正变换求和符号前的系数为1，而逆变换求和符号前的乘数因子为1/N。

由上述分析可知，正变换和逆变换的运算量几乎相等，因此，一般以DFT正变换为例进行运算量的估算。由于运算结构相同，故逆变换可以根据正变换的运算量进行准确估计。


为了提高适应性，一般假定序列x(n)在复数范围内取值。实际上，无论x(n)是实数还是复数，WnkN一般都是复数、因而X(k)（k=,0,1,2,…,N-1）一般也是复数。

根据DFT正变换公式，对于N点频域序列X(k)而言，每计算一个X(k)值，需进行N次复数乘法，即求和符号中x(n)与WnkN有N次复数乘，（N-1）次复数加； 而 k=0,1,2,…,N-1，即X(k)一共有N个值，因此 DFT正变换的全部运算共有N2次复数乘法，（N（N-1））次复数加法。

复数运算是以实数运算为基础来完成的，因此，DFT正变换可写为如下形式。


X（k）=∑N-1n=0x(n)WnkN=∑N-1n=0{Re［x(n)］+jIm［x(n)］}{Re［WnkN］+jIm［WnkN］}
=∑N-1n=0Re［x(n)］Re［WnkN］-Im［x(n)］Im［WnkN］+
j(Re［x(n)］Im［WnkN］+Im［x(n)］Re［WnkN］
（53）


根据式（53），可得： 

（1） 一次复数乘法=4次实数乘法+2次实数加法； 

（2） 一次复数加法=2次实数加法。

因此，每计算一个X(k)需4N次实数乘法以及2N（2N-1）次实数加法。所以N点的DFT总计算量包括4N2次实数乘法和2N（2N-1）次实数加法，如表51所示。


表51N点DFT运算量表



计算方式
复数乘法
复数加法
实数乘法
实数加法



运算量
N2
N2-N
4N2
4N2-2N




上述运算量估计是纯理论估计值，它包含了乘以常数1等特殊情况，与实际运算量会略有差异。这是因为式（51）中时域因子n和频域因子k取不同值时，WnkN可能等于1、-1或j，这时无须乘法运算，如： 


W0N=1，WN/2N=-1，WN/4N=-j


这些特殊情况会使实际运算量略小于运算量的理论值，当N越大时，这种特例所占的比例就越小。因此，为了统一和方便比较，DFT运算量估计一般不考虑上述特殊情况。

根据表51，完成一个长度为N点的DFT计算，乘法次数和加法次数都与N2成正比，当N很大时，运算量是非常大的。例如，若N=8，则DFT运算量为64次复数乘法； 若N=1024，则DFT运算量为1048576次复数乘法，如果系统对实时性要求很高，这样巨大的运算量对计算机的速度要求非常高。若采样点数N进一步增大，则运算量以平方关系快速增长，因此，需要对DFT的算法进行改进，以有效降低DFT算法的运算量，提高运算速度。



视频讲解


5.2.2降低运算量的途径

任何改进都必须从可改之处着手，从最有效之处改进，DFT算法的改进也不例外。根据5.2.1节的分析可知，影响DFT运算量的因素有如下两个。

（1） WnkN等于1、-1等常数时可适当降低运算量。

（2） 序列x(n)的长度N（即点数）对运算量有显著影响。

下面分析如何利用这两个因素减少运算量。

1． 利用WnkN减少运算量

利用WnkN等于1、-1等常数，以及WnkN的相关性质，可以适当降低运算量，主要包括如下特性。

（1） WnkN=1，-1，…等常数，如，W0N=1，WN/2N=-1，WN/4N=-j。

（2） WnkN的对称性。


（WnkN）*=W-nkN


（3） WnkN的周期性。


WnkN=W（n+N）kN=Wn（k+N）N


（4） WnkN的可约性。


WnkN=WmnkmN,WnkN=Wnk/mN/m


由此可得


Wn（N-k）N=W（N-n）kN=W-nkN,Wn/2N=-1,W(k+N/2)N=-WkN


利用上述特性，可以使DFT运算中某些项无须乘法运算，也可以使某些项合并运算。

2． 降低序列长度N减少运算量

由于运算量与N2成正比，因此，减小N对减少DFT运算量是最有效的，它将使得运算量以平方的速度快速下降。然而，N的减小是受限的，即受到奈奎斯特采样定理约束，序列长度N不能无原则缩短，工程上，一般应在满足采样定理的基础上，选择符合采样要求的最短序列。

DFT的运算量与N2成正比，N2是非线性函数，N减小时运算量显著下降，即点数少的序列DFT运算量比点数多的序列的运算量要小很多。因此，当点数N不能进一步减小时，降低运算量的思路是将N点的序列x(n)分解为2个N/2点的序列，这样两个短序列的复数乘法运算量均为


N22+N22=N22


经过一次分解，复数乘法运算量大约降低了一半。以此思路继续分解，每一个N/2点的序列又可以分解为2个N/4点的DFT，运算量可进一步减小，以此类推，对每一个新产生的短序列继续分解，直至分解到全部子序列均为2点为止。


DFT算法改进主要依据这一思路展开，因此快速算法理论上主要可以分成两大类，一种是按时间抽取(DecimationInTime，DIT)法，将长序列分解为短序列，另一种是按频率抽取(DecimationInFrequency，DIF)法，将长序列分解为短序列。


5．3按时间抽取（DIT）基2 FFT算法
◆

FFT快速算法的思路是将长序列分解为短序列，并根据DFT变换的奇、偶、虚、实的对称性及相关特性，对算法进行改进从而得到快速算法。FFT算法并没有在理论上对傅里叶变换有新的推进，但在推动傅里叶变换的应用方面却取得了巨大的进步。根据将序列分解为短序列的方法不同而产生了不同的FFT算法，典型的算法主要包括DIT基2（按时间抽取）和DIF基2（按频率抽取）算法。

5.3.1DIT算法原理

对序列x(n)按序号的奇、偶抽取得到两个子序列，并不断进行分组而实现的FFT算法称为DIT基2 FFT算法。该算法由库利图基首先提出，故又称为库利图基算法。

设序列x(n)的点数为N=2L，L为正整数。若序列点数不满足条件，可以从序列末端补零点，使之满足要求。N为2的整数幂的FFT称为基2 FFT。




视频讲解




1． 序列x(n)第一次分解

将长度为N=2L的序列x(n) (n=0，1，…，N-1)，根据时域序号n的奇、偶分为如下两组（两个子序列）： 



x1(r)=x(2r)x2(r)=x(2r+1)r=0,1,…,N2-1
（54）


于是，DFT正变换公式可以化为



X（k）=DFT［x(n)］=∑N-1n=0x(n)WnkN=∑N-1n=0n为偶数x(n)WnkN+∑N-1n=0n为奇数x(n)WnkN
=∑N2-1r=0x(2r)W2rkN+∑N2-1r=0x(2r+1)W(2r+1)kN
=∑N2-1r=0x1(r)(W2N)rk+WkN∑N2-1r=0x2(r)(W2N)rk
（55）


由于W2N=e-j2πN·2=WN/2，式（55）进一步可以化为


X(k)=∑N2-1r=0x1(r)WrkN/2+WkN∑N2-1r=0x2(r)WrkN/2（56）


由于x1(r)、x2(r)均为N/2点的时域序列，其对应的DFT如下： 



X1（k）=∑N2-1r=0x1(r)WrkN/2=∑N2-1r=0x(2r)WrkN/2
X2（k）=∑N2-1r=0x2(r)WrkN/2=∑N2-1r=0x(2r+1)WrkN/2
（57） 


因此，式（56）可表示为如下形式:


X(k)=X1(k)+WkNX2(k)（58）


X1(k）、X2(k）都是长度为N/2点的序列x1(r)、x2(r)的标准DFT正变换，至此，一个N点DFT已分解成两个N/2点的DFT。

（1） X1(k）、X2(k）都是N/2点的标准DFT正变换，因此，一个N点DFT已分解成两个N/2点的DFT。

（2） 式（58）中x1(r)、x2(r)以及X1(k）、X2(k）都是N/2点的序列，即r，k=0，1，…，N/2-1，而X(k）却有N点，因此，根据式（58）只能计算X(k）的前一半点的结果，而X(k）后一半点的结果需要根据DFT的周期性以及WnkN的周期性对式（58）进行改进。

根据WnkN周期性


WrkN/2=Wrk+N2N/2（59）


代入式（57）可得


X1N2+k=∑N2-1r=0x1(r)WrN2+kN/2=∑N2-1r=0x1(r)WrkN/2=X1(k)（510）


对于X2(k)，同样可得


X2N2+k=X2(k)（511）


由于WN/2N=-1，可得


WN2+kN=WN/2NWkN=-WkN（512）


因此，将式（510）~式（512）代入式（58），可得后一半点的公式如下： 


Xk+N2=X1k+N2+Wk+N2NX2k+N2
=X1(k)-WkNX2(k),k=0,1,…,N2-1（513）


由此可得计算全部X(k）值的算法如下： 



X（k）=X1(k)+WkNX2(k),k=0,1,…,N2-1
Xk+N2=X1(k)-WkNX2(k),k=0,1,…,N2-1
（514）


根据式（514），由长度为N/2点的X1(k）和X1(k），即可求出N点X(k）的全部值。



图51按时间抽取法蝶形运算流图


式（514）的运算规律可用图51的蝶形运算流图进行描述，式（514）中的乘法运算，对应于系数为WkN的蝶形支路。采用计算流图表示算法结构更简明扼要，也更加形象化，若支路上标有系数，则表示信号与支路系数相乘，若支路没有标注系数，则表示该支路的系数为1。



图51是一个基于DIT基2快速算法的一般蝶形运算图，若以k=0,1,2，…，N2-1代入图51，则可得到FFT快速运算的第一次分解蝶形运算图。为便于深入理解FFT快速算法原理，以N=23=8为例，详细分析图51所包含的全部碟形运算。

当k=0时，碟形运算如下： 


X(0)=X1(0)+W0NX2(0)
X(4)=X1(0)-W0NX2(0)


当k=1时，碟形运算如下： 


X(1)=X1(1)+W1NX2(1)
X(5)=X1(1)-W1NX2(1)


当k=2时，碟形运算如下： 


X(2)=X1(2)+W2NX2(2)
X(6)=X1(2)-W2NX2(2)


当k=3时，碟形运算如下： 


X(3)=X1(3)+W3NX2(3)
X(7)=X1(3)-W3NX2(3)


这就是图51或式（514）所包含的全部蝶形运算，对于上述表达式中的X1(k）、X2(k）(k=0,1,2,3)的计算，依然采用传统DFT公式进行计算。根据上面k=0,1,2，3的4个表达式以及式（514），可以得到8点DFT经过第一次按序号n的奇偶分组（分解）之后计算X(k）的过程,如图52所示。



图52一个N点DFT分解为两个N/2点DFT算法过程


图52中，标注DFT的两个方框表示依然采用DFT公式进行计算，即X1(0）~X1(3）采用式（57）中的第一式计算，X2(0）~X2(3）采用式（57）中的第二式计算，而最右端一列的X(0）~X(3）采用式（514）中的第一式计算，X(4）~X(7）采用式（514）中的第二式计算。可以看出，每一个蝶形运算仅需一次复数乘法X2(k）WkN及两次复数加法运算。

图52完整地表示了一个N点DFT分解成两个N/2点DFT的过程，该图完全等价地描述了式（514）的全部计算过程。根据第一次分解过程，若直接计算N/2点DFT，则每一个N/2点DFT有N22=N24次复数乘法，N2N2-1次复数加法，两个N/2点DFT共需2×N22=N22次复数乘法和NN2-1次复数加法。

此外，采用蝶形运算将两个N/2点DFT组合成一个N点DFT时，有N/2个蝶形运算，还有N/2次复数乘法及2×N/2=N次复数加法。因而第一次分解完成之后，复数运算量如下。

复数乘法运算量： 


N22+N2+N（N+1）2≈N22（次）


复数加法运算量： 


NN2-1+N=N22


由此可知，计算N点序列x(n)的频域序列X(k)，在经过第一次分解之后，运算量大约为直接采用DFT运算量的一半。



视频讲解


2． 序列x(n)第二次分解

N点序列分解为两个N/2点短序列可以降低运算量，以此类推，由于N=2L，N/2仍是偶数，对x1(r)、x2(r)两个N/2点子序列应继续按序号r的奇偶分别分解为两个N/4点子序列。

1） 对x1(r)分解

按第一次分解的原理，先将x1(r)分解为x3(r)、x4(r)。


x1(2l)=x3(l)
x1(2l+1)=x4(l),l=0,1,…,N4-1
（515）
X1(k)=∑N4-1l=0x1(2l)W2lkN/2+∑N4-1l=0x1(2l+1)W(2l+1)kN/2
=∑N4-1l=0x3(l)WlkN/4+WkN/2∑N4-1l=0x4(l)WlkN/4
=X3(k)+WkN/2X4(k),k=0,1,…,N4-1
（516）


式（516）可以求出X1(k)前一半点的值k=0,1,…,N4-1，对于X1(k)的后一半点，根据第一次分解的方法，可得


X1N4+k=X3（k）-WkN/2X4(k),k=0,1,…,N4-1


式中 



X3（k）=∑N4-1l=0x3(l)WlkN/4
X4（k）=∑N4-1l=0x4(l)WlkN/4
（517）


与第一次分解类似，对x1(r)的分解过程也可以用蝶形图进行形象的表述，图53（a）给出了当N=8时，将一个N/2点的DFT分解成两个N/4点DFT，并由这两个N/4点DFT组合成一个N/2点DFT的过程。



图53N/2点DFT分解为两个N/4点DFT


2） 对x2(r)分解

对x2(r）也可进行同样的分解： 



X2(k)=X5(k)+WkN/2X6(k)
X2N4+k=X5（k）-WkN/2X6(k),k=0,1,…,N4-1



式中



X5（k）=∑N4-1l=0x2(2l)WlkN/4=∑N4-1l=0x5(l)WlkN/4
X6（k）=∑N4-1l=0x2(2l+1)WlkN/4=∑N4-1l=0x6(l)WlkN/4
（518）


与x1(r)类似，x2(r)的分解过程也可以用蝶形图描述，图53（b）给出了当N=8时，序列x2(r)分解成两个N/4点的DFT，并由这两个N/4点DFT组合成一个N/2点DFT的过程。

至此，一个N点（N=8）DFT分解成了四个N4=2点的DFT。经过两次分解之后，由x(n)计算X(k)的运算过程如图54所示。为了统一支路上的乘数因子（又称为旋转因子，也称为支路系数），根据WnkN的可约性，对WnkN类型的支路系数进行了适当的运算，如WkN/2=W2kN。图54中，标注DFT的方框，无论序列长度是多少，均表示采用经典DFT计算。



图54一个N点DFT分解为四个N/4点DFT（N=8）


显然，序列x(n)经过第二次分解之后，由四个N/4点的DFT以及右边的两级蝶形组合运算来计算N点DFT，其运算量比只有第一次分解的运算量又降低了大约一半。



视频讲解


3． 序列x(n)第三次分解

以N=23=8为例，经过两次分解以后，图54左边一列为四个N/4点的DFT，4个子序列的长度均为2点，最左列每一框图都是2点DFT运算，其输出分别为X3(k)、X4(k)、X5(k)、X6(k)，k=0，1，其中X3(k）、X4(k)依据式（517）进行计算，X5(k）、X6(k)依据式（518）进行计算。

根据DFT的公式，展开2点DFT公式可以发现，完成2点DFT的计算只需1次复数乘法（由于W02=1，实际上无须乘运算）和2次加法运算，将图54左上第一个方框中的N/4点DFT按公式展开，形式如下： 


X3(0)=x(0)+W02x(4)
X3(1)=x(0)-W02x(4)
（519）


以N=8为例，根据式（519）可得2点长序列的蝶形算法流图，即计算X3(k)的蝶形运算流图如图55所示。



图552点DFT运算流图


当N=8时，图54左列其他3个DFT（分别计算X4(k）、X5(k）、X6(k）)也都是2点DFT，都可以将2点DFT框图直接展开，并以蝶形运算表示，由此可得出当N=8时FFT的完整碟形运算图，如图56所示。



图56N=8时（DIT）FFT运算流图


图56称为FFT快速运算的完整蝶形运算图（N=8时），FFT蝶形图由x(n）计算X(k）是从左至右按列（图中虚线）进行计算的，因此，当N=8时只需L=3级（2L=23）运算。

第1级运算。按照左边第一列4个蝶形图进行计算，序列x(n）经过分组排队后，以［x(0)，x(4)，x(2)，x(6)，x(1)，x(5) x(3)，x(7)］的顺序输入，按照蝶形运算规律，计算出第一级运算的输出值［X3(0），X3(1），X4(0），X4(1），X5(0），X5(1），X6(0），X6(1）］。

第2级运算。按第二列4个蝶形图进行计算，这时第一级的输出就是第二级的输入值，依据第二列的4个蝶形图（上下两组，每组两个蝶形图）进行计算，得出输出值［X1(0），X1(1），X1(2），X1(3），X2(0），X2(1），X2(2），X2(3）］。

第3级运算。按第三列4个蝶形图进行计算，根据输入［X1(0），X1(1），X1(2），X1(3），X2(0），X2(1），X2(2），X2(3）］按最右一列的4个蝶形图计算出8点频域序列X(k),k=0,1,2,…,N-1。

上述3级计算中，每一级计算均包含N/2次复数乘法运算，且每一个基本蝶形运算仅包含一个乘法，因此，若N=8，依据图56所示的FFT蝶形运算计算DFT，则仅需12次复数乘法运算（包括乘W0N=1）。

上述3级蝶形运算中，每一级基本蝶形图的系数具有明显的规律性。

（1） 第1级蝶形运算的系数为W0N=W08。

（2） 第2级蝶形运算的系数为W0N，W28。

（3） 第3级蝶形运算的系数为W0N，W18，W2N，W3N。

以此类推，大家也可以得出当N=16时，DIT基2 FFT快速算法完整的蝶形运算流图以及与流图相对应的FFT4级运算过程，同样的方法，也可以得出N=2L的完整FFT算法流程图以及FFT的L级运算过程。

DIT基2算法的每一次分解都是按输入序列x(n）的序号n是偶数还是奇数进行分解，因此称为按时间抽取法。

5.3.2DIT算法的运算量

根据DIT算法的FFT流图，当N=2L时，蝶形流图从左往右共有L级蝶形运算，每一级都包含N/2个基本蝶形运算单元，每个基本蝶形运算单元包含一次复数乘法、二次复数加法。因此，一个完整蝶形流图在从左至右的L级运算中，每一级都需要N/2次复数乘法和N次复数加法，这样L级完整蝶形运算的总运算量如下。

复数乘法次数： 


mT=N2L=N2log2N
（520）


复数加法次数：  


aT=NL=Nlog2N
（521）


通常情况下，考虑到W0N=1、WN/2N=-1和WN/4N等因素，实际复数乘法运算量一般略低于理论估计值。通过分析N点时域序列的FFT算法过程可以得出，发生W0N=1的情况共有N-11+2+4+…+2L-1=∑L-1i=02i=2L-1=N-1次，这些特殊系数都不用乘法运算，对于工程信号，一般N都比较大，特殊系数在运算量中占比较低。因此，对快速算法运算量的估计一般均按理论值，而不考虑系数为±1等特例。

由于计算机进行乘法运算所需时间比加法运算所需时间长很多，因此，也可以只用复数乘法次数表示运算量的大小。FFT算法与DFT算法的运算量比较如表52所示。


表52DFT算法与FFT算法运算量比较



N
N2
N2log2N
N2N2log2N



2414.0
41644.0
864125.4
16256328.0
3210248012.8
64409619221.4
1381638444836.6
25665536102464.0
5122621442304113.8
102410485765120204.8
2048419430411264372.4




DFT的复数乘法次数是N2，FFT复数乘法次数为N2log2N，直接计算DFT与用FFT算法的运算量（复数乘法次数）之比为


N2N2L=N2N2log2N=2Nlog2N（522）


根据式（520）和式（522）可以得出FFT算法与DFT算法的运算量与序列长度N的关系曲线，如图57所示，由该图可以看出，FFT算法在节省运算量上具有巨大的优越性，且N越大，FFT快速算法节约运量的效果就越明显，当N=1024时，DFT算法的运算量约为FFT运算量的205倍。



图57DFT与FFT乘法运算量比较


5.3.3DIT基2算法的特点

DIT基2算法具有一系列的特点，例如，根据快速算法的分解过程，可以很方便地得出N=4点和N=8点FFT算法流图的递推关系，同样也可以得出N=2L-1点与N=2L点FFT算法流图之间的递推关系。



视频讲解


1． 倒位序规律


在图56中，序列经过L=3级分解，对于N=8点的序列，FFT的输出序列X(k)为正常顺序排列，顺序为X（0），X（1），…，X（7），但输入序列x(n)经过L=3次排序以后，已经不是自然顺序，而是按x(0)，x(4)，x(2)，x(6)，x(1)，x(5) x(3)，x(7)的顺序排列。如果N=16，32，…，1024，则随着N的不断增大，输入信号x(n)被分组的次数越多，经过不断分组以后，x(n）的输入顺序也越复杂。若依据逐次奇偶排序的方法计算输入的顺序，则既耗时间，又容易产生错误。

当N很大时，在完成各次分组排队之后，x(n）的输入顺序有没有规律可循呢？答案是肯定的，这一规律就是倒位序规律。

为了探讨输入的排序规律，可以从x(n）的分组方法和二进制原理的内在关系入手。

回顾x(n)按序号n的偶和奇分组的过程，以N=8为例，n可以用三位二进制数表示，x(n)的序号可表示为（n2n1n0）2，第一次分组，序号n为偶数组成x1(r)，在上半部分； n为奇数组成x2(r)，在下半部分。

从二进制视角看，第一次分组完全由二进制数n的最低位n0决定，n0=0时n为偶数，n0=1时n为奇数。也就是说，x(n)按序号n的偶和奇分组实际上就是按序号二进制数的最低位n0进行分组。

第一次分组之后，得到两个长度为4点的子序列，这时n0已经确定。第二次分组是对4点子序列x1(r)、x2(r)进行分组，这时是对（n2n1n0）2中n2n1进行排队。因此，第二次分组对x1(r)、x2(r)按偶序号和奇序号排队，和第一次分组的道理一样，这时是对二进制数n2n1进行排队，低位是n1，因此由n1决定第二次分组，n1=0时r为偶数，n1=1时r为奇数。

同理，第三次分组的奇偶由最低位n2确定，若N=8，则只有n2，因此，n2=0对应于偶数，n2=1对应于奇数。

以N=8为例，按序号n的偶序号和奇序号进行分组的过程，可用如图58所示的二进制树状图进行描述和解释。



图58按序号分组形成倒位序的机制图


图58中，n0、n1、n2所在的列，正好对应了N=8点序列的3次分组排队的过程，而最左边的x（000），x（100），…，x（111）则是x(n)三次分组后的顺序，观察括号内的三位二进制序号，若从右往左看，则为序列的正常顺序； 而从左往右看，则为倒位序。

因此，倒位序就是DIT基2算法的输入顺序，由N=8点推而广之，若序号n的正常二进制数为n0n1n2…nL-1，则FFT算法的输入为nL-1…n2n1n0，即序号正常二进制的倒位序，根据倒位序规律，可得DIT基2 FFT算法的输入顺序如表53所示。


表53DIT基2 FFT算法倒位序规律表（N=8）



自 然 顺 序
倒位序

自然顺序（n）
二 进 制 数
倒位序二进制数
倒位序顺序（n^）



00000000
10011004
20100102
30111106
41000011
51011015
61100113
71111117



FFT实际运算中，一般先按自然顺序将输入序列x(n)存入存储单元，为了得到倒位序的顺序，可以通过程序实现，以N=8为例，若正常二进制数为（n2n1n0），倒位序二进制数则为（n0n1n2），用程序实现这一功能非常简单。

2．  原位运算

从图56可以看出，FFT运算具有明显的规律，Ｌ级运算中的每一级（每一列）运算都由N/2个蝶形图组成，每一个基本蝶形运算完成如下迭代运算： 


Xm(k)=Xm-1(k)+Xm-1(j)WrN
Xm(j)=Xm-1(k)-Xm-1(j)WrN
（523）


式中，m表示第m列，k、j表示FFT蝶形运算图的第k行和第j行。

因此，式（523）对应的基本蝶形运算如图59所示，包括一次复数乘和两次复数加。



图59按时间抽取蝶形运算结果


根据图56，第m-1列的两个节点k和j的节点变量进行蝶形运算后，得到结果为第m列的k、j两节点的节点变量。在算法的实现上，FFT蝶形运算是从左至右，按列进行运算。以N=8为例，即先由输入序列x(n)，分别采用第一级（列）的4个蝶形运算，求出8个数据； 然后以这8个数据作为输入，采用第二级的4个蝶形运算图，求出8个数据； 再采用第三列蝶形运算图，求出最终的X(k)的全部8个值。

每一个基本蝶形运算的两个输出值存入下一级蝶形运算图的两个输入存储单元。因此，在编写程序代码时，某一列的N个数据存入存储单元之后，进入下一列蝶形运算，同样得到N个运算结果，仍可存储在同一组存储器中，直至计算出最后结果，无须开设更多的存储单元，这种优化和节省内存的方式称为FFT的原位运算（inplace computation）。


采用原位运算进行优化编程仅需N个存储单元，代码设计只需输入序列x(n)(n=0，1，…，N-1)的N个存储单元，以及用于存储每一列蝶形运算支路系数WrN（r=0，1，…，N/2-1）的N/2个存储单元即可，降低了对内存资源的占用，提高了程序运行效率。

3． 蝶形运算两节点间的“距离”

DIT基2 FFT算法是以蝶形图为单元的运算结构，在由计算机实现通用算法时，需要对各节点定位，以图56所示的N=8点蝶形运算为例，序列x(n)的输入顺序是倒位序，X(k)输出顺序为自然顺序。在该蝶形图中，从左至右共有3列，m=1时为第1级（即第1列），每个蝶形图两节点间的“距离”为1； m=2时，第2级运算每个蝶形图两节点间的“距离”为2； 第3级每个蝶形图两节点间的“距离”为4。

由此可得，当N=2L时，FFT蝶形运算图输入为倒位序，输出正常顺序时，用d表示第m级运算每个蝶形图两节点间的“距离”，则d可用下式表示： 


d=2m-1


式中，m表示FFT完整蝶形运算流图中列的序号，即第m列。

4． WrN的确定

在FFT完整蝶形运算图中，第m级蝶形运算两节点间的“距离”为2m-1，因而式（523）可写成如下形式： 


Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WrN
Xm(k+2m-1)=Xm-1(k)-Xm-1(k+2m-1)WrN
（524）


观察图56，分析FFT蝶形运算流程可以发现，运算流图中的系数因子WrN，第一列仅一个因子W0N； 第二列有两个因子W0N，W2N； ……，以此类推，最后一列有N/2个WrN系数因子，分别为W0N，W1N，…，WN2-1N。

对于程序实现，通过数学推导可以得到r值的递推算法，具体步骤如下： 

（1） 将式（524）蝶形运算两节点的第一个节点标号k表示成L位二进制数； 

（2） 将二进制数左移L-m位，右边空位补零（即乘以2L-m）； 

（3） 将所得二进制数转换为十进制数，该数即为r值。

5.3.4DIT基2算法其他形式流图

图56所示的蝶形运算图描述了FFT快速算法的顺序结构和运算量，对于FFT算法流程图，若各节点的连接支路以及支路系数不变，则无论各节点和支路如何变形（改变空间位置），其算法和运算量均不会改变。即仅改变节点和支路的空间位置，蝶形图的运算结果不会改变，但由于支路物理位置改变，输入数据的提取方式（顺序）和输出的顺序可能会发生变化。以图56为例，在不改变支路连接关系和支路乘数因子（旋转因子）的条件下，可以得到多种其他形式的流图。

1． 输入自然顺序、输出倒位序

对于图56所示蝶形运算流图，按如下两个步骤改变支路的物理位置（不改变节点的链接支路）。

（1） 将第2根水平线和第5根水平线平移互换位置，即x(4)水平相连的所有节点和x(1)水平相连的所有节点（包括输入数据和输出数据节点）互换位置。

（2） 将第4根水平线和第7根水平线平移互换位置，即x(6)水平相连的所有节点和x(3)水平相连的所有节点互换位置，该互换也包括输入数据和输出数据节点。

图56所示蝶形运算图经过上述两个步骤的变化，可得如图510所示的蝶形运算流图。根据碟形运算的规律，图510与图56的蝶形图运算结构完全相同，运算量也完全相同。不同点主要体现在如下两方面： 



图510DIT基2 FFT算法输入自然顺序、输出倒位序的FFT流图


（1） 图510和图56的输入顺序不一样，图56输入为倒位序，输出为自然顺序，而图510输入为自然顺序，输出却是倒位序。

（2） 支路乘数因子的顺序不同，图56的最后一列从上到下，使用系数的顺序为W0N，W1N，W2N，W3N，前一列系数是后一列系数中的偶数幂系数（如W0N，W2N）； 而图510的最后一列使用系数的顺序为W0N，W2N，W1N，W3N，前一列系数是后一列系数的前一半，即仅使用前一半系数W0N，W2N，顺序为W0N，W0N，W2N，W2N。该运算流程图就是库利和图基最初得出的按时间抽取法算法流程图。


2． 其他形式的流图

如果在不改变图510各节点的支路连接数量和连接方式的前提下，继续改变某些节点和支路的空间位置，则可以得到输入和输出均为正常顺序的FFT蝶形运算图，如图511所示。该算法流程图的优点是输入和输出均为正常顺序，无须按倒位序排序。缺点是程序设计时不能进行原位运算，N个输入数据必须开设2N个以上的复数存储单元。



图511输入和输出均为正常顺序的FFT流图


FFT蝶形运算图除了图510和图511所示的变化之外，还可以在保持各节点连接支路的数量和连接方式不变的前提下，再做进一步的变化，即图56和图510，都可以进一步演变为其他形式，得到其他形式的FFT蝶形运算图。

【例51】用蝶形图计算FFT示例。

已知序列x(n)如下： 


x(n)=［1，2，1，1］


根据DIT基2快速算法原理，给出基于序列x(n)的DIT基2快速运算蝶形图，并对序列x(n)进行快速DFT计算（即FFT计算）。

解： 
序列x(n)为4点有限长序列，根据DIT基2算法原理，其完整蝶形图如图512所示。



图512DIT基2（N=4）FFT蝶形运算图


N=4=22，共有两级蝶形运算（见图512），支路系数（旋转因子）的数值已计算并标注在图512中，输入序列x(n)的数据为倒位序。

虚线左边为第一级蝶形运算，从上往下四个节点的计算结果如下： 


1+1=2，1-1=0，2+1=3，2-1=1


虚线右边为第二级蝶形运算，第二级运算结果如下： 


X(0)=2+3=5,X(1)=0-j=-j
X(2)=2-3=-1,X(3)=0+j=j


计算过程表明，4点蝶形运算中，由于W04=1,W14=-j，4点FFT实际上没有乘法运算，非常节省运算量。


5．4按频率抽取（DIF）基2 FFT算法
◆

在DIT基2算法中，对序列x(n)的分组采用按序号n的奇偶分解序列，得到FFT算法。如果对x(n)不按奇偶分组，而是直接将序列的前一半和后一半分为两组，得到两个新的子序列，然后对新的子序列继续均分，一直分解到所有子序列为2点为止，这样就得到了一种新的FFT快速算法，这就是按频率抽取FFT算法，即DIF基2 FFT算法。

该方法对x(n)按序号的中间位置对半均分，实际上形成了输出X(k)按频域序号k的奇、偶分组排序。该算法由桑德图基首先提出和完善，故又称为桑德图基算法。

5.4.1DIF算法原理

设序列x(n)的点数为N=2L，L为正整数，先将序列x(n)按前一半、后一半分为两个子序列。因此，DFT公式可进行如下形式的变化:


X(k)=∑N-1n=0x(n)WnkN=∑N2-1n=0x(n)WnkN+∑N-1n=N2x(n)WnkN
=∑N2-1n=0x(n)WnkN+∑N2-1n=0xn+N2Wn+N2kN
=∑N2-1n=0x(n)+xn+N2WNk/2NWnkN,k=0,1,…,N-1
（525）



由于WNk2N=（-1）k，代入式（525），可得


X(k)=∑N2-1n=0x(n)+(-1)kxn+N2WnkN,k=0,1,…,N-1（526）



若k为偶数，(-1)k=1； 若k为奇数，(-1)k=-1，因此，根据频域序号k的奇偶情况可将式（526）表示为


X(2r)=∑N2-1n=0x(n)+xn+N2W2nrN
=∑N2-1n=0x(n)+xn+N2WnrN/2（527)
X(2r+1)=∑N2-1n=0x(n)-xn+N2Wn(2r+1)N
=∑N2-1n=0x(n)-xn+N2WnNWnrN/2
（528)


令


x1(n)=x(n)+xn+N2
x2(n)=x(n)-xn+N2WnNn=0,1,…,N2-1
（529)



即x(n)前半区间某一序列值与间距为N/2的后半区间的对应序列值之和组成了新的序列x1(n)，x(n)前半区间某一序列值与间距为N/2的后半区间对应的序列值之差再与WnN相乘组成了新的序列x2(n)。新序列x1(n)和x2(n)均为N/2点序列。

于是，式 (527)和式（528）可简化为


X(2r)=∑N2-1n=0x1(n)WnrN/2
X(2r+1)=∑N2-1n=0x2(n)WnrN/2
r=0,1,…,N2-1
（530）


式(530)中的第一式是关于序列x1(n)的N/2点DFT运算，第二式为关于序列x2(n)的N/2点DFT运算。

式（529）所表示的运算关系可以用如图513所示的蝶形运算图表示。



图513DIF基本蝶形运算图


这是DIF基本蝶形运算图，与图51进行对比可以看出，它与DIT蝶形图类似，DIF的每一个基本蝶形运算仅包括一次复数乘法及两次复数加法。

根据图513，若N=23=8时，以n=0，1，2，3代入，可以得出DIF基2第一次分解的运算过程，结合式（530）可得第一次分解DIF蝶形运算如图514所示，N点DFT分解为左边的4个基本蝶形图和两个N/2点的DFT。第一次分解完成之后，其计算过程是从左往右按列计算，即第一列按4个蝶形图计算，右边DFT方框依然按DFT公式计算。



图514N点DFT分解为两个N/2点DFT（DIF N=8）


与DIT算法的推导过程类似，第一次分解以后，图514的左边包含了4个基本蝶形运算，即x(0)与x(4)、x(1)与x(5)、x(2)与x(6)、x(3)与x(7)组成了4个基本蝶形运算，每一个蝶形运算包含一次复数乘法。

至此，一个N点DFT分解成了两个N/2点DFT，若直接计算N/2点DFT，则每一个N/2点DFT的复数乘法次数如下： 


N22=N24


复数加法次数如下： 


N2N2-1


因此，图514中两个N/2点DFT共需N22次复数乘法和NN2-1次复数加法。此外，左边第一列N/2个蝶形运算共有N/2次复数乘法和2×N/2=N次复数加法。因而第一次分解全部完成之后，复数总运算量如下。

复数乘法次数： 


N22+N2+N（N+1）2≈N22


复数加法次数： 


NN2-1+N=N22


因此，DIF经过第一次分解，其运算量与DIT的第一次分解的运算量完全一样，约为直接采用DFT运算量的一半。

既然长序列分解为短序列可以降低运算量，那么接下来继续进行第二次分解，与DIT的分解类似，由于N=2L，故N/2仍是偶数，对x1(n)和x2(n)两个N/2点的序列按前一半后一半继续分解，得到两个N/4点的子序列。以x1(n)为例，按序号从中间均分为两组，每组N/4点。

于是，一个N/2点DFT被分解为两个N/4点DFT。当N=8时，分解过程如图515所示。



图515N点DFT分解为4个N/4 点DFT（DIF N=8）


对于N=2L序列，分解共需进行L次，即直到序列被分解为每一个子序列均为2点为止。2点DFT正好组成一个基本蝶形运算，而且乘数因子W0N＝１，无须乘法运算。当N=8时，DIF基2 FFT完整蝶形运算图如图516所示。



图516DIF基2 FFT算法流图（N=8）


图516称为DIF基2 FFT完整蝶形运算图（N=8时），与DIT快速算法类似，由x(n）计算X(k）是从左至右按列进行计算，因此，8点FFT共需要3级计算。

第1级运算，按左边第一列4个蝶形图进行计算，得到8个中间值； 第2级运算，以第1级的结果为输入，按第二列4个蝶形图计算，得到8个中间值； 以此类推进行第3级运算，直至完成最后一级运算，得出全部X（k）,k=0,1,2,…,N-1。

若N=8，上述3级计算中，每一级计算均包含N2=4次复数乘法运算，即每一个基本蝶形运算包含一个乘法，因此，依据图516所示的FFT蝶形运算计算DFT，仅需12次复数乘法运算（包括乘W0N=1）。

上述3级蝶形运算中，每一级基本蝶形图的系数也具有如下规律。

（1） 第1级蝶形运算的系数为W0N，W1N，W2N，W3N。

（2） 第2级蝶形运算的系数为W0N，W28。

（3） 第3级蝶形运算的系数均为W0N。

以此类推，也可以得出N=2L，DIF基2 FFT快速算法完整的蝶形运算流图以及FFT的L级运算过程。

5.4.2DIF基2算法的特点

DIF基2算法具有与DIT算法类似的特点，如根据序列分解过程及原理，可以很方便地得出N=4点和N=8点FFT流程图的递推关系，进而可得出N=2L-1点与N=2L点FFT算法流图之间的递推关系，以及DIF基2 FFT算法也可实现原位运算等。

1． DIF的原位运算

从图516可以看出，DIF与DIT的FFT算法结构相似，也具有明显的规律性，当N=8时，从左至右，共包括三级蝶形运算。若N=2L，则包括L级（列）运算，每级都包含N/2个基本蝶形图。每一个蝶形结构完成如下基本迭代运算： 


Xm(k)=Xm-1(k)+Xm-1(j)
Xm(j)=［Xm-1(k)-Xm-1(j)］WrN
（531）


式中，m表示DIF基2完整蝶形运算的第m列，k、j分别表示第k行和第j行，式（531）可以用如图517所示的蝶形图表示，共包括一次复数乘法和两次复数加法。



图517DIF蝶形运算结构


在编程实现DIF的FFT算法时，同样可以实现原位运算，即仅需N个存储单元和N/2个乘数因子存储单元即可实现从输入到输出的运算过程。

2． 两节点间的“距离”

与DIT基2 FFT算法类似，编程实现算法时，也需要对各节点进行定位，从图516可以看出，其输入数据为正常顺序，输出为倒位序。该蝶形图从左至右共有3级，m=1时为第1级（即第1列），每个蝶形图两节点间的“距离”为4； m=2时为第2级（第2列），每个蝶形图两节点间的“距离”为2； 第3级每个蝶形图两节点间的“距离”为1。

若N=2L，用d表示第m级运算蝶形图两节点间的“距离”，则d可用下式表示： 

d=2L-m=N2m（532）


式中，N=2L，m表示FFT完整蝶形运算流图中列的顺序，即第m列。

3． WrN的计算

对DIF的完整蝶形图进行第m级运算时，基本蝶形运算两节点间的“距离”为2L-m，因而第m级的蝶形计算可表示为


Xm(k)=Xm-1(k)+Xm-1k+N2m
Xmk+N2m=Xm-1(k)-Xm-1k+N2mWrN

（533）


用代码实现DIF基2 算法时，对于m=0，1，…，L的不同级运算，可以得出支路乘数因子WrN中指数r的递推算法，具体步骤如下： 

（1） 将式（533）蝶形运算两节点中的第一个节点标号k表示为L位二进制数； 

（2） 该二进制数左移（m-1）位，即乘以2m-1，移位时右边空位补零； 

（3） 将所得二进制数转为十进制数，即为r值。

5.4.3DIT与DIF的比较

1. DIT与DIF的不同点

（1） DIT基2 FFT的输入序列是倒位序，输出是自然顺序，DIF正好相反，但这不是绝对的，无论是DIT还是DIF，都可以通过改变蝶形运算图支路的空间位置实现输入或输出的重新排列，两者的输入或输出顺序都可以改变为自然顺序或倒位序。

（2） DIT的基本蝶形(见图51)与DIF的基本蝶形(见图517)不同，DIT是先乘法运算，再进行代数和运算，而DIF则是先求代数和，再乘法运算。这是DIT与DIF本质的不同点。

2． DIT与DIF的相同点

（1） 两者的运算量相同，即都有L级（列）运算，每一级运算均包括N/2个基本蝶形运算，总共需要mT=N2log2N次复数乘法和aT=Nlog2N次复数加法。

（2） DIT基2 FFT和DIF基2 FFT都可通过编程实现原位运算。

（3） DIT与DIF的基本蝶形运算流图互为转置关系。若将DIT的基本蝶形运算流图进行转置，就得到DIF的基本蝶形图，反之亦然。转置是指将流图的所有支路都改变方向，输入和输出互换。经转置运算，由图51可以得到图517，反之也可以由图517得到图51。

由此可知，根据输入和输出是按自然顺序还是按倒位序排列的不同，有四种原位运算的FFT算法流图，DIT与DIF各两种，DIT与DIF算法的特点如表54所示。


表54DIT与DIF算法的特点（N=2L）



按时间抽取（DIT）

输入自然顺序、输出倒位序
输入倒位序、输出自然顺序




m级蝶形图两节点间的距离
2L-m=N2m
2m-1
第m级蝶形运算公式
Xm(k)=Xm-1(k)+Xm-1k+N2mWrN

Xmk+N2m=Xm-1(k)-Xm-1k+N2mWrN
Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WrN

Xm(k+2m-1)=Xm-1(k)-Xm-1(k+2m-1)WrN

r值计算（WrN）
1． k表示为L位二进制数

2． 右移L-m位，左边空出位补0

3． 补0后倒位序排列即为r值
1． k表示为L位二进制数

2． 左移L-m位，右边空出位补0

3． 补0后的结果即为r值


按频率抽取（DIF）

输入自然顺序、输出倒位序
输入倒位序，输出自然顺序




m级蝶形图两节点间的距离
2L-m=N2m
2m-1


续表



按频率抽取（DIF）

输入自然顺序、输出倒位序
输入倒位序，输出自然顺序




第m级蝶形运算公式
Xm(k)=Xm-1(k)+Xm-1k+N2mWrN

Xmk+N2m=

Xm-1(k)-Xm-1k+N2mWrN
Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WrN

Xm(k+2m-1)=［Xm-1(k)-Xm-1(k+2m-1)］WrN

r值计算（WrN）
1． k表示为L位二进制数

2． 左移m-1位，右边空出位补0

3． 补0后的结果即为r值
1． k表示为L位二进制数

2． 右移m-1位，左边空出位补0

3． 补0后倒位序排列即为r值



5.4.4IDFT快速算法

离散傅里叶逆变换的快速运算称为IFFT。DFT和IDFT的算法结构类似，差异非常小。

DFT正变换:


X(k)=∑N-1n=0x(n)WnkN,k=0,1,…,N-1（534）


DFT逆变换:


x(n)=IDFT［X(k)］=1N∑N-1k=0X(k)W-nkN,k=0,1,…,N-1（535）


正变换和逆变换的差别主要表现在两方面，一是正逆变换中的WN指数的符号相反； 二是逆变换算法包含常数因子1/N，而正变换的因子为1。

由于正变换和逆变换算法的内核要素相同，因此，无论是DIT基2 FFT还是DIF基2 FFT算法，若解决了逆变换与正变换的两个不同点，则正变换的快速算法也能用于逆变换。

解决上述指数符号和系数的差异问题，有如下两种方法： 

（1） 研究FFT快速算法运算流图，解决指数符号问题，实现IFFT算法； 

（2） 研究DFT和IDFT公式，利用复数性质改逆变换公式，实现IFFT运算。

1． 基于运算流图的IFFT

（1） 指数符号问题。

由于正变换和逆变换的算法结构基本相同，以DIT基2 FFT算法为例，根据DIT基2快速算法原理，只需将正变换快速算法蝶形运算中的支路增益系数WnkN置换为W-nkN即可。类似地，对于DIF基2快速算法，对其每一级蝶形做同样的改进之后也可以用于IFFT运算。

（2） 系数1/N的处理。

在乘数因子WnkN解决之后，可以将所得的运算结果均乘以常数1/N，或者在每一级蝶形运算的各支路上增加支路系数1/2，则可以解决逆变换系数1/N的问题。

也就是说，对正变换的FFT运算进行上述改进之后，按时间抽取和按频率抽取的FFT算法均可用于逆变换的计算。在运算中注意将按时间抽取的流图（见图56）或者按频率抽取的流图（见图516）中的支路乘数因子WnkN改变为W-nkN，并在每列（级）运算中乘以1/2，12L=1N。以N=8点的DIT基2 FFT为例，IFFT算法流图如图518所示。



图518DIT逆变换快速算法完整蝶形图（N=8）


2． 基于公式改进的IFFT

根据上述IFFT算法的实现过程可知，对正变换FFT的代码进行相应的修改即可用于逆变换的实现。这种方法虽然工作量小、实现方便，但需要修改FFT正变换程序的内部代码，即修改WnkN的指数符号以及对每一级运算乘以1/2才能实现逆变换的快速算法。

显然，修改既有程序内部的代码对程序员的要求较高，稍有不慎，容易引起其他问题。如果不修改DFT正变换的内部代码，直接将正变换程序用于逆变换运算，则可以从DFT的公式入手，通过对变换公式进行适当变换，利用复数的性质，使逆变换与正变换的因子WnkN的符号相同，这样就无须修改FFT正变换程序的内部代码，也可以直接将正变换程序用于逆变换计算。

先对DFT逆变换进行共轭运算


x*(n)=1N∑N-1k=0X*(k)WnkN


再次进行共轭运算


x(n)=1N∑N-1k=0X*(k)WnkN*=1N{DFT［X*(k)］}*（536）


此时，式（536）的WnkN指数的符号与正变换完全一致，因此，若先对X(k）进行共轭运算，就可以直接应用正变换快速算法的程序进行逆变换计算，然后将变换结果再进行一次共轭运算，并乘以1/N，就完成了逆变换的快速运算。

该方法的优点是无须对正变换FFT程序的内部代码进行修改，只需运算前先对逆变换的输入数据进行共轭运算，再将计算结果进行一次共轭运算并除以N，就完成了求逆变换x(n)的全部计算。


5.5基4 FFT算法
◆


基2 FFT 算法的基本方法是先将N点序列一分为二，一个N点DFT被分解为两个N/2点DFT的组合，再将N/2点DFT又一分为二，分解为4个N/4点的DFT，直至全部分解为2点DFT为止。基4 FFT的思路与基2 FFT的思路基本类似，不同点是前者将N点序列一分为四，再对各子序列继续一分为四，直至分解到全部子序列为4点为止。

5.5.1基4 FFT原理

根据基2 FFT算法原理，对于基4 FFT算法，序列x(n)的点数应满足N=4L。基4 FFT算法的基本思路是将N点序列一分为四，即一个N点DFT分解为四个N/4点DFT的线性组合，再将每一个N/4点DFT又一分为四，分解为4个N/16点的DFT，直至全部分解为4点DFT为止。而对于4点DFT运算，则采用基2 FFT运算，并以此作为基4 FFT最后一次分解的基本运算单元。因此，类似于DIT基2的分解过程，DIT基4快速算法先将长度为N=4L的序列x(n)，按如下方式分为四个子序列： 


x0(r)=x(4r)
x1(r)=x(4r+1)
x2(r)=x(4r+2)
x3(r)=x(4r+3)，r=0,1,…,N4-1
（537）


于是，DFT公式可以改变为如下形式:


X(k)=DFT［x(n)］=∑N-1n=0x(n)WnkN
=∑N-1n=0((n))4=0x(n)WnkN+∑N-1n=0((n))4=1x(n)WnkN+
∑N-1n=0((n))4=2x(n)WnkN+∑N-1n=0((n))4=3x(n)WnkN
=∑N4-1r=0x(4r)W4rkN+∑N4-1r=0x(4r+1)W（4r+1）kN+∑N4-1r=0x(4r+2)W（4r+2）kN+
∑N4-1r=0x(4r+3)W（4r+3）kN
=∑N4-1r=0x0(r)W4rkN+∑N4-1r=0x1(r)W(4r+1)kN+∑N4-1r=0x2(r)W(4r+2)kN+∑N4-1r=0x3(r)W(4r+3)kN



式中，((n))4表示n除以4的余数。根据WN的可约性，上式可以进一步化为


X(k)=∑N4-1r=0x0(r)WrkN/4+WkN∑N4-1r=0x1(r)WrkN/4+
W2kN∑N4-1r=0x2(r)WrkN/4+W3kN∑N4-1r=0x3(r)WrkN/4
（538）


根据N/4点DFT公式，式（538）可以表示为


X(k)=X0(k)+WkNX1(k)+W2kNX2(k)+W3kNX3(k),0≤k≤N4-1（539）


式中，X0(k)、X1(k)、X2(k)、X3(k)分别为N/4点时域序列x0(r)、x1(r)、x2(r)、x3(r)对应的DFT，即


X0(k)=∑N4-1r=0x0(r)WrkN/4
，X1(k)=∑N4-1r=0x1(r)WrkN/4
X2(k)=∑N4-1r=0x2(r)WrkN/4
，X3(k)=∑N4-1r=0x3(r)WrkN/4


类似DIT基2算法的第一次分解，式（539）只能计算出X(k)前四分之一的值，X(k)其余值的计算需根据WN的周期性以及有限长序列隐含的周期性，由式（539）直接得出： 


Xk+N4=X0(k)+Wk+N4NX1(k)+W2k+N4NX2(k)+W3k+N4NX3(k)
Xk+N2=X0(k)+Wk+N2NX1(k)+W2k+N2NX2(k)+W3k+N2NX3(k)
Xk+3N4=X0(k)+Wk+3N4NX1(k)+W2k+3N4NX2(k)+W3k+3N4NX3(k)
0≤k≤N4-1


根据WN的特性，上述公式可以简化为如下形式： 


Xk+N4=X0(k)-jWkNX1(k)-W2kNX2(k)+jW3kNX3(k)（540）
Xk+N2=X0(k)-WkNX1(k)+W2kNX2(k)-W3kNX3(k)（541）
Xk+3N4=X0(k)+jWkNX1(k)-W2kNX23(k)-jW3kNX3(k)（542）
0≤k≤N4-1



由式（539）和式（540）~式（542）可得如图519所示的DIT基4快速算法基本蝶形运算图。




图519DIT基4快速算法基本蝶形图



为了直观和深入理解基4 FFT算法原理，现以N=16为例分析基4快速运算流图的规律。由于N=16=4L，即L=2，根据式（539）和式（540）~式（542）对序列x(n)进行第一次分解，分别以k=0,1,2,3代入，可得当N=16时第一次分解的蝶形运算图，如图520所示。经过第一次分解之后，输出为正常顺序，输入顺序为四进制倒位序。



图520DIT基4快速算法第一次分解蝶形图



由于N=16=42仅需2次分解，经过第一次分解之后，对于4个N4点的DFT，各子序列长度均为4点，正好是DIT基4快速运算的基本单元。根据基4算法的原则，4点长的子序列DFT采用基2 FFT运算，这是基4 FFT的最后一次分解。根据图519或式（539）和式（540）~式（542）的运算规律，这时k只有一个取值，即k=0，由于需按基2运算分组，因此需要按倒位序排列。以图520左上角第一个4点DFT为例，其运算可展开如下： 


X0（0）X0（1）X0（2）X0（3）=
W04W04W04W04
W04W24W14W34
W04W04W24W24
W04W24W34W14
x（0）x（2）x（1）x（3）


由于W04=1、W14=-j、W24=-1、W34=j，并考虑到2点DFT仍需按偶序号和奇序号排列，因此矩阵表达式可以进一步简化为如下形式： 


X0（0）X0（1）X0（2）X0（3）=
1111
1-1j-j
1-11-1
1-1j-j
x（0）x（2）x（1）x（3）
（543）




图521基4 FFT基本运算的信号流图


根据式（543），可得第一列左上角4点基4快速算法的蝶形运算过程如图521所示。


蝶形图的计算是从左往右按列依次运算，根据图521可知，基4 FFT的第一级运算具有如下特点： 

（1） 一个基4算法包含有2点DFT蝶形结构。

（2） 第一级4点基4 FFT运算完全不需要复数乘法运算，虽然流图中有支路系数-j，但实际上这是交换复数的实部和虚部，不用乘法运算也可完成。

根据第一次分解和第二次分解过程可以得出当N=16时DIT基4 FFT的完整蝶形运算流图如图522所示，序列x(n)经过两次分解之后，输出为正常顺序，输入为四进制倒位序，图中虚线框内为一个4点基4 FFT流图的基本单元。



图522DIT基4 FFT流图



5.5.2基4 FFT运算量

根据基4 FFT算法原理，对于第一级基本的4点FFT，实际上无须乘法运算。根据图522所示的算法流图，支路有乘旋转因子才包含复数乘法运算，每一个4点DFT仅有3次需乘旋转因子，另一乘W0N=1。每一级基4 FFT包含N/4个4点DFT，因此，每一级共包含3×N/4次复数乘法运算，由于N=4L，共有L级运算，L=(log2N)/2，第一级运算的因子为±1和j，实际上没有乘法运算，因此，基4 FFT算法总的复数乘法次数如下。

（1） 复数乘法次数。


mT=34N(L-1)=34N12log2N-1≈38Nlog2N,L1（544）


由此可知，基4 FFT与基2 FFT相比乘法运算量更少，运算效率更高。由于一些特定的乘数因子无须相乘，如当N=16时，W016=1、W416=-j、W816=-1等，这些乘数因子实际上无须乘法运算，故根据图522，当N=16时，基4 FFT仅需8次复数乘法运算。

对于复数加法运算量，根据运算流图522可知，求和节点和基2 FFT的数量相等，由此可得，基4 FFT所需的复数加法次数和基2 FFT完全相等。

（2） 复数加法次数。


aT=Nlog2N（545）



5.6分裂基FFT算法
◆

基2 FFT算法的复数乘法运算量为N2log2N，为了进一步寻找更高效的算法，1984年，有学者在基2和基4算法的基础上提出了分裂基算法（splitradix）。根据现有公开文献和资料，若序列点数满足N=2L，分裂基FFT算法被认为是目前运算量最优的FFT算法，这是因为该算法所包含的复数乘法次数最少，并可实现原位运算。

5.6.1分裂基算法原理

分裂基算法又称为基2基4混合基算法，它既与基2算法相关，也与基4算法相关。无论是按时间抽取还是按频率抽取，在基2 FFT快速算法的计算过程中（见图52和图54），每一级运算中的奇序号（DIT观察输入序号，DIF观察输出序号）均乘以一个旋转因子，偶序号则不用乘旋转因子。由于基4 FFT比基2 FFT更高效，故分裂基算法的基本方法是对偶序号采用基-2 FFT算法，对奇序号采用基4 FFT算法，从而进一步降低复数乘法的运算量。

与DIT和DIF基2 FFT算法类似，分裂基FFT算法也要求N=2L（L为正整数）。

DFT正变换公式如下： 


X（k）=∑N-1n=0x(n)WnkN,0≤k≤N-1


根据分裂基算法思想，一个N点的DFT 应分解为一个N/2点的DFT和两个N/4点的DFT。这种分解既有基2的部分，又有基4的部分，故称为分裂基分解。因此，时域序列x(n)按序号n分解为如下三个子序列:


x1(r)=x(2r),0≤r≤N2-1
x2(l)=x(4r+1)
x3(l)=x(4l+3),0≤r≤N4-1


序列分解以后，DFT公式可以改变为如下形式:


X(k)=∑N2-1r=0x(2r)W2rkN+∑N4-1r=0x(4r+1)W(4r+1)kN+∑N4-1r=0x(4r+3)W(4r+3)kN
=∑N2-1r=0x1(r)WrkN/2+WkN∑N4-1r=0x2(r)WrkN/4+W3kN∑N4-1r=0x3(l)WrkN/4
=X1(k)+WkNX2(k)+W3kNX3(k)
（546）


式中


X1(k)=∑N2-1r=0x1(r)WrkN/2=∑N2-1r=0x(2r)WrkN/2
X2(k)=∑N4-1r=0x2(r)WrkN/4=∑N4-1r=0x(4r+1)WrkN/4
X3(k)=∑N4-1r=0x3(r)WrkN/4=∑N4-1r=0x(4r+3)WrkN/4


X(k)为N点DFT，而X1(k)为x(n)中由偶序号组成的N/2点DFT； X2(k)、X3(k)均为x(n)中由奇序号组成的N/4点DFT，根据WnkN、X1（k）、X2（k）和X3（k）所隐含的周期特性可得


X1(k)=X1k+N2（547）
X2(k)=X2k+N4（548）
X3(k)=X3k+3N4（549）



将X(k)分为如下四个区间:



X(k)=X1(k)+WkNX2(k)+W3kNX3(k)
Xk+N4=X1k+N4-j
WkNX2(k)+jW3kNX3(k)
Xk+N2=X1(k)-WkNX2(k)-W3kNX3(k)
Xk+34N=X1k+N4+jWkNX2(k)-jW3kNX3(k)
（550）
0≤k≤N4-1



式（550）的运算关系可表示为如图523所示的分裂基蝶形运算图。



图523DIT分裂基FFT算法的第一级运算流图


由此可得分裂基FFT快速算法的基本蝶形运算单位可表示为如图524所示的形式。



图524分裂基FFT算法的基本蝶形运算


与基2算法类似，以N=16为例，基于式（550）可得出DIT分裂基FFT算法第一次分解的蝶形运算过程，如图525所示。



图525DIT分裂基（N=16）FFT算法第一级分解过程


从图525可以看出，当N=16时，第一级分解得到了图中虚线框内的4个分裂基、虚线框左方的8点DFT（左上）和两个4点DFT（左下）。

在第一级分解的基础上，根据式（550）继续进行下一级分解，即对虚线框左上的8点DFT以及左下的两个4点DFT继续进行分解，直至分解到分裂基的最小单元。对虚线框左上方的8点DFT，根据式（550）或图523可知，对X1(k)进行第二级分解包含2个分裂基以及一个4点DFT和两个2点DFT，而X2(k)和X3(k)均为4点DFT，每一个4点DFT正好包含一个分裂基单元和一个2点DFT。

由此可得出当N=16时，完整分裂基运算流程图共包含9个分裂基，如图526所示。由于图526分裂基FFT算法是根据时域序号n的奇偶对序列进行分解的，因此，其输入和输出顺序与DIT基2 FFT算法一样，输入为正常顺序，输出为倒位序。与基2快速算法类似，分裂基快速算法也包括按频率抽取的分裂基快速算法。



图526DIT分裂基（N=16）FFT算法完整蝶形运算图


5.6.2分裂基算法的运算量

对于基2和基4算法，无论是基于时间抽取还是基于频率抽取，其快速算法的运算量均与算法的基本蝶形单元数量有关。分裂基算法的运算量也类似，对于N=2L点的分裂基算法的运算量，其复数乘法次数与基本分裂基单元的数量密切相关，每一个分裂基蝶形结构包含两次复数乘法运算，而复数加法的运算量与蝶形图节点总数相关。若序列的点数相同，则基2、基4和分裂基算法的节点数相同，故复数加法运算量与基2算法相同，不同FFT算法的运算量的差别主要体现在复数乘法运算量上，因此，估计分裂基算法总运量实际上就是计算基本分裂基单元的数量。

根据第一次分解后的蝶形运算流图（见图525）可知，当N=2L=16时，第一级分解得到了2L-2=22个分裂基，同时还有2L-1=8点的DFT和两个2L-2=4点的DFT。

设N=2L时，分裂基的数量为sL，因此，可得如下关于分裂基数量SL的递推方程： 


sL=2L-2+sL-1+2sL-2


由于4点构成一个基本分裂基单元，由此可得完整的递推关系表达式为


sL=2L-2+sL-1+2sL-2,L≥2（551）


初始条件为


s0=s1=0,s2=1



因此，当N=2L时，采用分裂基算法的复数乘法次数为


mT=2sL


即


mT=2L-1+2sL-1+4sL-2,L≥2（552）


若初始条件为


s0=s1=0,s2=1


根据式（551）及式（552），可得4点以上分裂基算法复数乘法运算量如表55所示，由于支路系数为W0N、WN/2N时并不需要进行乘法运算，因此，实际上乘法运算量比表55略少。


表55分裂基FFT算法复数乘法次数



N

481632641282565121024
L
2345678910
SL
13923571353137111593

mF
24184611427062614223186



表55可以查询1024点以下的分裂基算法的复数乘法运算量，对于N>1024点的序列x(n)，可以由式（552）递推求得复数乘法次数。从表55可以发现，分裂基算法的复数乘法次数比13Nlog2N更少。以N=128为例，13Nlog2N=299，而mT=270<299； 基2快速算法的乘法次数为N2Nlog2N=448，基4算法的复数乘法次数为38Nlog2N=336，均比分裂基算法的运算量大一些。


5.7线性调频z变换
◆

FFT算法可以快速计算出序列x(n)的离散傅里叶变换，即有限长序列x(n)的z变换X(z)在z平面单位圆上N个等间隔抽样值。尽管z变换和DFT应用非常广泛，对推动离散频域分析起了非常重要的作用，但z变换在面对如下问题时仍存在明显的不足。

（1） 工程应用中有时可能仅对信号的某一频段感兴趣，即只需要计算单位圆上某一局部频段的频谱值，例如，为提高窄带信号的频率分辨率，应在窄带范围内进行高密度抽样，而在窄带范围之外可以不抽样。由于DFT是基于单位圆上的等间隔抽样，故为了满足窄带部分需求应大幅增加频域的抽样点数，但这样做不仅增大了计算量； 进而影响实时性，而且浪费资源。

（2） 对非单位圆上进行抽样，例如语音信号处理需要了解z变换极点处的复频率，而对于稳定系统，极点位置可能与单位圆有一定的距离，因此，仅在单位圆上进行抽样难以了解极点处复频率，因而客观上需要在单位圆内部的某一曲线上进行抽样。

线性调频z变换（Chirpz变换，CZT）的提出正是为了满足上述需求，它是一种可沿z平面任意螺线进行抽样的频谱分析方法，并可利用FFT原理进行快速计算。

5.7.1基本原理


已知序列x(n)是点数为N的有限长序列，则x(n)的z变换为


X(n)=∑N-1n=0x(n)z-n


为满足z变量在z平面的非单位圆路径上进行抽样，可以选择z平面的一条螺旋线进行等间隔（角度）的抽样。设抽样点zk为


zk=AW-k,k=0,1,…,M-1（553）


M为所分析的复频域的抽样点数，M可以大于或等于N，A和W都是任意复数，为满足抽样路径的灵活多样性，令


A=A0ejθ0
W=W0e-j0
（554）


式中，A0和W0为正实数，将式（554）代入式（553）得


zk=A0ejθ0W-k0ejk0=A0W-k0ej(θ0+k0)（555）


把k=0,1,2,…,M-1代入式（555），可得


z0=A0ejθ0
z1=A0W-10ej(θ0+0)
︙
zk=A0W-k0ej(θ0+k0)
︙
zM-1=A0W-(M-1)0ej［θ0+(M-1)0］


抽样值zk在如图527所示的z平面的螺旋线上进行。



图527CZT在z平面螺旋线上的抽样点


根据上述分析及式（553）可知： 

（1） A0表示螺旋线的第一个抽样点（起始抽样点）到坐标原点的距离，即抽样点z0的向量半径长度，由于稳定系统的极点在单位圆内，因此A0≤1。

（2） θ0表示第一个抽样点z0的相位角，可以为任意正负值。

（3） 0表示两相邻抽样点之间的相角差，若0>0时，表示抽样沿逆时针方向进行； 0<0，表示抽样沿顺时针方向进行。


（4） W0表示螺旋线的伸展率，W0>1时，螺旋线随着k的增加向圆心收缩； W0<1时，螺旋线随k的增加远离圆心向外伸展； W0=1表示以A0为半径的一段圆弧，若A0=1则该圆弧是单位圆的一段。

若0=2πN，M=N,A=A0ejθ0=1,W0=1，则可得出W=e-j2πN，螺旋线演变为单位圆，式（553）演变为N点DFT运算，即N个抽样点等间隔均匀分布在单位圆上。这时若取0为非零任意值，则DFT的结果是单位圆上某一段圆弧的频谱。

将式（553）代入z变换表达式，可得


X(zk)=∑N-1n=0x(n)z-nzk=∑N-1n=0x(n)A-nWnk,0≤k≤M-1（556）


5.7.2CZT快速算法

1． 快速算法原理

若对式（556）直接进行计算，则计算过程和方法均与直接计算DFT类似，可获得频域M个抽样点的离散频谱。

式（556）的运算量也可以准确估计，共包含NM次复数乘法和（N-1）M次复数加法。显然，若N与M大致相当，则乘法次数（NM）与DFT的N2为同阶量级，若采样点数很大，即N与M均较大，则运算量很大，将影响到CZT的运算速度，因此，应考虑采用快速运算。

布鲁斯坦（Bluestein）提出了实现CZT快速算法的基本思路，可以将该运算转换为卷积运算形式，从而可以根据FFT算法原理实现CZT快速算法，提高运算速度。

布鲁斯坦快速算法原理利用了如下代数恒等式： 


nk=12［n2+k2-(k-n)2］


将该恒等式代入式（556），可得


X(zk)=∑N-1n=0x(n)A-nWn22W-(k-n)22Wk22
=Wk22∑N-1n=0x(n)A-nWn22W-(k-n)22

（557）


为利用快速卷积对式（557）进行计算，设


g（n）=x(n)A-nWn22，n=0,1,…,N-1（558）
h(n)=W-n22（559）


则式（557）可得到以下类似卷积形式的表达式： 


X(zk)=Wk22∑N-1n=0g(n)h(k-n),k=0,1,…,M-1（560）


该式表明，抽样点zk处的z变换可通过求g(n)与h(n)的卷积结果再乘以系数Wk22实现。因此，式（560）可表示为


X(zk)=Wk22［g(k)*h(k)］,k=0,1,…,M-1（561）


由此可得，计算X(zk)的流程如图528所示。



图528Chirpz变换运算流程


序列h(n)=W-n22可以理解为频率随离散时间（n）呈线性增长的复指数序列，该信号形式是雷达系统中的一种重要信号类型，称为线性调频信号(Chirp signal)，因此，该变换称为线性调频z变换。

2． 快速算法实现步骤

根据式（561）计算卷积，先分析序列g(n)和h(n)的长度和区间范围，由于x(n)是N点序列，n=0,1,2,…,N-1，因此g(n)也是Ｎ点序列。而h(n)=W-n22为无限长偶对称序列，但由于z平面上的抽样点只有M个，根据式（560），若将h(n)视为单位抽样响应，则


k=0,1,2,…,M-1


h(n)是非因果线性系统，h(n)在n=-(N-1)~(M-1)范围内取值，如图529所示，h(n)为 L=N+M-1点有限长序列。



图529序列h(n)的区间范围


根据序列自变量的取值范围，g(n)*h(n)卷积序列的长度为2N+M-2，根据圆周卷积代替线性卷积的条件，圆周卷积长度应不小于2N+M-2。

由于仅需计算X(zk)的前M个值，因此，为了降低运算量，可将圆周卷积的点数L取为


L≥N+M-1


L应为2的整数幂，取满足条件的最小L值即可。根据运算要求，对序列h(n)应补零值点，使序列长度等于L，即从n=M开始补L-（N+M-1）个零点，当n=L-N时，对序列h(n)以L为周期进行周期延拓，对g(n)同样在末端补零至L点。h(n) 和g(n)均为L点序列之后，就可以进行CZT快速运算了，一般包括如下五个步骤。

（1） 计算FFT［g(n)］。

根据线性调频计算原理、圆周卷积代替线性卷积的条件以及FFT计算对点数的要求，选择满足L≥N+M-1和L取2的整数幂的最小整数，然后对序列g（n）=x(n)A-nWn22的末端按如下方式补零值点，使序列g(n)成为L点有限长序列:


g（n）=A-nWn22x(n),0≤n≤N-1
0,N≤n≤L-1
（562）


并计算序列g(n)的L点FFT运算，得L点G(k)，即


G(k)=FFT［g(n)］


（2） 计算FFT［h(n)］。

对序列h(n)补零点成为L点有限长序列，在n=0~（M-1）范围内取h(n)=W-n22，在n=M~（L-N）范围内h(n)可取任意值，或者直接取零值，在n=（L-N+1）~（L-1）范围内，h(n)=W-n22并进行周期延拓，则延拓序列为W-m22，即


h(n)=
W-n22,0≤n≤M-1
0(或任意值),M≤n≤L-N
W-（L-n）22,L-N+1≤n≤L-1
（563）


根据序列h(n)的表达式可知，h(n)是序列W-m22以L为周期进行周期延拓的主值序列，如图530所示，对序列h(n)进行L点FFT运算，即


H（k）=FFT［h（n）］





图530序列h(n)延拓后的主值序列


（3） 计算频域乘积。

计算频域序列H(k)和G(k)的乘积，即计算Y(k)=H(k)G(k)，得到L点频域序列Y(k)。

（4） 计算圆周卷积。

计算Y(k)的L点傅里叶逆变换，即对Y(k)进行IFFT计算，得到h(n)与g(n)的圆周卷积序列为


y(n)=h(n)g(n)=1L∑L-1k=0H(k)G(k)ej2πLkn


圆周卷积序列y(n)的前M个值等于原h(n)与g(n)的线性卷积，当n≥M时，序列y(n)的值无意义，可以不予计算。

（5） 计算X（zk）。


X（zk）=Wk22y(k),0≤k≤M-1（564）


3． MATLAB实现CZT快速算法

MATLAB提供了实现CZT快速算法函数czt，调用格式如下： 


y=czt(x，M，W，A)


该函数的功能是计算长度为M点的序列x(n)的Chirpz变换X（zk），A为频率抽样的起始点位置，W为采样轮廓线上各抽样点之间的比率，y用于存储CZT的变换值X（zk）。

5.7.3CZT算法运算量

当N和M很大时，采用CZT快速算法计算X（zk）比直接求X（zk）能有效节省运算量，CZT快速算法所需的运算量主要包括如下几项。

（1） 求g(n)的系数。补零点时获得L点序列g(n)=(A-nWn22)x(n)，由于仅有N点非零值，复数乘法次数为N次，而g(n)的系数A-nWn22可用如下递推法求解： 


g(n)=(A-nWn22)x(n)=Cnx(n)


式中


Cn=A-nWn22=Cn-1Dn-1
Dn=WnW12A-1=WnD0=WDn-1


初始条件为


C0=1,D0=W12A-1=W1/20A0e-j02+θ0


因此，根据式（585）可求出N个系数Cn，共包括2N次复数乘法。

（2） 求h(n)的系数。补零点得到L点序列h(n)，而h(n)由W-n22在-(N-1)≤n≤M-1段内的序列值构成，由于W-n22具有偶对称性，故若N>M，则只需求出0≤n≤N-1范围内的N点序列值即可。和步骤（1）类似，W-n22的值同样采用递推法求解，共需2N次复数乘法。

（3） 计算G(k)H(k)。需两次L点FFT运算，复数乘法运算量为2×12Llog2L。

（4） 计算圆周卷积y(n)。需一次L点IFFT运算，复数乘法运算量为12Llog2L。

（5） 计算Y(k)=G(k)H(k)。需L次复数乘法。

（6） 计算X(zk)。X(zk)=Wk22y(k)(0≤k≤M-1)，需M次复数乘法。

因此，计算CZT全部值所需的复数乘法次数为


mT=32Llog2L+5N+L+M（565)


而直接计算X(zk)的复数乘法次数为NM次，当N、M均为较大数时，采用基于FFT的CZT快速算法比直接计算的运算量要小很多，且N、M越大，优势越明显，例如若N=M=64，CZT快速算法与直接算法的乘法运算量相差无几，若N=M=128，则快速运算开始显示出优势。

CZT的算法非常灵活，时域输入序列x(n)的点数N可以和频域序列的点数M不相等，且N和M均可以是质数或其他任意正整数。抽样间隔0可以根据需要任意选取，即频率分辨率可以根据需要进行调整。计算z变换的螺旋线也可以根据需要进行选择。第一个采样点的选取也无特殊限制，可以根据应用需要自由选择，即可以从任意复频率开始对输入数据进行分析，可以方便地对窄带高频信号进行高密度抽样分析。相比DFT而言，CZT是进行频谱细化分析的一种方法，对于某些特定应用，FFT不能精确反映信号的局部频谱特性，但采用CZT算法可以获得更精确的频谱特性。


5.8线性卷积与线性相关的FFT算法
◆

有限长序列傅里叶变换（DFT）的快速算法是很多其他快速算法的基础，库利图基发表的傅里叶变换快速算法的论文，极大地推动了傅里叶变换走向工程应用。实际上，FFT不仅是一种快速算法，而且是一种思考问题的方法，它不限于DFT快速运算的实现，将FFT的思维方法应用于其他算法的快速实现也是FFT思想的核心。

5.8.1线性卷积的FFT算法

线性卷积运算是信号处理的基本运算，应用非常广泛。例如，对于FIR滤波器，其输出等于滤波器的单位冲激响应h(n)与输入信号x(n)的线性卷积。

设输入序列x(n)为N点，单位冲激响应序列h(n)为M点，则输出y(n)为


y(n)=x(n)*h(n)=∑M-1m=0h(m)x(n-m)


y(n)为L=N+M-1点有限长序列。

1． 线性卷积的运算量

根据线性卷积的计算公式，序列x(n)的每一个值都必须和序列h(n)的全部值进行乘法运算，因此，线性卷积的乘法的总运算量为NM次，以MCon表示乘法次数，则有


MCon=NM（566）


对于线性相位FIR滤波器而言，由于满足


h(n)=± h(M-1-n)


所以线性相位数字滤波器的网络结构通常是h(n)与h(M-1-n)共用放大器，因此，乘法运算次数可以减少一半，即


MCon=NM2（567）


由上可知，线性卷积算法复数乘法的运算量为MCon=NM，虽然在线性相位数字滤波器中运算量为MCon=NM2，当两序列长度接近时，如N与M接近或者N=M时，则运算量近似为N2。由于计算N点序列x(n)的DFT运算量为N2，即线性卷积和DFT运算量相等，都非常大，因此，对于线性卷积，采用快速算法具有很好的效果。

2． 线性卷积的快速算法

线性卷积快速算法的基本思路是以FFT算法为基础进行快速运算，即用圆周卷积来代替该线性卷积，然后用FFT算法进行计算，如图531所示。



图531线性卷积快速算法结构


当序列x(n)和h(n)末端补零值点至L≥N+M-1时，则圆周卷积可以代替线性卷积，因此，需要在序列x(n)和h(n)末端补零点，使L≥N+M-1并取L为2的整数幂，此时可以采用FFT算法，即


x(n)=x(n),0≤n≤N-1
0,N≤n≤L-1
h(n)=h(n),0≤n≤M-1
0,M≤n≤L-1


补零后计算圆周卷积为


y(n)=x(n)h(n)


由于L≥N+M-1时，圆周卷积与线性卷积结果一致。因此，可得线性卷积快速运算的基本步骤如下。

（1） 计算H(k)=FFT［h(n)］。

（2） 计算X(k)=FFT［x(n)］。

（3） 计算Y(k)=X(k)H(k)。

（4） 计算y(n)=IFFT［Y(k)］。

步骤（1）、（2）、（4）均采用FFT进行计算，步骤（3）为乘法运算，因此，线性卷积快速算法的乘法运算量包括三次FFT运算和N次乘法运算，总运算量如下： 


mT=32Llog2L+L=L1+32log2L（568）


虽然运算量包含3个L点的FFT运算，但当N和M较大时（即L比较大），依然具有很好的快速运算效果。现在以线性相位FIR数字滤波器为例，分析直接计算线性卷积和用FFT计算线性卷积这两种方法的乘法次数。

令k为mCon与mT的比值，则


k=mConmT=NM2L1+32log2L（569）


若k>1时，说明采用线性卷积快速算法能节省运算量。L=N+M-1，工程应用中，一般x(n)和h(n)的点数较大，即N和M较大，因而L≈N+M，根据代数知识，M和N越接近时，乘积NM越大，分两种情况讨论。

1） x(n)与h(n)点数接近

若M=N，对提升比值k有利，这时L=2N-1≈2N，则


k=NM2L1+32log2L≈N24L1+32log2(2N)
=N4（2.5+1.5log2N）=N10+6log2N


根据序列长度不同，运算量比值如表56所示


表56线性卷积直接计算与快速算法运算量比值（N=M）



N=M



8
32
64
128
256
512
1024
2048
4096

k



0.286
0.80
1.39
2.46
4.41
8
14.62
26.95
49.95



根据表56，若N=M=8、32时，快速算法的运算量反而大于直接卷积，此时不宜采用快速算法。若N=M=64，快速算法比直接卷积效果略好，因此N=M=64是节省运算量的临界值。若N=512，快速算法的运算量仅为直接卷积的1/8； 若N=4096，快速算法的运算量仅有直接卷积的2%。

【例52】FFT计算快速卷积示例。

已知序列x(n)和h(n)如下： 


x(n)=0.6n,0≤n≤15,h(n)=R15(n)


试计算y(n)=x(n)*h(n)。

解： 
用FFT计算线性卷积，应注意圆周卷积代替线性卷积的条件为


L≥N1+N2-1
N1=length［x(n)］，N2=length［h(n)］


计算线性卷积的MATLAB程序如下： 


clc;clear all;close all

n=［0:1:15］;xn=0.6.^n;n1=length(n);

m=［0:1:15］;n2=length(m);hn=ones(1,n2);

L=n1+n2-1;

Xk=fft(xn,L);Hk=fft(hn,L);Yk=Xk.*Hk;

yn=ifft(Yk,L)

figure(1);stem(abs(yn))

xlabel('n');ylabel('y(n)=x(n)*h(n)');

title('FFT实现快速卷积运算');

axis(［0,32,0,1.1*max(yn)］);



程序计算结果如下： 


yn =Columns 1 through 12

1.00001.60001.96002.17602.30562.38342.43002.4580
2.47482.48492.49092.4946

Columns 13 through 24

2.49672.49802.49882.49931.49930.89930.53930.3233
0.19370.11590.06930.0413

Columns 25 through 31

0.02450.01440.00840.00470.00260.00130.0005



FFT快速卷积结果的波形如图532 所示。



图532FFT计算线性卷积的结果


2） x(n)与h(n)点数相差很大

若序列x(n)的点数很多，即


NM


则有


L=N+M-1≈N


可得


k=NM2L1+32log2L≈NM2N1+32log2N=M2+3log2N（570）



当L与M相差很大时，k值将快速下降，线性卷积快速算法的优越性就难以体现，这时，可以采用分段卷积方法进行卷积快速运算。

应用中，有时会遇到一个短序列与一个长序列的卷积计算问题，显然，这时不宜直接应用上述卷积快速算法，否则浪费资源。可以考虑将长序列分解为与短序列点数相当或相等的若干段，分别计算每一段的卷积结果，然后用合适的方法将各段的计算结果合并，得到最终计算结果。由于分段后点数相当，故对每一段的卷积运算均可采用快速算法。对x(n)分段进行卷积运算的方法包括重叠相加法和重叠保留法两种。

（1） 重叠相加法。

重叠相加法（overlapadd）是指在计算分段卷积时需要对各输出段的重叠部分进行相加，因而称为重叠相加法。

设x(n)为长序列，h(n)为短序列，点数为M。将长序列x(n)分解为互不重叠的若干段，每段为N点，N与M相等或相差不大，令xi(n)表示x(n)的第i段： 


xi(n)=x(n),iN≤n≤(i+1)N-1
0,其他,i=0,1,…
（571）


则输入序列x(n)可表示为


x(n)=∑ixi(n)（572）
yi(n)=∑ixi(n)*h(n)



因此，x(n)与h(n)的线性卷积等于各xi(n)与h(n)的线性卷积之和，即


y(n)=x(n)*h（n）=∑ixi(n)*h(n)（573）


求和的每一卷积项都可采用线性卷积的快速算法进行计算。由于每一段的卷积结果yi(n)=∑ixi(n)*h(n)的长度为L≥N+M-1，L向上取值并符合FFT算法要求。

应用快速算法计算每一段的卷积前，应先对xi(n)及h(n)末端补零点，补到L点，以满足圆周卷积代替线性卷积yi(n)的条件，方便进行快速卷积运算： 


yi(n)=xi(n)h(n)


因为xi(n)为L点，而yi(n)为N=L+M-1点，故相邻两段卷积计算的结果将有（M-1）点产生重叠，即上一段卷积结果的后（M-1）点和下一段卷积运算结果的前（M-1）个点发生重叠。根据式（573），应将重叠部分相加，再与无重叠部分合并为y(n)的输出结果。

因此，采用FFT实现重叠相加法的步骤如下： 

① 用FFT计算H(k)=DFT［h(n)］;

② 用FFT计算Xi(k)=DFT［xi(n)］;

③ 计算乘积运算Yi(k)=Xi(k)H(k);

④ 用IFFT计算yi(n)=IDFT［Yi(k)］； 

⑤ 将各段卷积结果yi(n)相加，注意重叠区间的长度，y(n)=∑iyi(n)。

【例53】重叠相加法计算线性卷积示例。

已知序列x(n)和h(n)如下： 


x(n)=3n-1,0≤n≤15,h(n)=［1,0,-1］


试采用重叠相加法计算y(n)=x(n)h(n)。

解： 
对于两序列长度相差很大的线性卷积计算，可以采用重叠相加法。

MATLAB提供了基于FFT的重叠相加法计算线性卷积的函数fftfilt.m，该函数计算圆周卷积时采用快速傅里叶变换FFT。函数调用格式如下： 


y=fftfilt(h,x,L)


x和h分别表示序列x(n)和h(n)，L表示分段长度，该参数可用默认值，调用格式如下： 


y=fftfilt(h,x)


分段长度L为默认值时，系统自动对长序列选择一个合适的分段长度。

计算线性卷积的MATLAB程序如下： 


clc;clear all;close all

n=0:15;

x=3*n-1;h=［1,0,-1］;

y=fftfilt(h,x)


本题x(n)和h(n)的序列长度相差较大，为了加深对重叠相加法的理解，程序设置了分段长度为3和4两种情况，计算结果相同。

程序计算的卷积结果如下： 


y = Columns 1 through 11

-1.00002.00006.00006.00006.00006.00006.00006.0000
6.00006.00006.0000

Columns 12 through 16

6.00006.00006.00006.00006.0000



（2） 重叠保留法。

重叠保留法（overlapsave）与重叠相加法的相同之处是都需要先对序列x(n)进行分段，分为长度为N点的若干短序列。不同之处是，对分段之后的短序列的处理方式不同，重叠保留法的计算步骤如下。

① 分段。重叠保留法对短序列不进行补零操作，而是在xi（n）的开始部分补上前一段xi-1（n）保留下来的（M-1）个输入序列的值，组成长度为L=N+M-1点的第i段序列xi（n）。若L=N+M-1不满足2的整数幂，则可在每段序列末端补零值点以达到要求。

② 首段单独处理。分段的第一段x0（n）由于没有上一段，应特殊处理，这时可以直接将第一段x0（n）的前（M-1）个值赋为零。

③ 计算各段线性卷积。采用FFT计算h(n)和xi（n）的线性卷积，这时每段圆周卷积计算结果的前（M-1）点的值存在误差，不等于线性卷积的值，应舍去，保留后面的N-M+1点的值。

④ 修正后连接成卷积序列y(n)。完成每一段卷积的计算并修正到正确值，最后将修正后正确的卷积结果yi(n)连接成一个输出序列y(n)。

由于该方法在计算时，每一组（每一段）均由前一段保留下来的（M-1）个点和新的（N-M+1）个点所组成，故称为重叠保留法。

【例54】重叠保留法计算线性卷积示例。

已知序列x(n)和h(n)如下： 


x(n)=6n-1,0≤n≤15,h(n)=［1,0,-1］


试采用重叠保留法计算y(n)=x(n)*h(n)。

解： 
重叠保留法计算线性卷积的MATLAB程序如下： 


clc;close all;clear all;

n=0:12;x=6*n-1;h=［1,0,-1］;

L_x=length(x);L_h=length(h);

N=6;N=2^(ceil(log10(N)/log10(2)));

M=L_h-1;L=N-M;

h_fft=fft(h,N);

k=ceil(L_x/L);

x=［zeros(1,M),x,zeros(1,N-1)］;

y=zeros(k,N);

for i=0:k-1

xk=fft(x(i*L+1:i*L+N));

y(i+1,:)=real(ifft(xk.*h_fft));

end

y=y(:,L_h:N)'; %舍弃前面M个值

yn=(y(:))'


计算的卷积结果如下： 


yn =-1.00005.0000 12.0000 12.000012.0000 12.0000 12.0000 12.0000 

12.000012.000012.000012.0000 12.0000-65.0000 -71.0000



5.8.2线性相关的FFT算法

线性相关包括自相关与互相关，它们在信息与通信、统计和信号处理等众多领域得到了广泛应用。互相关函数和自相关函数都是信号分析与处理的重要概念，互相关函数表示两个不同时间序列之间的相关程度，它描述了信号x(n)与y(n)在任意两个不同时刻取值之间的相关程度，两个不同信号既可以是确定性信号，也可以是随机信号。一般情况下，信号相关的概念通常指两个信号的互相关。

应用FFT计算相关，即利用圆周相关来计算线性相关，常称之为快速相关计算。与线性卷积的快速算法类似，也应采用末端补零点的方法避免混叠失真。

设序列x(n)长度为N点，y(n)长度为M点，线性相关定义如下:


rxy(n)=∑M-1m=0x(n+m)y*(m)（574）


用圆周相关代替线性相关应满足L≥N+M-1，且L为2的整数幂，因此，需对序列末端按如下方式补零： 


x(n)=x(n),0≤n≤N-1
0,N≤n≤L-1
y(n)=y(n),0≤n≤M-1
0,M≤n≤L-1


快速相关计算的步骤如下： 

（1） 求L点FFT，X(k)=DFT［x(n)］； 

（2） 求L点FFT，Y(k)=DFT［y(n)］； 

（3） 求乘积，Rxy(k)=X(k)Y*(k)； 

（4） 求L点IFFT，rxy(n)=IFFT［Rxy(k)］。

由计算步骤可知，采用FFT算法计算线性相关的运算量与采用FFT计算线性卷积的运算量是相等的，而且可以应用已有的FFT算法程序进行计算。

【例55】FFT计算线性相关示例。

已知序列x(n)和y(n)如下： 


x(n)=［3,5,7,0,-1,3,5］,y(n)=R15(n)


试计算x(n)与y(n)的线性相关。

解： 
用FFT计算线性相关，应先对序列x(n)和y(n)的末端补零值点，直至长度满足条件： 


L≥N1+N2-1
N1=length［x(n)］，N2=length［y(n)］


计算线性卷积的MATLAB程序如下： 


clc;clear all;close all

xn=［3,5,7,0,-1,3,5］;n1=length(xn);

noise=randn(1,n1);

yn=［5,7,0,-1,3,5,3］;

n1=length(xn);n2=length(yn);

L=n1+n2-1;

Xk=fft(xn,L);Yk=fft(yn,L);

Rk=conj(Xk).*Yk;Rxy=ifft(Rk,L)

figure(1);n=length(Rxy);stem(abs(Rxy))

xlabel('n');ylabel('序列相关');

title('FFT实现快速相关计算');

axis(［0,16,0,1.1*max(Rxy)］);



线性相关的计算结果如下： 


Rxy =77.0000 18.0000 13.0000 47.0000 55.0000 30.00009.0000 25.0000 

50.000016.0000 -12.000047.0000109.0000


线性相关的波形如图533所示。



图533FFT计算线性相关的结果



习题
◆

1. 已知X(k)和Y(k)是两个N点实序列x(n)和y(n)的DFT，即


X（k）=DFT［x(n)］=∑N-1n=0x(n)z-n
Y（k）=DFT［y(n)］=∑N-1n=0y(n)z-n


为了提高计算效率，试设计用一个N点IDFT/IFFT计算出x(n)和y(n)的算法。

2. 试根据DIT算法原理，绘出当N=8时，DIT基2快速算法完整的蝶形运算图，其中输入为自然顺序，输出为倒位序，并标出正确的支路系数。

3. 若计算机运算一次复数乘法的平均时间为2μs，计算一次复数加法的平均时间为0.2μs，用该计算机求1024点DFT，试计算完成下列运算所需要的时间： 

（1） 使用DFT计算需要多少时间？

（2） 使用FFT计算需要多少时间？

4. 试根据DIT算法原理，绘出当N=16时，DIT基2快速算法完整的蝶形运算图，其中输入为倒位序，输出为自然顺序，并标出正确的支路系数。

5. 若有限长序列x(n)的长度N=64，试写出DIT基2 FFT算法（输入倒位序、输出自然顺序）完整蝶形图中第三级蝶形运算中的全部支路系数。

6. 试根据DIF算法原理，绘出当N=16时，DIF基2快速算法完整的蝶形运算图，其中输入为自然顺序，输出为倒位序，并标出正确的支路系数。

7. 已知x(n)为8点有限长序列，其定义如下： 


x(n)=1,0≤n≤7
0,其他n


若z平面螺线路径参数如下： 


A0=0.8，W0=1.2,θ0=π3,0=π10



试完成如下计算： 

（1） 采用CZT算法求其前10点的复频谱X（zk）； 

（2） 根据CZT实现原理，绘出zk的路径图。

8. CZT算法可用于计算N点有限长序列h(n)在z平面实轴上各zk点的z变换H(z)，使

（a） zk=ak，k=0，1，…，N-1，a为实数，a≠0

（b） zk=ak，k=0，1，…，N-1，a为实数，a≠±1

（c） (a)和(b)都可以

（d） (a)和(b)都不可以

上述说法中哪一说法是正确的？并简要分析原因。

9. 对信号x(t)=e－0.1t,t≥0进行频谱分析。

（1） 根据傅里叶变换求出其频谱； 

（2） 若以T=0.75s的抽样间隔对x(t)采样，求离散信号频谱的重复周期Ωs。

10. 已知序列h(n) 为N=8点的有限长序列，H(ejω)=DTFT［h(n)］。如果要计算H(ejω)在频率ωk=2π64k2(k=0，1，…，7)共8个频率抽样点处的值，要求采用绿色计算，即不能采用先算出超过8个抽样点，再舍弃一些点的计算方法，试分析采用CZT算法是否可行？并简要说明理由。

11. 已知FIR滤波器的单位抽样响应h(n)为M=60的有限长序列，若用该滤波器对一串很长的数据进行滤波，并采用基于重叠保留法的FFT实现这一滤波功能，离散傅里叶变换为128点。为了实现该滤波功能： （1）各输入段必须重叠多少个抽样点？（2）数据如何分段？（3）从每一段卷积结果中取出多少个点才能使这些值连接在一起得到正确的卷积序列（滤波器的输出）？（4）若对于长数据已分段为每段120点，应如何处理？

12. 已知序列x(n)和h(n)如下： 


x(n)={2,8,2,1,8,2,1,2,1,6},h(n)={-1,0,1}


试用基于重叠相加法的FFT算法计算序列x(n)和h(n)的线性卷积。

13. 已知序列x(n)和h(n)如下： 


x(n)={2,1,2,1,8,2,1,6,2}，h(n)=x(n-3)


试用FFT计算序列x(n)和h(n)的互相关序列。