第5章
CHAPTER 5


快速傅里叶变换






DFT和IDFT的求和可以简单地看成是在一个域中取出N个数并且在另一个域中产生N个(复)数的一种计算方法。虽然DFT是一种重要的数字信号处理工具，能实现信号时域和频域的转换，但是几乎很少使用。因为直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量太大，所以在快速傅里叶变换(Fast Fourier Transform，FFT)出现以前，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。
由于FFT的输出与DFT的输出相同，但运算量却少很多，所以被广泛用于计算机计算。本章将介绍FFT的经典算法并给出一些应用举例，内容包括： 
DFT的运算量及基本对策； 
按时域抽选(DITFFT)的基2FFT算法； 
FFT算法与DFT运算量的对比； 
按频域抽选的基2FFT算法； 
IDFT的高效算法； 
进一步减少运算量； 
FFT应用举例。






5.1DFT的运算量及基本对策
下面分析一下，为什么直接计算DFT，当N较大时，计算量会很大。
设x(n)为有限长序列，其N点DFT为
X(k)=∑N-1n=0x(n)WnkNk=0,1,…,N-1(5.1)
因为序列x(n)不一定是实数序列，所以考虑x(n)为复数序列的一般情况。若直接按式(5.1)计算，对任一个k值，比如求k=0时X(k)的值，
X(k=0)=x(0)W0×0N+x(1)W1×0N+x(2)W2×0N+…+x(N-1)W(N-1)×0N
则需要N次复数乘法和(N-1)次复数加法。因此，计算出X(k)全部的N个值，共需要N2次复数乘法和N(N-1)次复数加法运算。当N1时，N(N-1)≈N2。这样看来，对于一个N点DFT的复数乘法和复数加法的运算次数大概都为N2次。DFT运算量如表5.1所示。


表5.1DFT的运算量



复数乘法的次数复数加法的次数

求一个X(k)NN-1
求所有的X(k)N·NN·(N-1)

从上面的分析可以知道，当N较大时，DFT的运算量将非常大。例如N=1024时，N2=1048576，这对于实时信号处理滤波器的计算速度来说，将是难以实现的。因此，必须将DFT的运算量尽可能地减少，才能使DFT在各种科学和工程计算中得到真正的应用。
从表5.1可以看出，N点DFT的复数乘法次数与N2有关，如果能把N点DFT分解为几个较短点数的DFT，将使计算复数乘法的次数大大减少。
通过分析DFT变换公式中WmN的表达式，可以发现旋转因子WmN具以下这些性质： 
(1) 周期性
W(n+lN)·kN=Wn·(k+lN)N=WnkNl为整数(5.2)
(2) 共轭对称性
(WnkN)=W-nkN(5.3)
(3) 可约性
WnkN=WmnkmN=WnkmNmm为整数，Nm为整数(5.4)
由上面这些性质，还可得出一些特殊值
WN/2N=-1,Wk+N2N=-WkN,W(N-k)nN=W(N-n)kN=W-knN
利用这些性质，可以将DFT中的一些项合并。
因此，FFT算法减少运算量的基本途径就是不断地把长序列的DFT分解成几个短序列的DFT，并利用旋转因子WmN的周期性和对称性来减少DFT的运算次数。基2FFT算法是FFT算法中最简单最常用的，该算法分为两类： 时域抽选法FFT(DecimationInTime Fast Fourier Transform, DITFFT)； 频域抽选法FFT(DecimationInFrequency Fast Fourier Transform, DIFFFT)。
【随堂练习】
试证明旋转因子WnkN的下列性质： ①W(n+lN)·kN=Wn·(k+lN)N=WnkNl为整数； ②(WnkN)=W-nkN； ③WnkN=WmnkmN=WnkmNmm为整数，Nm为整数。







5.2按时间抽选的基2FFT算法
按时间抽选(DecimationInTime，DIT)的基2FFT算法可以简写为DITFFT算法，也称为库利图基算法。设序列x(n)长度为N，且满足N=2M，M为自然数。按n的奇偶性把x(n)分成以下两组长度为N2点的子序列： 
x1(r)=x(2r)r=0,1,…,N2-1(5.5)

