第5章 快速傅里叶变换（FFT）

离散傅里叶变换（DFT）在信号的频谱分析，系统的分析、设计与实现方面起着非常重要的作用，一个关键的原因就是计算DFT有许多高效快速的算法，快速傅里叶变换（Fast Fourier Transform，FFT）算法就是其中之一。需要强调的是，FFT并不是一种新的变换算法，只是DFT的一种快速实现算法。
由于直接计算DFT的运算量太大，很多年一直没有得到实际应用。直到1965年，库利（J.W.Cooley）和图基（J.W.Tukey）在论文An Algorithm for the Machine Calculation of Complex Fourier Series中提出快速算法FFT，使得DFT的运算大为简化，从而使DFT真正得到广泛应用。FFT算法的横空出世，破解了DFT运算量偏大的历史性难题，促进了数字信号处理突飞猛进的发展，因此往往把1965年视作数字信号处理元年。
本章首先对DFT运算复杂度进行分析，找出DFT运算的瓶颈所在和提速的可能性，随后介绍将DFT在时域逐渐分解为较小点数DFT运算的方法，即按时间抽取（Decimation in Time）的FFT算法，在DITFFT算法的基础上介绍按频率抽取（Decimation in Frequency）的FFT算法，最后介绍FFT算法在工程实现中的一些经验技巧。
5.1DFT运算复杂度分析
本节主要用乘法和加法次数来衡量DFT的运算复杂度，这是因为在计算机或者DSP芯片中，乘法和加法次数直接决定了算法执行效率和速度。
5.1.1DFT的运算瓶颈

设x(n)为N点有限长序列，其N点的DFT正变换为
X(k)=∑N-1n=0x(n)WnkN,k=0,1,…,N-1（5.1.1）

式中，WnkNΔe-j2πNnk表示旋转因子。
DFT反变换为
x(n)=1N∑N-1k=0X(k)W-nkN,n=0,1,…,N-1（5.1.2）
从DFT正变换和反变换公式可以看出，二者差别在于旋转因子指数符号相反，以及差一个常数因子1/N。因此DFT正变换和反变换运算量是完全相同的，只需讨论DFT正变换的运算量即可。
从DFT正变换公式可以看出，每计算一个X(k)，需要进行N次复数乘法和N-1次复数加法。一共有N个X(k)需要计算，就需要N2次复数乘法和N(N-1)次复数加法，这就是完成整个DFT正变换需要的计算量。
当N1时，直接计算DFT，复数乘法次数和复数加法次数都与N2成正比。随着N的增大，运算次数会急剧增加。例如，当N=8时，需要64次复数乘法，当N=1024时，就需要1048576次复数乘法，复数乘法次数达到了百万量级。如果信号要求实时处理，对运算速度的要求实在太高了。

5.1.2DFT的运算特点
为了提高DFT的运算效率，研究人员对其进行了大量研究，发现在式（5.1.1）中其实存在大量的冗余运算，并且这些冗余运算主要来自于旋转因子的反复相乘。
假定x(n)为4点的有限长序列，其4点的DFT结果可以写为
X(0)=x(0)W04+x(1)W04+x(2)W04+x(3)W04
X(1)=x(0)W04+x(1)W14+x(2)W24+x(3)W34
X(2)=x(0)W04+x(1)W24+x(2)W44+x(3)W64
X(3)=x(0)W04+x(1)W34+x(2)W64+x(3)W94（5.1.3）
如果严格按照式（5.1.1）的定义进行运算，则需要16次复数乘法和12次复数加法。
旋转因子WnkN具有以下4个特性： 
共轭对称性(WnkN)=W-nkN（5.1.4）
周期性WnkN=W(N+n)kN=Wn(N+k)N（5.1.5）
可约性WnkN=WmnkmN=Wnk/mN/m（5.1.6）
特殊值W0N=1,WN/2N=-1,WN/4N=-j,Wk+N/2N=-WkN（5.1.7）

