1、第五章第五章 快速傅立叶变换快速傅立叶变换(FFT)5.1 引言引言5.2 直接计算直接计算DFT的特点及减少运算量的基本途径的特点及减少运算量的基本途径5.3 基基2FFT算法算法5.1 引言引言影响数字信号处理发展的最主要因素之一是处理速度。DFT使计算机在频域处理信号成为可能，但是当N很大时，直接计算N点DFT的计算量非常大。快速傅立叶变换(FFT：Fast Fourier Transform)可使实现DFT的运算量下降几个数量级，从而使数字信号处理的速度大大提高。自从1965年第一篇DFT快速算法的论文发表以来，人们已经研究出多种FFT算法，它们的复杂度和运算效率各不相同。本章主要介绍
2、最基本的基2FFT算法及其编程方法。5.2 直接计算直接计算DFT的特点及减少运算量的基本途径的特点及减少运算量的基本途径长度为N的序列x(n)的N点DFT为 的周期性：10)()()(NnknNnWnxnxDFTkX点k=0,1,N-1mNmNjlNmNjlNmNWeeW2)(2(5.2.1)mNW 的对称性：(5.2.2)(5.2.3)mNNmNWW2mNWmNmNNWW*222222NNjNjNWeeW5.3 基基2FFT算法算法5.3.1 DIT-FFT算法序列x(n)的N(N=2M)点DFT为knNNnNWnxnxDFTkX10)()()(点k=0,1,N-1将上面的和式按n的奇偶性
3、分解为令x1(l)=x(2l),x2(l)=x(2l+1)。因为W2klN=WklN/2,所以上式可写成)12(12/0212/0)12()2()()()(lkNNllkNNlknNnknNnWlxWlxWnxWnxkX奇数偶数klNNlkNklNNlWlxWWlxkX2/12022/12/01)()()(5.3.1)(5.3.1)式说明，按n的奇偶性将x(n)分解为两个N/2长的序列x1(l)和x2(l)，则N点DFT可分解为两个N/2点DFT来计算。用X1(k)和X2(k)分别表示将以上两式代入(5.3.1)式，并利用 和X1(k)、X2(k)的隐含周期性可得到：)3.3.5(,12,.,
4、1,0,)()()()2.3.5(,12,.,1,0,)()()(12/02/22/2212/02/12/11NkWlxlxDFTkXNkWlxlxDFTkXNlklNNNlklNN点点kNNkNWW2这样，就将N点DFT的计算分解为计算两个N/2点离散傅立叶变换X1(k)和X2(k)。为了将如上分解过程用运算流图表示，以便估计其运算量，观察运算规律，总结编程方法，先介绍一种表示上式的蝶形图。蝶形图及其运算功能如图5.3.1所示。12,.1,0)()()2()()()(2121NkkXWkXNkXkXWkXkXkNkN图5.3.1 蝶形运算图ABCA CBA CB 根据图5.3.2可以求得第一
5、次分解后的运算量。图5.3.2包括两个N/2点DFT和N/2个蝶形，每个 点DFT需要 次复数乘法和 次复数加法运算，每个蝶形只有一次复数乘法运算和两次复数加法运算。所以，总的复数乘法次数为总的复数加法次数为2N2)2(N2)12(NN2|)1(222)2(212NNNNNN22222)12(2NNNN图 5.3.2 8点DFT一次时域抽取分解运算流图N/2点DFTWN0N/2点DFTWN1WN2WN3x(0)X1(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)X(0)X(1)X(2)X(3)X(4)X(5)X
6、(6)X(7)N=8点DIT-FFT的运算流图如图5.3.3(a)所示。根据WkN/m=WkmN，将图5.3.3(a)转换成如图5.3.3(b)所示的标准形式的运算流图。（a）图5.3.3 8点DIT-FFT运算流图04NWx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)
7、A(7)04NW04NW04NW12NW02NW12NW02NW0NW1NW2NW3NWx(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)0NW1NW2NW3NW2NW0NW2NW0NW0NW0NW0NW0NW（b）图5.3.3 8点DIT-FFT运算流图5.3.2
8、 DIT-FFT的运算效率DIT-FFT的运算效率指直接计算DFT的运算量与DIT-FFT的运算量之比。由图5.3.3可见，N=2M时，其DIT-FFT运算流图由M级蝶形构成，每级有N/2个蝶形。因此，每级需要N/2次复数乘法运算和N次复数加法运算，M级形共需复数乘法次数CM(2)和复数加法次数CA(2)分别为 (5.3.5)(5.3.6)NNMNCM2log22)2(NNMNCA2log)2(直接计算N点DFT的复数乘法次数为N2,复数加法次数为(N-1)。当N1时,N2/CM(2)1，所以N越大，DIT-FFT运算效率越高。DIT-FFT算法与DFT所需乘法次数与N的关系曲线如图5.3.4
9、所示。例如，N=210=1024时，DIT-FFT的运算效率为而当N=211=2048时，8.20410210241024)2(22MCNFFTDITDFT的乘法次数的乘法次数37.372112048222)2(22MNMNNCNM图5.3.4 DIT-FFT与DFT所需乘法次数比较曲线N(取 样 点 数)1024512256128640直 接 计 算 DFTDIT-FFT算 法64128 2565121024乘 法 次 数(1000)5.3.3 DIT-FFT的运算规律及编程思想1.原位计算 由图5.3.3可以看出，DIT-FFT的运算过程很有规律。2.旋转因子的变化规律观察图5.3.3(a
10、)不难发现，第L级共有2L-1个不同的旋转因子。N=23=8时的各级旋转因子表示如下：L=1时，J=0L=2时，J=0,1L=3时，J=0,1,2,3JJNpNLWWW24JJNpNLWWW22JJNpNLWWW2对N=2M的一般情况，第L级的旋转因子为 J=0,1,2,2L-1-1由于 2L=2M2L-M=N2L-M所以 J=0,1,2,2L-1-1，(5.3.7)(5.3.8)这样，就可按(5.3.7)式和(5.3.8)式确定第L级运算的旋转因子(实际编程序时，L为最外层循环变量)。JpNLWW2LMMLJNJNpNWWW22LMJp23.蝶形运算规律设序列x(n)经时域抽选(倒序)后，存
11、入数组X中。如果蝶形运算的两个输入数据相距B个点，应用原位计算，则蝶形运算可表示成如下形式：XL(J)XL-1(J)+XL-1(J+B)WpN XL(J+B)XL-1(J)-XL-1(J+B)WpN式中 p=J2M-L;J=0,1,2L-1-1;L=1,2,M设式中，下标R表示取实部，I表示取虚部，)()()(11JXjJXJXRLIIIRRRIIIRRRIRLRIIIRRTJXBJXTJXBJXTJXJXTJXJXJjXJXJXpNBJXpNBJXTpNBJXpNBJXT)()()()()()()()()()()(2sin)(2cos)(2sin)(2cos)(IRpNLjTTWBJXT)(
12、14.编程思想及程序框图仔细观察图 5.3.3，还可归纳出一些编程序有用的运算规律：第L级中，每个蝶形的两个输入数据相距B=2L-1个点；同一旋转因子对应着间隔为2L点的2M-L个蝶形。总结上述运算规律，便可采用下述运算方法。先从输入端(第1级)开始，逐级进行，共进行M级运算。当进行第L级运算时，依次求出2L-1个不同的旋转因子，每求出一个旋转因子，就计算完它对应的所有2M-L个蝶形。这样，我们可用三重循环程序实现DIT-FFT运算，程序框图如图 5.3.5 所示。图5.3.5 DIT-FFT运算和程序框图开 始送入 x(n)，MN2 M倒 序L1,MJ0,B 1P2 M LJk J,N1,2
13、LTkXWBkXkXBkXWBkXkXTpNpN)()()()()()(输 出结 束12LB5.序列的倒序DIT-FFT算法的输入序列的排序看起来似乎很乱，但仔细分析就会发现这种倒序是很有规律的。由于N=2M，因此顺序数可用M位二进制数(nM-1 nM-2n1n0)表示。M次偶奇时域抽选过程如图5.3.6所示。第一次按最低位n0的0和1将x(n)分解为偶奇两组，第二次又按次低位n1的0、1值分别对偶奇组分解；依次类推，第M次按nM-1位分解，最后所得二进制倒序数如图5.3.6所示。表5.3.1列出了N=8时以二进制数表示的顺序数和倒序数。图5.3.6 形成例序的树状图(N=23)0101010
14、1010101(n2n1n0)200004261537100010110001101011111表 5.3.1 顺序和倒序二进制数对照表 形成倒序J后，将原存储器中存放的输入序列重新按倒序排列。设原输入序列x(n)先按自然顺序存入数组A中。例如，对N=8,A(0)，A(1)，A(2)，A(7)中依次存放着x(0),x(1),x(7)。对x(n)的重新排序(倒序)规律如图5.3.7所示。倒序的程序框图如图5.3.8所示，图5.3.8中的虚线框内是完成计算倒序值的运算流程图。图 5.3.7 倒序规律x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)A(0)A(1)A(2)A(3)A(
15、4)A(5)A(6)A(7)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)图 5.3.8 倒序程序框图221NNLHJNLHI1,N1I J？T)J(A)J(X)I(A)I(XTJ K？LHK KJJ2KKKJJYNNY5.3.4 DIF-FFT在基2快速算法中，频域抽取法FFT也是一种常用的快速算法，简称DIF-FFT。设序列x(n)长度为N=2M,首先将x(n)前后对半分开，得到两个子序列，其DFT可表示为如下形式：12/012/02/12/0)2/(12/012/12/010)2()()2()()()()(
16、)()(NnknNNnkNNNnNnkNNnknNNNnknNNnknNNnknNWNnxWnxWNnxWnxWnxWnxWnxnxDFTkX式中 WkN/2N=(-1)k=将X(k)分解成偶数组与奇数组，当k取偶数(k=2r,r=0,1,N/2-1)时，1 k=偶数-1 k=奇数rnNNnrnNNnWNnxnxWNnxnxrX22/12/0212/0)2()()2()()2(5.3.9)当k取奇数(k=2r+1,r=0,1,N/2-1)时令nrNnNNnrnNNnWWNnxnxWNnxnxrX2/212/0)12(12/0)2()()2()()12(5.3.10)12/,.2,1,0)2()
17、()()2()()(221NnWNnxnxnxNnxnxnxnN将x1(n)和x2(n)分别代入(5.3.9)式和(5.3.10)式，可得 x1(n)、x2(n)和x(n)的关系也可用图5.3.9所示的蝶形运算流图符号表示。图5.3.10表示N=8时一次分解的运算流图。rnNNnrnNNnWnxrXWnxrX2/12/022/12/01)()12()()2(5.3.11)图 5.3.9 DIF-FFT蝶形运算流图符号x(n)2(Nnx)2()()(1NnxnxnxnNWNnxnxnx)2()()(2nNW图5.3.10 DIF-FFT一次分解运算流图(N=8)N/2点DFTWN0N/2点DFT
18、WN1WN2WN3X(0)x1(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)x1(1)x1(2)x1(3)x2(0)x2(1)x2(2)x2(3)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)图5.3.11示出了N=8时二次分解运算流图。这样继续分解下去，经过M-1次分解，最后分解为2M-1个2点DFT，2点DFT就是一个基本蝶形运算流图。当N=8时，经两次分解，便分解为四个2点DFT,如图5.3.11所示。N=8的完整DIF-FFT运算流图如图5.3.12所示。图5.3.11 DIF-FFT二次分解运算流图(N=8)N/4点DFTWN0WN1WN2WN3x(
19、0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)WN0WN2WN0WN2N/4点DFTN/4点DFTN/4点DFT图5.3.12 DIF-FFT运算流图(N=8)WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0X(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)5.3.5 IDFT的高效算法上述FFT算法流图也可以用于离散傅立叶逆变换(IDFT:Inverse Discrete Fourier Transform)。比较
20、DFT和IDFT的运算公式：由DIF-FFT运算流图改成的DIT-IFFT运算流图如图5.3.13所示。rnNNkrnNNnWkXNnxIDEFnxWnxnxDEFkX1012/0)(1)()()()()(图 5.3.13 DIT-IFFT运算流图WN0WN1WN2WN3WN0WN0N1x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN2WN2N1N1N1N1N1N1N1在实际中，有时为了防止运算过程中发生溢出，将1/N分配到每一级蝶形运算中。由于1/N=(1/2)M,因此每级的每个蝶形输出支路均有一相乘因子1/
21、2，这种运算的蝶形流图如图5.3.14所示。由图可知，乘法次数比图5.3.13增加了(N/2)(M-1)次。图5.3.14 DIT-IFFT运算流图(防止溢出)WN02121x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)212121WN121WN221WN3212121WN021WN2212121WN021WN22121WN02121WN0212121WN021WN021如果希望直接调用FFT子程序计算IFFT，则可用下面的方法：对上式两边同时取共轭，得knNNkknNNkWkXNnxWkXNnx1010)(1)()(1)(由于 因此*)(1)(1)(*10kXDFTNWkXNnxknNNk