x2(r)=x(2r+1)r=0,1,…,N2-1(5.6)
则x(n)的N点DFT为
X(k)=∑n=偶数x(n)WnkN+∑n=奇数x(n)WnkN

=∑N2-1r=0x(2r)W2krN+∑N2-1r=0x(2r+1)Wk(2r+1)N

=∑N2-1r=0x1(r)W2krN+WkN∑N2-1r=0x2(r)W2krN(5.7)
利用旋转因子WnkN的可约性，即 W2krN=e-j2πN2kr=e-j2πN/2kr=WkrN/2，上式可以表示为
X(k)=∑N2-1r=0x1(r)WkrN/2+WkN∑N2-1r=0x2(r)WkrN/2 

=X1(k)+WkNX2(k)k=0,1,…,N-1(5.8)
其中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，即
X1(k)=∑N2-1r=0x1(r)WkrN/2=DFT［x1(r)］N/2(5.9)

X2(k)=∑N2-1r=0x2(r)WkrN/2=DFT［x2(r)］N/2(5.10)
并且X1(k)和X2(k)均以N/2为周期。因为X(k)有N个点，用式(5.8)计算得到的只是X(k)的前面一半项数的值，若要用X1(k)和X2(k)表达出所有X(k)的值，还需利用旋转因子WnkN的性质Wk+N2N=-WkN，把X(k)表示为前后两部分： 
前半部分X(k)k=0,1,…,N2-1可表示为
X(k)=X1(k)+WkNX2(k)k=0,1,…,N2-1(5.11)
后半部分X(k)k=N2,…,N-1可表示为
Xk+N2=X1k+N2+Wk+N2NX2k+N2
=X1(k)-WkNX2(k)k=0,1,…,N2-1(5.12)



图5.1蝶形运算符号

这样，只需求出两个长度为N/2序列x1(r)和x2(r)的N/2点DFT X1(k)和X2(k)的值，就可以求出X(k)在0≤k≤N-1内的所有值，运算量大大减少。
式(5.11)和式(5.12)的运算可用图5.1所示的流图符号表示，称为蝶形运算符号。采用这种图示法，经过一次奇偶抽取分解后，一个8点DFT运算可用图5.2表示。图中，N=23=8,X(0)~X(3)由式(5.11)给出，而X(4)~X(7)则由式(5.12)给出。



图5.28点DFT一次时域抽取分解运算流图


由图5.1可见，要完成一个蝶形运算，需要一次复数乘法和两次复数加法运算。由图5.2容易看出，经过一次分解后，计算1个N点DFT共需要计算两个N/2点DFT和N/2个蝶形运算。而计算一个N/2点DFT需要(N/2)2次复数乘法和(N/2)(N/2-1)次复数加法。所以，按图5.2计算N点DFT时，总的复数乘法次数为
2N22+N2=N(N+1)2N1≈N22(5.13)
复数加法次数为
NN2-1+2N2=N22(5.14)
由此可见，仅仅经过一次分解，就使运算量减少近一半。既然这样分解对减少DFT的运算量是有效的，且N=2M，N/2仍然是偶数，故可以对N/2点DFT再作进一步分解。表5.2列出了经过一次分解与不进行分解DFT运算量的比较。


表5.2经过一次分解与不进行分解DFT运算量的比较




不分解X(k)
(长度为N)
一次分解
X1(k)长度为N2X2(k)长度为N2
DFT里的运算
复数乘法N2次N22次N22次

复数加法N(N-1)次N2N2-1次N2N2-1次
蝶形里的运算
复数乘法—N2次

复数加法—N2+N2=N次
合计
复数乘法N2次2N22+N2=N(N+1)2N1≈N22次

复数加法N(N-1)次2×N2N2-1+2×N2=N22次

与第一次分解相同，将x1(r)按奇偶分解成两个N/4点的子序列x3(l)和x4(l)，即
x3(l)=x1(2l)l=0,1,…,N4-1(5.15)

x4(l)=x1(2l+1)l=0,1,…,N4-1(5.16)
X1(k)又可表示为
X1(k)=∑N4-1l=0x1(2l)W2klN/2+∑N4-1l=0x1(2l+1)Wk(2l+1)N/2