如果充分利用旋转因子的这4个特点，则式（5.1.7）可化解为
X(0)=x(0)+x(2)+x(1)+x(3)
X(1)=x(0)-x(2)+x(1)-x(3)W14
X(2)=x(0)+x(2)-x(1)+x(3)
X(3)=x(0)-x(2)-x(1)-x(3)W14（5.1.8）
从式（5.1.8）可以看出，计算4点的DFT结果，此时仅需要1次复数乘法和8次复数加法，相比按照定义计算的运算量明显降低。
根据式（5.1.1）的定义，DFT的运算量与N2呈正比，显然N越小越有利，因此很自然就希望将大点数的DFT分解为若干个小点数的DFT来计算。在这种分解过程中，充分利用旋转因子的4个特性，为提高DFT运算速度提供了可能性，FFT算法就是基于这种基本思路发展起来的。
为了能够不断地进行分解，FFT算法要求DFT的运算点数N=2L，L为正整数。这种N为2的整数次幂的FFT算法，称作基2 FFT算法。按照运算特点，基2 FFT算法可分为按时间抽取（DIT）和按频率抽取（DIF）两种方案。
5.2按时间抽取（DIT）的FFT算法
5.2.1算法原理

设x(n)为N点长序列，N=2L，L为正整数。按照n的奇偶取值将x(n)分为两组： 
x0(r)=x(2r)

x1(r)=x(2r+1)（5.2.1）
其中r=0,1,…,N/2-1。
则序列x(n)的N点DFT结果为
X(k)=∑N-1n=0x(n)WnkN

=∑N2-1r=0x(2r)W2rkN+∑N2-1r=0x(2r+1)W(2r+1)kN

=∑N2-1r=0x0(r)W2rkN+WkN∑N2-1r=0x1(r)W2rkN（5.2.2）
根据旋转因子的可约性，即W2rkN=WrkN/2，可得
X(k)=∑N2-1r=0x0(r)WrkN/2+WkN∑N2-1r=0x1(r)WrkN/2

=X0(k)+WkNX1(k)（5.2.3）
其中，k的取值范围为0≤k≤N/2-1，X0(k)和X1(k)分别是x0(r)和x1(r)的N/2点DFT结果。
从式（5.2.3）的结果可以看出，一个N点的DFT可以由两个N/2点的DFT结果组合得到，但这种方法只得到X(k)的前一半结果。要想得到X(k)的后一半结果，还需用到旋转因子的周期特性，即WrkN/2=Wr(k+N/2)N/2，这样可以得到
X0k+N2=∑N2-1r=0x0(r)Wr(k+N/2)N/2=∑N2-1r=0x0(r)WrkN/2=X0(k)（5.2.4）
同理可得，
X1k+N2=X1(k)（5.2.5）
式（5.2.4）和式（5.2.5）说明，后半部分k值（N/2≤k+N/2≤N-1）对应的X0(k)和X1(k)，分别与前半部分k值（0≤k≤N/2-1）对应的X0(k)和X1(k)相等。
结合式（5.2.3）、式（5.2.4）和式（5.2.5）的结论，再利用旋转因子的特殊值，即Wk+N/2N=WN/2NWkN=-WkN，就可得出后半部分的X(k)，即
Xk+N2=X0(k)-WkNX1(k)（5.2.6）
从式（5.2.3）和式（5.2.6）可以看出，只需要计算前半部分的X0(k)和X1(k)，就可以得到全部范围内的X(k)，这样就大大降低了运算量。
可用图5.2.1来表示式（5.2.3）和式（5.2.6）的运算流图，因为该流图形如蝴蝶，故常称作“蝶形运算”。一个蝶形运算单元需要1次复数乘法和2次复数加法。如果支路上没有标出系数，则该支路系数默认为1。


图5.2.1蝶形运算单元



图5.2.2给出了将N点DFT分解为两个N/2点DFT的蝶形运算流图，以N=8为例。


图5.2.2按时间抽取，将N点DFT分解为两个N/2点DFT（N=8）



由5.1.1节的分析可知，直接计算一个N/2点的DFT运算，需要N22次复数乘法和N2N2-1次复数加法，那么计算两个N/2点的DFT运算，则需要2×N22=N22次复数乘法和NN2-1次复数加法。此外，根据图5.2.2，把两个N/2点的DFT结果合成一个N点的DFT结果，还需要N/2次蝶形运算，即还需要N2次复数乘法和2×N2=N次复数加法。因此，按照图5.2.2的分解方式，总共需要N22+N2次复数乘法和NN2-1+N=N22次复数加法，运算量相比直接计算N点DFT减少将近一半。
因为N=2L，故N/2仍是偶数，很自然地就想到还可以把每个N/2点DFT分解为两个N/4点DFT来计算，分解和组合的方式与前面介绍的思路完全一致。
按照r的奇偶取值将x0(r)分为两组： 
x2(m)=x0(2m)

x3(m)=x0(2m+1)（5.2.7）

其中m=0,1,…,N/4-1。
同样可以得到，
X0(k)=∑N4-1m=0x0(2m)W2mkN/2+∑N4-1m=0x0(2m+1)W(2m+1)kN/2

=∑N4-1m=0x2(m)WmkN/4+WkN/2∑N4-1m=0x3(m)WmkN/4

=X2(k)+WkN/2X3(k)（5.2.8）
其中k的取值范围为0≤k≤N/4-1，X2(k)和X3(k)分别是x2(m)和x3(m)的N/4点DFT结果。
同理可得后半部分的X0(k)，即N/4≤k+N/4≤N/2-1，
X0k+N4=X2(k)-WkN/2X3(k)（5.2.9）
图5.2.3给出N=8时，将一个N/2点DFT分解为两个N/4点DFT的蝶形运算流图。


图5.2.3按时间抽取，将N/2点DFT分解为两个N/4点DFT（N=8）



对x1(r)也按照r的奇偶取值分为两组，
x4(m)=x1(2m)

x5(m)=x1(2m+1)（5.2.10）
同理得到
X1(k)=X4(k)+WkN/2X5(k)（5.2.11）
X1k+N4=X4(k)-WkN/2X5(k)（5.2.12）
其中k的取值范围仍然为0≤k≤N/4-1，X4(k)和X5(k)分别是x4(m)和x5(m)的N/4点DFT结果。
根据旋转因子的可约性，将系数统一为WkN/2=W2kN，则N=8点的DFT可分解为4个N/4点的DFT，如图5.2.4所示。


图5.2.4按时间抽取，将N点DFT分解为4个N/4点DFT（N=8）



同理可知，用4个N/4点的DFT来计算N点的DFT，相比分解成两个N/2点DFT的方式，运算量又减少将近一半。
如此不断分解下去，一直分解到基本运算单元为2点DFT为止。对于N=8，图5.2.4已经是“分解到底”的情况，此时将“N/4点DFT”模块用蝶形运算单元来代替（图5.2.1），可得完整的8点DITFFT流图，如图5.2.5所示。


图5.2.5N=8点的基2 DITFFT流图


由于每一步分解都是按照输入序列的奇偶顺序进行分组的，这就是该方法称作“按时间抽取”的原因。
5.2.2算法特点
从前面的分析可以看出，基2 DITFFT算法的推导过程虽然看起来有些复杂，但实现起来却非常有规律，主要有以下这些特点。
1. 同址运算
同址运算也称作原位运算，是指输入和输出在运算过程中占用相同的存储地址，从而可以节约存储空间，降低硬件成本，FFT运算就是典型的同址运算。
从图5.2.5可以看出，FFT运算由多级蝶形运算构成。若最初的输入数据保存在A0～AN-1的N个存储单元中，第1级运算完成之后的结果仍可以保存在A0～AN-1存储单元，原来保存在同样位置的输入数据被覆盖，相当于对数据进行了迭代更新。第1级的输出为第2级的输入，以此类推，最终的FFT结果也保存在A0～AN-1存储单元。也就是说，中间过程的所有结果都没有（不需要）被保存，存储数据空间始终为N个存储单元。
此外，还需要N/2个存储单元用于存放蝶形运算支路的系数，这些系数随着蝶形运算的迭代更新即可，也不需要额外的存储单元。
2. 输入倒序，输出顺序
从图5.2.5可以看出，输入数据是按照x(0)、x(4)、x(2)、x(6)、x(1)、x(5)、x(3)、x(7)送入DITFFT流图进行后续运算的。乍一看，这一顺序杂乱无章，其实是很有规律的，对应关系就是二进制的倒位序关系。比如十进制数n=6，用二进制表示为110，将二进制的110高低位顺序反转得到011，二进制的011转换为十进制为n=3，这意味着x(6)需要调整到x(3)的位置上去。通过这种变换关系，x(3)也需要调整到x(6)的位置上去，因此调整输入数据位置的过程是成对互换，并不需要额外的存储单元。
这种位置互换关系的根本原因在于“每一步分解都是按照输入序列的奇偶顺序进行分组的”，每次分组都是将二进制最低位进行了移动，分到最后刚好将所有位置全部进行了反转。按位反转在计算机中很容易实现，这也是基2 FFT算法受到广泛欢迎的重要原因之一。
表5.1以N=8为例，给出了输入数据调整位置前后的对应关系。对于一般情况N=2L，调整位置的过程是完全一致的。