=∑N4-1l=0x3(l)WklN/4+WkN/2∑N4-1l=0x4(l)WklN/4

=X3(k)+WkN/2X4(k)k=0,1,…,N2-1(5.17)
其中
X3(k)=∑N4-1l=0x3(l)WklN/4=DFT［x3(l)］N/4(5.18)

X4(k)=∑N4-1l=0x4(l)WklN/4=DFT［x4(l)］N/4(5.19)
同理，由X3(k)和X4(k)的周期性和WmN/2的对称性Wm+N4N/2=-WmN/2最后得到： 
X1(k)=X3(k)+WkN/2X4(k)k=0,1,…,N4-1(5.20)

X1k+N4=X3(k)-WkN/2X4(k)k=0,1,…,N4-1(5.21)

用同样的方法可计算出
X2(k)=X5(k)+WkN/2X6(k)k=0,1,…,N4-1(5.22)

X2k+N4=X5(k)-WkN/2X6(k)k=0,1,…,N4-1(5.23)
其中
X5(k)=∑N4-1l=0x5(l)WklN/4=DFT［x5(l)］N/4(5.24)

X6(k)=∑N4-1l=0x6(l)WklN/4=DFT［x6(l)］N/4(5.25)

x5(l)=x2(2l)l=0,1,…,N4-1(5.26)

x6(l)=x2(2l+1)l=0,1,…,N4-1(5.27)
这样，经过第二次分解，又将一个N/2点DFT分解为2个N/4点DFT以及式(5.20)和式(5.21)所示的N/4个蝶形运算，如图5.3所示。以此类推，经过M次分解，最后将N点DFT分解N个1点DFT和M级蝶形运算(每级有N/2个蝶形运算)，而1点DFT就是时域序列本身。一个完整的8点DITFFT运算流图如图5.4所示。图中用到关系式WkN/m=WmkN。图中输入序列虽不是顺序排列，但其排列也是有规律的。图中的数组A用于存放输入序列和每级运算结果。


图5.3将一个N点DFT按时域抽取分解为4个N/4的分解运算流图(N=8)




图5.48点DITFFT运算流图


在利用时域抽取法的过程中,若序列x(n)长度不满足N=2M,则可以在序列后面补零,使序列长度满足。补零后，时域点数增加，但有效数值不变，不会增加信号的信息，也不会提高DFT频谱的准确性，故频谱X(ejω)不变，只是频谱的抽样点数增加，抽样点的位置有所改变，频率间隔减小而已。
【例5.1】假定现有8点按时间抽取的FFT芯片，如何利用这些芯片计算一个24点DFT？
解： 一个24点DFT的定义是
X(k)=∑23n=0x(n)Wnk24
对x(n)以因子3进行分解，就可以将这个DFT分解为如下的3个8点DFT： 
X(k)=∑7n=0x(3n)W3nk24+∑7n=0x(3n+1)W(3n+1)k24+∑7n=0x(3n+2)W(3n+2)k24
=∑7n=0x(3n)Wnk8+Wk24∑7n=0x(3n+1)Wnk8+W2k24∑7n=0x(3n+2)Wnk8
所以，如果组成3个序列
x1(n)=x(3n)n=0,1,…,7

x2(n)=x(3n+1)n=0,1,…,7