表5.1自然顺序与倒位序的对应关系（N=8）


自然顺序
（十进制）自然顺序
（二进制）倒位序
（二进制）倒位序
（十进制）


00000000
10011004
20100102
30111106
41000011
51011015
61100113
71111117


3. 蝶形运算
从图5.2.5可以看出，对于N=2L点长序列，其DITFFT流图中共有L级蝶形运算，每级中都有N/2个蝶形运算单元。如果序列长度不是2的整数次幂，就补零达到这个要求。
在DITFFT流图中，蝶形运算类型（即支路系数和数据点间隔）随着迭代次数成倍增加。以N=8为例，在第1级蝶形运算中，只有一种支路系数，即W08，并且参与蝶形运算的两个数据点间隔为1。在第2级蝶形运算中，有两种支路系数，即W08和W28，并且参与蝶形运算的两个数据点间隔为2。在第3级蝶形运算中，有四种支路系数，即W08、W18、W28和W38，并且参与蝶形运算的两个数据点间隔为4。
每一级蝶形运算的上半部分没有旋转因子，只有下半部分有旋转因子，并且这些旋转因子的取值也是非常有规律的。旋转因子WrN完全由变量r决定，而r的个数和取值完全由蝶形运算的级数i决定，其中i=1,2,…,L，N=2L。r的个数为2i-1，取值为r=N2ip，p=0,1,…,(2i-1-1)，每级旋转因子的具体取值如表5.2所示。


表5.2DITFFT流图中旋转因子WrN的取值规律


第i级r个数r取值旋转因子作用


i=11个r=02点DFT

i=22个r=0,N44点→2点

i=34个r=0,N8,2N8,3N88点→4点

…………


i=L2L-1个r=0,N2L,…,(2L-1-2)N2L,(2L-1-1)N2LN点→N/2点

5.2.3运算量分析
对于N=2L点长序列，共有L级蝶形运算，每级中都有N/2个蝶形运算单元，故DITFFT流图中共有N2L个蝶形运算单元。每个蝶形运算单元需要1次复数乘法和2次复数加法，因此，实现全部N点的FFT需要N2L=N2log2N次复数乘法和2N2L=Nlog2N次复数加法。
从5.1.1节的分析可知，直接计算N点的DFT，需要N2次复数乘法和N(N-1)次复数加法。由于计算机中，乘法运算更耗时间，因此主要考虑复数乘法次数，表5.3将直接计算DFT和FFT的运算量进行了对比。可以看出，与直接计算DFT相比，FFT算法在运算量上有数量级的下降，并且N越大，FFT算法优势越明显。


表5.3直接计算DFT和FFT算法的复数乘法次数对比


N直接计算DFTFFT算法


241
4164
32102480
12816384448
256655361024
5122621442304
102410485765120
2048419430411264
40961677721624576

例5.1利用DITFFT流图，计算x(n)={1,3,5,2}0的DFT结果。
解： N=4点的DITFFT运算流图如图5.2.6所示。


图5.2.6N=4点DITFFT流图


代入x(n)={1,3,5,2}0，以及W0N=1，W1N=W14=-j，可得运算过程如图5.2.7所示。


图5.2.7例题5.1运算过程



从图5.2.7的输出结果可以看出，DFT结果为X(k)={11,-4-j,1,-4+j}0。
5.3按频率抽取（DIF）的FFT算法
对于N=2L的情况，另外一种常用的FFT算法是按照输出序列X(k)的奇偶顺序进行分解的，称作按频率抽取（Decimation in Frequency）的FFT算法。
先把输入序列x(n)按照前后对半分开，则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)WnkN+WN2kN∑N2-1n=0xn+N2WnkN（5.3.1）
其中k的取值范围为0≤k≤N-1，且WN2kN=WN2Nk=(-1)k。
式（5.3.1）右边并不是两个N/2点DFT之和的形式，因为旋转因子是WnkN，而不是WnkN/2，继续整理可得
X(k)=∑N2-1n=0x(n)+(-1)kxn+N2WnkN（5.3.2）
根据k的奇偶取值，可以把X(k)分成两部分，
X(2r)=∑N2-1n=0x(n)+xn+N2W2rnN