x3(n)=x(3n+2)n=0,1,…,7
并利用8点FFT芯片计算出它们所对应的X1(k)、X2(k)和X3(k)，那么x(n)的24点DFT就可以通过将8点FFT的输出组合起来求得，组合的方法如下： 
X(k)=X1(k)+Wk24X2(k)+W2k24X3(k)
【随堂练习】
假定现有5点和3点按时间抽取FFT芯片，如何利用这些芯片计算15点DFT？
5.3DITFFT算法与直接计算DFT运算量的对比
根据表5.1所列，对于每个点，标准的DFT要计算N次复数乘法和(N-1)次复数加法。那么，如果要计算所有N个点，总共要N2次复数乘法和N(N-1)次复数加法。运算的次数常用来评价该运算的难易程度。因此，DFT的运算难度系数与N2成正比。根据上节DITFFT算法的分解过程及图5.4可知，当N=2M时，其运算流图最多可分解为M级蝶形，每一级都有N/2个蝶形运算构成。因此，FFT的每一级运算包括需要N/2次复数乘法和N次复数加法(每个蝶形运算需要两次复数加法)。所以，M级运算总的复数乘法的次数为
CM=N2·M=N2log2N(5.28)
复数加法的次数为
CA=N·M=Nlog2N(5.29)
可以看出，DITFFT的运算难度系数与Nlog2N成正比。
当N1时，N2N2log2N，所以DITFFT算法比直接计算DFT的复数乘法的运算次数大大减少。例如，N=210=1024时，
N2N21og2N=10485765120=204.8
这样，就使复数乘法的运算效率提高200多倍。同理，DITFFT的复数加法的运算效率也可以提高100多倍。
由于DITFFT运算如此高效，它的计算几乎总是在采样点数为2的整数次幂的基础上进行。当信号的采样点数不是2的整数次幂时，需要在信号的末尾补零。根据式(3.2)，多加的零值不会影响信号的DTFT。又因为DFT或FFT只是DTFT的在频域的采样形式，补零对信号的DFT或FFT特性没有影响。

图5.5为DITFFT算法和直接计算DFT所需复数乘法次数CM与变换点数N的关系曲线。由此图更加直观地看出FFT算法的优越性，显然，N越大时，优越性就越明显。



图5.5DITFFT算法与直接计算DFT运算效率比较曲线

5.4按频率抽选的基2FFT算法
由于DFT和IDFT的变换前后时域和频域序列的长度相同，求解表达式形式相近，因此，对照时域抽取法，又产生了基2FFT算法中的另一种按频域抽选(DecimationInFrequency，DIF)的基2FFT算法(DIFFFT)，又称为桑德图基算法。频域抽取法的实现过程如下： 






设序列x(n)长度为N=2M，首先将x(n)分成前后两半，其DFT表示为如下形式： 
X(k)=DFT［x(n)］=∑N-1n=0x(n)WknN

=∑N2-1n=0x(n)WknN+∑N-1n=N2x(n)WknN

=∑N2-1n=0x(n)WknN+∑N2-1n=0xn+N2Wkn+N2

=∑N2-1n=0x(n)+WkN2Nxn+N2WknN(5.30)
式中
WkN2N=(-1)k=1k为偶数

-1k为奇数

将X(k)分成了偶数组和奇数组两组。当k取偶数k=2m,m=0,1,…,N2-1时，
X(2m)=∑N2-1n=0x(n)+xn+N2W2mnN

=∑N2-1n=0x(n)+xn+N2WmnN/2(5.31)

当k取奇数k=2m+1,m=0,1,…,N2-1时，
X(2m+1)=∑N2-1n=0x(n)-xn+N2W(2m+1)nN

=∑N2-1n=0x(n)-xn+N2WnNWmnN/2(5.32)
令
x1(n)=x(n)+xn+N2
x2(n)=x(n)-xn+N2WnNn=0,1,…,N2-1
将x1(n)和x2(n)分别代入式(5.31)和式(5.32)中，可得
X(2m)=∑N2-1n=0x1(n)WmnN/2
X(2m+1)=∑N2-1n=0x2(n)WmnN/2(5.33)
式(5.33)表明，X(k)按奇偶k值分为两组，其偶数组是x1(n)的N/2点DFT，奇数组是x2(n)的N/2点DFT。x1(n)、x2(n)和x(n)之间的关系也可用图5.6所示的蝶形运算流图符号表示。图5.7表示N=8时第一次分解的运算流图。



图5.6DIFFFT 蝶形运算流图符号




图5.7DIFFFT 第一次分解运算流图(N=8)



由于N=2M，可将N/2点DFT再分成偶数组和奇数组，每个N/2点的DFT又可分成两个N/4点DFT，其输入序列分别是x1(n)和x2(n)各自按上下对半分开形成的4个子序列。图5.8表示了N=8时第二次分解运算流图，经过两次分解，便分解为4个两点DFT。图5.9表示了序列长度N=8时的完整DIFFFT运算流图。



图5.8DIFFFT 第二次分解运算流图(N=8)





图5.9DIFFFT第三次分解运算流图(N=8)








若序列x(n)长度为N=2M，经过M-1次分解，最后分解为N/2个两点DFT。在DIFFFT中，两点DFT就是一个基本蝶形运算。
这种算法是对X(k)进行奇偶抽取分解，因此被称为频域抽取法FFT。观察图5.4和图5.9可知，DIFFFT算法与DITFFT算法有下列相同点： 可以原位计算； 共有M级运算； 每级共有N/2个蝶形运算； 两种算法的运算次数也相同。
不同点是： 
(1) 输入和输出的顺序不同。DIFFFT算法输入序列为自然顺序，而输出为倒序排列； DITFFT算法输入序列为倒序排列，而输出为自然排列。
(2) 蝶形运算不同。DITFFT蝶形先相乘求DFT再相加(减)； 而DIFFFT蝶形先加(减)再相乘求DFT。
5.5IDFT的高效算法
DITFFT算法和DIFFFT算法也可以用于计算IDFT。比较DFT和IDFT的运算公式： 
X(k)=DFT［x(n)］=∑N-1n=0WknN

x(n)=IDFT［X(k)］=1N∑N-1k=0X(k)W-knN
可以得到求解IFFT的第一种算法： 将DFT表达式中的变换核WknN换成W-knN，最后再乘以1/N，就是IDFT运算公式。此外，DFT的流图输入的是x(n)，输出是X(k)； 而IDFT的流图输入的是X(k)，输出的是x(n)。因此，原来的DITFFT改为IFFT后，称为DIFIFFT更合适； DIFFFT改为IFFT后，称为DITIFFT更合适。以DIFFFT流图为例，可得到图5.10所示的IFFT流图。



图5.10N=8基2IFFT流图


第二种算法可以直接利用FFT程序。由于IDFT计算公式为
x(n)=IDFT［X(k)］=1N∑N-1k=0X(k)W-knN
两边取共轭，可得
x*(n)=1N∑N-1k=0X(k)W-knN*
=1N∑N-1k=0X*(k)WknN
=1NDFT［X*(k)］

x(n)=1N∑N-1k=0X*(k)WknN*=1N{DFT［X*(k)］}*(5.34)
所以，求解IFFT可以先将X(k)取复共轭，然后直接调用FFT子程序，进行DFT运算，最后将所得结果取复共轭并乘以1/N，便得到序列x(n)。这种方法虽然用了两次取复共轭运算，但可以与FFT共用同一程序，因此用起来也比较方便。
第三种算法也是利用现有的FFT程序。令
p(n)=∑N-1k=0X(k)WnkN
则有

x(n)=1Np(N-n)=1N∑N-1k=0X(k)W(N-n)kN=1N∑N-1k=0X(k)W-nkN(5.35)
所以这种共用FFT程序的步骤是： 
(1) 利用FFT程序由X(k)求出p(n)； 
(2) 计算1Np(N-n)即为x(n)［注意n=0时p(N)=p(0),x(0)=p(0)/N］。
5.6进一步减少运算量
通过观察以上DITFFT和DIFFFT算法，工程师们发现还可以进一步减少运算量，下面简单介绍几种方法。
5.6.1多类蝶形单元运算
由DITFFT运算流图可知，N=2M点FFT共需要MN/2次复数乘法。由旋转因子
WpN=WJN·2L-M=WJ·2M-LNJ=0,1,2,…,2L-1-1(5.36)
其中，
p=J·2M-L(5.37)
当L=1时，只有一种旋转因子W0N=1，所以第一级不需要乘法运算。当L=2时，共有两个旋转因子W0N=1和WN/4N=-j，因此，第二级也不需要乘法运算。在DFT中，又称其值为±1和±j的旋转因子为无关紧要的旋转因子，如W0N=1，WN/2N，WN/4N等。
因此，除去第一、二两级后，所需复数乘法次数应是
CM=N2(M-2)(5.38)