=∑N2-1n=0x(n)+xn+N2WnrN/2（5.3.3）
X(2r+1)=∑N2-1n=0x(n)-xn+N2W(2r+1)nN

=∑N2-1n=0x(n)-xn+N2WnNWnrN/2（5.3.4）
其中r的取值范围为0≤r≤N/2-1。
在此设
x0(n)=x(n)+xn+N2
x1(n)=x(n)-xn+N2WnN（5.3.5）
则式（5.3.3）和式（5.3.4）分别表示x0(n)和x1(n)的N/2点DFT结果，即
X(2r)=DFT［x0(n)］=∑N2-1n=0x0(n)WnrN/2（5.3.6）
X(2r+1)=DFT［x1(n)］=∑N2-1n=0x1(n)WnrN/2（5.3.7）
此时，N点DFT结果X(k)按照k的奇偶取值分成了两个N/2点的DFT来计算。图5.3.1给出了将N点DFT按频率抽取分解为两个N/2点DFT的蝶形运算流图。


图5.3.1按频率抽取，将N点DFT分解为两个N/2点DFT（N=8）


与按时间抽取的思路类似，还可以把N/2点的DFT分解为两个N/4点的DFT，分解依据仍然按照频率顺序的奇偶取值，一直分解到基本运算单元为2点DFT为止。N=8点的完整DIFFFT流图，如图5.3.2所示。


图5.3.2N=8点的基2 DIFFFT流图



从前面的分析可知，基2 DITFFT算法是根据输入序列的奇偶顺序不断进行分组，基2 DIFFFT算法是根据输出序列的奇偶顺序不断进行分组，这两种计算DFT的算法思路和运算流图非常相似，主要有以下特点。
（1） DIFFFT算法与DITFFT算法运算量相同。
如果把FFT算法看作是一个系统，那么只需要把DITFFT流图中输入和输出进行交换，再把所有箭头反向，所有数据就会反向“流动”，新的流动方式即DIFFFT流图。DIFFFT流图与DITFFT流图的这种倒置关系，意味着二者运算量完全相同，即都是N2log2N次复数乘法和Nlog2N次复数加法。
（2） 同址运算。
与DITFFT算法类似，DIFFFT运算也满足同址运算的特点，即每一级运算数据都可以保存在同样地址的N个存储单元中，且蝶形运算支路的系数也只需要N/2个存储单元即可。
（3） 输入顺序，输出倒序。
DIFFFT流图与DITFFT流图呈倒置关系，二者的输入和输出进行了交换，当然输入和输出位置关系的特点也进行了交换。
（4） 蝶形运算。
与DITFFT算法类似，DIFFFT运算也满足蝶形运算的特点。区别在于： DITFFT流图中蝶形运算类型随着迭代次数成倍增加，DIFFFT流图中蝶形运算类型随着迭代次数成倍减少。DIFFFT流图的这个特点从图5.3.2中也很容易看出来。
因为DIFFFT流图与DITFFT流图呈倒置关系，即将DITFFT流图中输入和输出进行交换，再把所有箭头反向就可得到DIFFFT流图，而旋转因子的位置和取值并不会发生变化，因此不用专门去考虑DIFFFT流图中的旋转因子。
为便于大家学习掌握DITFFT算法与DIFFFT算法，表5.4对它们的特点进行了总结。


表5.4DITFFT算法与DIFFFT算法特点


DITFFT算法DIFFFT算法



同址运算同址运算
输入倒序，输出顺序输入顺序，输出倒序
蝶形运算蝶形运算

例5.2请利用DIFFFT流图重做例5.1。
解： N=4点的DIFFFT运算流图如图5.3.3所示。


图5.3.3N=4点 DIFFFT流图



代入x(n)={1,3,5,2}0，以及W0N=1，W1N=W14=-j，可得运算过程如图5.3.4所示。


图5.3.4例5.2运算过程