再来找一下各级中的无关紧要旋转因子。当L=2时，有两个无关紧要的旋转因子W0N和WN/4N，因为同一旋转因子对应着2M-L=N2L蝶形运算，所以第二级共有N22=N4个蝶形不需要复数乘法运算。当L≥3时，第L级的两个无关紧要的旋转因子减少复数乘法的次数为2N2L=N2L-1。因此，从L=3至L=M共减少的复数乘法次数为
∑ML=3N2L-1=2N∑ML=312L=N2-2(5.39)
这样，DITFFT的复数乘法次数为
CM=N2(M-2)-N2-2=N2(M-3)+2(5.40)
进一步观察FFT中存在一些特殊的复数运算以进一步减少复数乘法次数。一般实现一次复数乘法运算需要四次实数乘，两次实数加。但对WN/8N=(1-j)22这一特殊复数，任一复数(x+jy)与其相乘时，有以下等式： 
(1-j)22(x+jy)=22(x+jy-jx+y)
=22［(x+y)-j(x-y)］=R+jI
其中，R=2/2(x+y),I=-2/2(x-y)=2/2(y-x)，它只需要两次实数加和两次实数乘就可实现。这样，WN/8N对应的每个蝶形节省两次实数乘。在DITFFT运算流图中，从L=3至L=M(M>3)级，每级都包含旋转因子WN/8N，第L级中，WN/8N对应N2L个蝶形运算。因此从第三级至最后一级，旋转因子WN/8N节省的实数乘次数与式(5.40)相同。所以从实数乘运算考虑，计算N=2M点DITFFT所需实数乘法次数为
RM=4N2(M-3)+2-N2-2=N2M-132+10(5.41)
在基2FFT程序中，将蝶形单元运算分为N 类，包含了所有旋转因子的称为一类蝶形单元运算； 去掉WmN=±1的旋转因子的称为二类蝶形单元运算； 再去掉WmN=±j的旋转因子的称为三类蝶形单元运算； 如果再判断处理WmN=(1-j)22，则称为四类蝶形单元运算。后三种运算称为多类蝶形单元运算。显然，蝶形单元类型越多，编程的复杂程度越高，但当N较大时，可大大减少乘法运算量。例如，N=4096时，三类蝶形单元运算的乘法次数为一类蝶形单元运算量的75%。
对于其他旋转因子WmN=cos2πmN-jsin2πmN，由于这些是需要正弦和余弦函数值的计算，因此计算量很大。可以在FFT程序开始前，预先计算出WmN，m=0,1,…,N2-1，存放在数组中，作为旋转因子表，在程序执行过程中直接查表，这样可使运算速度大大提高，但占用的内存较多。
5.6.2减少实序列的FFT的方法
对于许多信号，如数字语音信号，数据x(n)是实数序列。如果把x(n)看成一个虚部为零的复序列进行计算，这就增加了存储量和运算时间。根据实序列的FFT的特点可以有多种FFT运算量的方法，其中之一是用N/2点FFT计算一个N点实序列的DFT。以下介绍这种方法。
设x(n)为N点实序列，取x(n)的偶数点和奇数点分别作为新构造序列y(n)的实部和虚部，即
x1(n)=x(2n)n=0,1,…,N2-1
x2(n)=x(2n+1)n=0,1,…,N2-1

y(n)=x1(n)+jx2(n)n=0,1,…,N2-1
对y(n)进行N/2点FFT，输出Y(k)，则
X1(k)=DFT［x1(n)］=Yep(k)k=0,1,…,N2-1