从图5.3.4的输出结果可以看出，DFT结果为X(k)={11,-4-j,1,-4+j}0，注意此时需要将输出结果进行倒序排列。
5.4FFT算法的应用技巧
前面主要以数学公式和流图的形式介绍了基2 DITFFT算法和基2 DIFFFT算法，但在FFT算法的硬件和软件实现中还有一些实际问题需要考虑。此外，在实际应用中如果能够充分利用FFT算法的特点，还可以进一步提高处理效率。
1. 查表法计算旋转因子
根据定义，旋转因子WrN=e-j2πNr需要通过指数函数计算得到。如果在FFT运算中旋转因子是实时产生的，就会对FFT算法的运算速度产生较大影响。因为在计算机中指数运算耗时较多，并且在FFT流图中旋转因子会反复参与运算，使得指数运算反复进行。
一个有效的解决方案就是通过查表法来计算旋转因子，即事先将所有N/2个旋转因子计算出来，存储在特定位置成为一个旋转因子表。在FFT算法运行时，通过查表来得到所需的旋转因子取值，这样就可使运算速度大为提高。
2. 悄无声息的倒序操作
以DITFFT算法为例，“输入倒序，输出顺序”为该算法的特点之一，但如果需要事先将输入数据x(n)进行倒序操作，这是非常不方便的。
在实际的FFT函数或处理器中，不必考虑倒序操作，计算机会通过变址寻址的方式去寻找相应位置的输入数据。因此对于FFT函数或处理器而言，都是“输入顺序，输出顺序”的，使用者不必理会内部采用的是DITFFT流图还是DIFFFT流图，只需要按照自然顺序将时域序列x(n)输入即可，输出的也是按照自然顺序排列的频域样本序列X(k)。
3. 用FFT实现IFFT
DFT算法可以通过FFT算法快速实现，其反变换IDFT算法也有相应的IFFT算法来快速实现。并不需要专门去推导和研究IFFT算法，因为用FFT算法就可以实现IFFT算法。
由式（3.2.7）可知，IDFT的定义为
x(n)=1N∑N-1k=0X(k)ej2πNkn（5.4.1）
其中n的取值范围为0≤n≤N-1。
仔细观察DFT和IDFT公式，二者的最大区别在于旋转因子指数符号相反，互为共轭关系。因此，对式（5.4.1）两边取共轭，可得
x(n)=1N∑N-1k=0X(k)e-j2πNkn=1NDFT［X(k)］（5.4.2）
DFT运算可由FFT运算快速实现，对式（5.4.2）两边再次取共轭可得x(n)，即可用FFT算法实现IFFT算法，
x(n)=1N{FFT［X(k)］}（5.4.3）
4. 实序列的FFT
FFT流图中的旋转因子都为复数，因此无论输入的x(n)是复数还是实数，在第1级蝶形运算之后都是复数运算。现实世界中的许多信号都是实信号，如果按照前面介绍的方法进行FFT运算，相当于把实信号看作虚部为零的复信号，这样就会“浪费”x(n)的虚部信息，降低了运算效率。
提高运算效率的基本思路就是把虚部信息利用起来。对于两个实序列，可以把其中一个实序列作为实部，另外一个实序列作为虚部来进行FFT运算； 对于一个较长的实序列，可以参照DITFFT的思想，把偶数序号部分作为实部，奇数序号部分作为虚部来进行FFT运算。
第一种情况，在例3.6中就已经得到解决。对于实序列x1(n)和x2(n)，长度都是N，将x1(n)和x2(n)组合成一个N点的复序列y(n)=x1(n)+jx2(n)。
根据例3.6的结论可知，x1(n)和x2(n)的N点FFT结果为
X1(k)=12［Y(k)+Y(N-k)］
X2(k)=12j［Y(k)-Y(N-k)］（5.4.4）
其中k的取值范围为0≤k≤N-1，Y(k)表示复序列y(n)的FFT结果。为了能利用基2的FFT算法，要求序列长度N=2L，L为正整数。如果不满足，则补零满足这个条件。
下面分析用组合复数的方法所需要的运算量（仅讨论复数乘法次数）。X1(k)和X2(k)可以通过式（5.4.4）由蝶形运算实现，1次蝶形运算需要1次复数乘法，故N点的X1(k)和X2(k)需要N次复数乘法。此外，计算N点长复序列y(n)的FFT需要N2log2N次复数乘法，故组合复数的方法总共需要N2log2N+N次复数乘法。如果直接计算实序列x1(n)和x2(n)的N点FFT结果，则需要2×N2log2N=Nlog2N次复数乘法。表5.5给出了直接法和组合复数法计算两个N点长实序列的FFT所需复数乘法次数。