X2(k)=DFT［x2(n)］=-jYop(k)k=0,1,…,N2-1
根据DITFFT及式(4.16)和(4.17)，可得到X(k)的前N2+1个值： 
X(k)=X1(k)+WkNX2(k)k=0,1,…,N2(5.42)
式中，X1N2=X1(0),X2N2=X2(0)。由于x(n)为实序列，因此X(k)具有共轭对称性，X(k)的另外N/2点的值为
X(N-k)=X*(k)k=0,1,…,N2-1
下面来计算按以上方法运算速度的提高程度。计算N/2点FFT的复乘次数为N4(M-1)，计算式(5.41)的复乘次数为N2，因此用这种算法，计算X(k)所需复数乘法次数为N4(M-1)+N2=N4(M+1)。相对一般的N点FFT算法，上述算法的运算效率为η=N2MN4(M+1)=2MM+1，例如当N=2M=210时，η=2011，运算速度提高了82%。
自1965年基2FFT算法出现以来，经过人们不断研究探索，现在已提出了多种快速算法。例如分裂基FFT算法、离散哈特莱变换(DHT)、基4FFT、基8FFT、基rFFT、混合基FFT，以及进一步减少运算量的途径等内容，它们对研究新的快速算法都是很有用的。这里就不再赘述，相关内容请读者阅读相关的文献。
5.7FFT应用举例
在MATLAB信号处理工具箱中提供了函数fft和ifft以进行快速傅里叶变换和逆快速傅里叶变换。快速傅里叶变换函数fft的一种调用形式为： 
y=fft(x)(5.43)
式中，x是序列，y是序列的快速傅里叶变换，x可以为一向量或矩阵，若x为向量，则y为相同长度的向量； 若x为矩阵，则y是对矩阵的每一列向量进行FFT运算。
快速傅里叶变换函数fft的另一种调用形式为： 
y=fft(x,N)(5.44)
式中，N表示函数执行N点的FFT。若x为向量且长度小于N，则函数将x补零到长度N，若向量x的长度大于N，则函数截断x使之长度为N。
函数fft是用机器语言写成，不是用MATLAB命令写成，因此执行起来非常快。并且它是作为一种混合基算法写成的，如果运算的长度为2n，则就能使用一个高速的基2FFT算法。如果运算的长度不是2n，则fft执行一种称为混合基的FFT算法，计算速度慢。
应用ifft函数进行逆快速傅里叶变换，它与fft具用同样的特性，这里不再赘述。以下通过例题来说明用fft函数得到的频谱的特点。
【例5.2】已知信号x(t)=0.5sin(2πf1t)+2sin(2πf2t)，其中f1=20Hz,f2=50Hz，采用频率为200Hz,以N表示数据个数，FFT采样点数用L表示，试分别绘制出在N=128点和N=1024 点的幅频图。程序如下： 

fs=200

N=128

n=0:N-1

t=n/fs

x=0.5*sin(2*pi*20*t)+2*sin(2*pi*50*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,1,1)

plot(f,m1)

title('N=128')

grid on

fs=200

N=1024

n=0:N-1

t=n/fs

x=0.5*sin(2*pi*20*t)+2*sin(2*pi*50*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,1,2)

plot(f,m1)

title('N=1024')

grid on



图5.11对信号分别做128点和1024点的快速傅里叶变换幅频图

运行结果如图5.11所示,显然，整个频谱图以采样频率(200Hz)的一半(100Hz)为对称轴。因此，在利用fft函数对信号作谱分析时，取零频率到采样频率的一半即可。另一方面，N=128和N=1024均能很好地分辨两种频率成分： 20Hz和50Hz。根据式(4.70)，N=128比N=1024的频谱分辨率低，因此N=1024的幅频图中更尖端。最后，幅频图的振幅大小与采样点数直接相关，如果要得到真实振幅，则需将变换后的结果乘以2/N。为此，可将上述作图命令修改为： plot(f(1∶N/2),m1(1∶N/2)*2/N)，重新得到的幅频图如图5.12所示。


图5.12对幅度进行了修订且只取采样频率的一半的快速傅里叶变换幅频图


【例5.3】已知信号x(t)=0.5sin(2πf1t)+2sin(2πf2t)，其中f1=15Hz,f2=40Hz，采用频率为100Hz,以N表示数据个数，FFT采样点数用L表示，在下列情况下绘制其幅频图。
(1) N=32，L=32； 
(2) N=32，L=128； 
(3) N=136，L=128； 
(4) N=136，L=512。
试分析所用数据长度不同时对傅里叶变换结果的影响。程序如下： 

fs=100

l1=32

N=32

n=0:l1-1

t=n/fs

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,2,1)

plot(f(1:N/2),m1(1:N/2)*2/N)