表5.5计算两个N点长实序列的FFT所需复数乘法次数


N直接法组合复数法



128896576
25620481280
51246082816
1024102406144
20482252813312
40964915228672

第二种情况，对于一个较长的实序列，同样可以利用组合复数法来降低运算量。
对于N=2L点的实序列x(n)，按照n的奇偶取值拆分成两个短序列，即
x0(r)=x(2r)

x1(r)=x(2r+1)（5.4.5）
其中r的取值范围为0≤r≤N/2-1。
将x0(r)和x1(r)组合成一个复序列y(r)=x0(r)+jx1(r)，该复序列的长度为N/2。根据前面结论可知，x0(r)和x1(r)的FFT结果为
X0(k)=12［Y(k)+Y(N-k)］
X1(k)=12j［Y(k)-Y(N-k)］（5.4.6）
参照DITFFT的思想，通过蝶形运算可将X0(k)和X1(k)组合得到X(k)，即
X(k)=X0(k)+WkNX1(k)
Xk+N2=X0(k)-WkNX1(k)（5.4.7）
其中k的取值范围为0≤k≤N/2-1。
下面分析把长实数序列拆开组合成短复数序列，再进行FFT运算需要的复数乘法次数。X(k)可以通过式（5.4.7）由蝶形运算实现，N/2点长的X(k)只需要N/2次蝶形运算，故计算X(k)需要N/2次复数乘法。根据式（5.4.6），X0(k)和X1(k)也可通过蝶形运算实现，N/2点长的X0(k)和X1(k)需要N/2次蝶形运算，即N/2次复数乘法。此外，计算N/2点长复数序列y(n)需要N/22log2N2=N4(log2N-1)次复数乘法，故总共需要N4(log2N-1)+N2+N2次复数乘法。直接计算N点长实序列x(n)的FFT，需要N2log2N次复数乘法。表5.6对这两种方法的运算量进行了对比。


表5.6计算一个N点长实序列的FFT所需复数乘法次数


N直接法拆分组合复数法



128448320
2561024704
51223041536
102451203328
2048112647168
40962457615360

习题
1. 如果通用计算机计算一次复数乘法需要40ns，计算一次复数加法需要5ns，用它来计算512点的DFT，请问直接计算需要多少时间？用FFT算法来计算需要多少时间？
2. 已知序列x(n)和y(n)都为100点长，用FFT来计算这两个序列的线性卷积，假设计算机性能同第1题，请问需要多少时间？
3. 已知1024点序列x(n)的DFT结果为X(k)，假设计算机性能同第1题，请问用FFT来计算IDFT结果需要多少时间？
4. 一个实时卷积器是用FFT的重叠保留法分段处理的，假定系统的单位脉冲响应长度为128，每段1024点FFT的运算时间为0.2s，一次复数乘法的时间为200μs，不计数据采集、存储和加法的运算时间，请问： 
（1） 该系统的采样频率最高可以达到多少？
（2） 若两路信号同时卷积，则采样频率最高是多少？
5. 设序列x(n)是一个M点的有限长序列，0≤n≤M-1，x(n)的z变换为X(z)。现将X(z)在单位圆上进行N点的等间隔采样，试分别针对N≤M和N>M这两种情况，找出只用一个N点FFT就能计算X(z)的N个采样值的方法。
6. 已知序列x(n)={1,5,3,4}0
（1） 试给出DITFFT的信号流图，注意标出节点系数； 
（2） 利用流图计算输出结果X(k)； 
（3） 计算X(k)的幅度谱和相位谱。
7. 已知X(k)={11,-4-j,1,-4+j}0（例5.1），请用DITFFT的信号流图来计算x(n)。
8. 用重叠保留法来完成以下滤波功能，其中h(n)长度为31，信号x(n)的长度为19000，请利用512点的FFT算法来实现滤波运算。
9. 同上题的数据，请采用重叠相加法来实现滤波运算。
10. 请编写重叠相加法的MATLAB函数，并用来验证例4.6的结果。
11. 请编写重叠保留法的MATLAB函数，并用来验证例4.7的结果。
12. 请编写利用FFT实现IFFT运算的MATLAB函数，并与MATLAB函数ifft进行比较，验证X(k)={11,-4-j,1,-4+j}0的IFFT结果。