title('N=32 L=32')

grid on

l2=128

N=32

n=0:l2-1

t=n/fs

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,2,2)

plot(f(1:N/2),m1(1:N/2)*2/N)

title('N=32 L=128')

grid on

l3=128

N=136

n=0:l3-1

t=n/fs

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,2,3)

plot(f(1:N/2),m1(1:N/2)*2/N)

title('N=136 L=128')

grid on

l4=512

N=136

n=0:l4-1

t=n/fs

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)

y=fft(x,N)

m1=abs(y)

f=(0:N-1)*fs/N

subplot(2,2,4)

plot(f(1:N/2),m1(1:N/2)*2/N)

title('N=136 L=512')

grid on

运行结果如图5.13所示。


图5.13不同采样个数和fft 计算数据个数对所得频谱的影响


现对以上结果进行分析： 
(1) N=32，L=32时，频率分辨率较低； 但由于数据个数与FFT采用的数据个数相等，因此，不需要添加零而导致其他的频率成分。将振幅乘以2/N后，得到了绝对大小的振幅； 
(2) N=32，L=128时，FFT运算需加零，致使振幅中出现了很多其他的成分，其振幅的幅度也由于加零而明显减少； 
(3) N=136，L=128时，FFT运算需要截断128个数据，这时频率分辨率较高，将振幅乘以2/N后，得到了绝对大小的振幅； 
(4) N=136，L=512时，FFT运算需加零，致使振幅中出现了很多其他的成分，其振幅的幅度也由于加零而明显减少； 但是由于含有信号的数据个数足够多，因此，其振幅谱分辨率仍较高。
【例5.4】运用fft函数对信号x=0.5sin(2π·3·n·dt)+cos(2π·10·n·dt) 进行512点的FFT运算，并滤去此信号中频率为8~15Hz的波。采样间隔dt=0.02。绘出滤波前后的振幅谱以及滤波后的时域信号。程序如下： 

dt=0.02

N=512

n=0:N-1

t=n*dt

f=n/(N*dt)

f1=3

f2=10

x=0.5*sin(2*pi*f1*dt)+cos(2*pi*f2*dt)

subplot(2,2,1)

plot(t,x)

y=fft(x,N)

subplot(2,2,2)

plot(f,abs(y)*2/N)

f3=8

f4=15

yy=zeros(1,length(y))

for m=0:N-1

if(m/(N*dt)>f3&m/(N*dt)<f4|m/(N*dt)>(1/dt-f4)&m/(N*dt)<(1/dt-f3))

yy(m+1)=0

else

yy(m+1)=y(m+1)

end

end

subplot(2,2,4)

plot(f,abs(yy)*2/N)

x1=ifft(yy)

subplot(2,2,3)

plot(t,real(x1))

运行结果如图5.14所示。


图5.14运用fft进行信号的滤波

习题
1. 设一个信号x(n)={0,1,2,3,4,5,6,7}。
(1) 用时域抽取法计算其FFT； 
(2) 用(1)中的结果，计算并画出信号的幅度频谱和相位频谱。
2. 对一个长度为1024点的信号： 
(1) 分别进行DFT和FFT计算，需要多少次复数乘和复数加？
(2) DFT的运算量是FFT的多少倍？
3. 试证明下面关于FFT的旋转因子的等式： 
(1) W02=W04=W08； 
(2) W14=W28。
4. 信号以8kHz进行采样，进行512点的FFT。
(1) 求FFT的频率间隔； 
(2) 将信号补零为4096个采样点，再次计算FFT，频率间隔又是多少？
5. 如果一台通用计算机的速度为平均每次复乘用时40ns，每次复加用时5ns，用它来计算512点的DFT［x(n)］，问直接计算需要多少时间？用FFT运算需要多少时间？若做128点快速卷积运算，问所需最少时间是多少？
6. 设x(n)=［1,2,1,2,1］,h(n)=［1,2,2,1］。
(1) 在时域求y(n)=x(n)h(n)； 
(2) 用FFT流图法来求y(n)，即求出H(k),Y(k)=H(k)X(k)，y(n)=IFFT［Y(k)］。