1、1第四章第四章 快速傅立叶变换快速傅立叶变换（FFT）Fast Fourier Transform 数字信号处理2 4.1 引言引言 DFT 是数字信号处理的基础，是数字信号处理的基础，是是 一种重要变换。一种重要变换。快速算法的引出，使数字信号处快速算法的引出，使数字信号处 理技术得到广泛应用。理技术得到广泛应用。各种快速算法不断被采用各种快速算法不断被采用3T)(txa)(txaDA/)(nx)(nxDFT)()(fXkXa数字计算机N足够大.1.1 DFT提供了计算机处理的可能性 模拟信号的频谱分析 4.1.2 DFT的运算量 10)()()(NnnkNWnxnxDFTkXk=0,1,2
2、,N1计算机运算时：0k0)1(0100)1()1()0()0(NNNNWNxWxWxX1k 0111(1)1(1)(0)(1)(1)NNNNXxWxWx NW 2k 0212(1)2(2)(0)(1)(1)NNNNXxWxWx NW 1Nk0111(1)1(1)(0)(1)(1)NNNNNNNX NxWx Wx NW N项 N个 计算一个 N点DFT，共需次复乘。2N以做一次复乘1s计，若N=4096，所需时间为由于计算量大，且要求相当大的内存，难以实现实时处理，限制了DFT的应用。s16777216)4096(2s175 长期以来，人们一直在寻求一种能提高DFT运算速度的方法。FFT便是C
3、ooley&Tukey 在1965 年提出的的快速算法，它可以使运算速度提高几百倍，从而使数字信号处理学科成为一个新兴的应用学科。4.1.3 FFT算法的设计思想 1利用nkNW的特点NjNeW2具有 1）周期性 kNnNnkNWW)()(NknNW2）共轭性 knNNnkNWW)()(3）对称性 kNNkNWW)2(4）恒等性122NNNNWW5）可约性nNnNWW22ReImj8N0NW1NW2NW3NW4NW5NW6NW7NW)(kNnNW2把N 点DFT化为几组点数较少的DFT运算 N点DFT运算的复乘次数为2N次，若将N点DFT化为2 组，则复乘次数约为22222NN次。74.2 基
4、基 2 抽取抽取FFT 算法算法（the Decimation-In-Time Radix-2 FFT Algorithm）FFT分为两大类：分为两大类：时域抽取时域抽取FFT（Decimation-In-Time FFT，简称，简称DIT-FFT）频域抽取频域抽取FFT（Decimation-In-Frequency FFT，简称简称DIF-FFT）8 4.2.1 直接计算直接计算DFT的问题及改进的途径的问题及改进的途径 k=0,1,N-110)()(NnknNWnxkX10)(1)(NkknNWkXNnxn=0,1,N-1DFT及及IDFT的定义的定义9可见，可见，DFT 与与 IDFT
5、 的计算成本基本相同。的计算成本基本相同。直接计算直接计算N点点DFT 时：时：对应一个对应一个k需要需要N次复数乘和（次复数乘和（N-1）次）次复数加；复数加；对所有对所有N个个k值，则需要值，则需要 N复数乘和复数乘和N(N-1)次复数加次复数加。其中其中:一次复数乘需要一次复数乘需要4次实数乘和次实数乘和2次实数加方能次实数加方能完成。当完成。当N较大时，运算量很大。较大时，运算量很大。10 例如：当例如：当N=8时，时，DFT需要需要64次复数乘；次复数乘；当当N=1024时，时，DFT所需复数乘为所需复数乘为1048576次，即一百多万次复数乘。次，即一百多万次复数乘。为了减少运算次
6、数，为了减少运算次数，改进计算的方法：改进计算的方法：1）利用旋转因子的对称性、周期性和可约性；）利用旋转因子的对称性、周期性和可约性；旋转因子旋转因子（twiddle factor）2）使长序列变短。）使长序列变短。NjNeW/2114.2.2 基基2时域抽取法原理时域抽取法原理 （库利（库利-图基算法）图基算法）设序列设序列x（n）的长度为）的长度为N，且，且M为自然数为自然数N-point DFT,N is evenMN2)(log2NM 1,.,1,0,)()(10NkWnxkXNnknN12将其一分为二，分成偶数和奇数序列项将其一分为二，分成偶数和奇数序列项（the even-ind
7、exed and odd-indexed terms）则则N/2点的序列为：点的序列为：even:x1(r)=x(2r),r=0,1,N/2-1 odd:x2(r)=x(2r+1),r=0,1,N/2-113偶数序列偶数序列 the range:02nN-2 (N is even)0n(N/2)-1奇数序列奇数序列 the range:12n+1N-1 (N is even)02nN-2 0n(N/2)-114 奇数偶数nknNnknNWnxWnxkX)()()(120)12(120)2()12()2(NrrkNNrrkNWrxWrx1202212021)()(NrkrNkNNrkrNWrxW
8、Wrx则则x（n）的）的DFT：15由于由于所以所以krNkrNjkrNjkrNWeeW2/2/2222/2 1/2 112/2/20012()()()()()NNkrkkrNNNrrkNX kx r WWx r WX kW Xkk=0,1,N-116n上式说明,按n的奇偶性将 分解为两个N/2长的序列 和 ,则N点DFT可分解为两个N/2点的DFT来计算.)(nx)(1nx)(2nx17 其中其中 N/2-point DFT:k=0,1,N/2-1 )()()(112/02/11rxDFTWrxkXNrkrN12/022/22)()()(NrkrNrxDFTWrxkXk=0,1,N/2-11
9、8因此，可写出两个因此，可写出两个N/2点的方程点的方程:)()()(21kXWkXkXkN)2()2()2(2)2(1NkXWNkXNkXNkN12,.,1,0Nk19 而而120)2(2)()2(11NrrNkNWrxNkXkNNkNWW222)()()2(1112021kXWrxNkXNrkrN20同理同理而而1222jNNjNNeeW)()2(22kXNkX21所以所以)()()2()()()(2121kXWkXNkXkXWkXkXkNkN12,.1,0Nk22表示上述算法可用蝶形结（表示上述算法可用蝶形结（butterfly）kNW)()(21kXWkXkN)()(21kXWkXkN
10、23 Example:如如N=4 x(n)=x0,x1,x2,x3 even:x0,x2 even:x0 odd:x2 odd:x1,x3 even:x1 odd:x324 bit reversal/shuffling process x x0 0 x x2 2x x1 1x x3 3x x0 0 x x2 2x x1 1x x2 2x x0 0 x x1 1x x2 2x x3 3混序或反混序或反序序码位倒置码位倒置分成四个分成四个1点的序列点的序列 25 the butterfly（蝶形运算）（蝶形运算）04W14W02W02W1点序列的点序列的DFT就是序列本身，即不用计算就是序列本身，
11、即不用计算-1-1-1-126 如如N4，则，则将将 x1(r)再按再按r的奇偶进一步分解成两个的奇偶进一步分解成两个N/4点长的子序列：点长的子序列：)12()()2()(1413lxlxlxlx14,.,1,0Nl2714/0)12(2/114/022/11)12()2()(NllkNNlklNWlxWlxkX12,.,1,0Nk)()(42/3kXWkXkN14/04/42/14/04/3)()(NlklNkNNlklNWlxWWlx28 其中其中)()()(314/04/33lxDFTWlxkXNlklN)()()(414/04/44lxDFTWlxkXNlklN14,.,1,0Nk2
12、9由由X3（k）和）和X4（k）的周期性和）的周期性和 的对称性，得的对称性，得kNNkNWW2/4/2/)()()4()()()(42/3142/31kXWkXNkXkXWkXkXkNkN14,.,1,0Nk30同理，得同理，得)()()4()()()(62/5262/52kXWkXNkXkXWkXkXkNkN14,.,1,0Nk31 其中其中)12()()2()(2625lxlxlxlx)()()(514/04/55lxDFTWlXkXNlklN)()()(614/04/66lxDFTWlXkXNlklN32 8点点DIT-FFT运算流图运算流图x(0)x(0)x(4)x(4)x(2)x(
13、2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)X X3 3(0)(0)X X3 3(1)(1)X X4 4(0)(0)X X4 4(1)(1)X X5 5(0)(0)X X5 5(1)(1)X X6 6(0)(0)X X6 6(1)(1)X X1 1(0)(0)X X1 1(1)(1)X X1 1(2)(2)X X1 1(3)(3)X X2 2(0)(0)X X2 2(1)(1)X X2 2(2)(2)X X2 2(3)(3)X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)02
14、W02W02W02W04W04W14W14W08W28W18W38W33 4.2.3 IDFT的高效算法的高效算法10)(1)()(NkknNWkXNkXIDFTnx这样这样IFFT可与可与FFT用同一子程序实现。用同一子程序实现。*10*)(1NkknNWkXN*)(1kXDFTN34IDFT的运算规律的运算规律X*(0)X*(0)X*(4)X*(4)X*(2)X*(2)X*(6)X*(6)X*(1)X*(1)X*(5)X*(5)X*(3)X*(3)X*(7)X*(7)x x3 3(0)(0)x x3 3(1)(1)x x4 4(0)(0)x x4 4(1)(1)x x5 5(0)(0)x
15、x5 5(1)(1)x x6 6(0)(0)x x6 6(1)(1)x x1 1(0)(0)x x1 1(1)(1)x x1 1(2)(2)x x1 1(3)(3)x x2 2(0)(0)x x2 2(1)(1)x x2 2(2)(2)x x2 2(3)(3)x(0)*x(0)*x(1)*x(1)*x(2)*x(2)*x(3)*x(3)*x(4)*x(4)*x(5)*x(5)*x(6)*x(6)*x(7)*x(7)*02W02W02W02W04W04W14W14W08W28W18W38W1/N35 求求IFFT，也可用，也可用DIT-FFT的流程来实现。的流程来实现。x(0)x(0)x(4)x
16、(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)0NW1NW2NW3NW0NW0NW0NW0NW0NW0NW2NW2NW1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N1/N36 Example:Determine the x(n)by IDFT.X=5,-1,1,-1 Solution:n=0,1,2,3*304*)(41)(1)(kknkXWXDFTNnx371)(41
17、)(41)0(3210*3004XXXXXWxkk1)15(41)(41)1(*304jjXWxkkk2)(41)2(30*24kkkXWx38 所以所以 x(n)=1,1,2,1 1)(41)3(*3034kkkXWx39%the program in MATLAB:X=input(Please input X=);N=length(X);x=fft(conj(X),N);x=conj(x/N)40 Example:用用FFT算法计算下面信号的算法计算下面信号的 8点点DFT；x(n)=4 3 2 0 1 2 3 1 然后，再用然后，再用 IFFT恢复原信号。恢复原信号。41 solutio
18、n:4 4-3-32 20 0-1-1-2-23 31 14 42 2-1-13 3-3-30 0-2-21 14 4-1-12 23 3-3-3-2-20 01 13 35 55 5-1-1-5-5-1-11 1-1-11 1-j-j1 1-j-j8 85+j5+j-2-25-j5-j-4-4-1+j-1+j-6-6-1-j-1-j1 1(1-j)/(1-j)/-j-j-(1+j)/-(1+j)/4 45+j+j5+j+j-2+j6-2+j65-j+j5-j+j12125+j-j5+j-j-2-j6-2-j65-j-j5-j-jshufflingshufflingDFT mergingDFT
19、 mergingX X0 0X X1 1X X2 2X X3 3X X4 4X X5 5X X6 6X X7 7 2 2 2 2 2 242 IFFT(X)=1/NFFT(X*)*4 45 5-j j-j j-2 2-j j6 65 5+j j-j j1 12 25 5-j j+j j-2 2+j j6 65 5+j j+j js sh hu uf ff fl li in ng gD DF FT T m me er rg gi in ng g4 4-2 2-j j6 61 12 2-2 2+j j6 65 5-j j-j j5 5+j j-j j5 5-j j+j j5 5+j j+j j4
20、41 12 2-2 2-j j6 6-2 2+j j6 65 5-j j-j j5 5-j j+j j5 5+j j-j j5 5+j j+j j1 16 6-8 8-4 4-1 12 2j j1 10 0-j j2 2-j j2 21 10 0+j j2 2-j j2 21 12 2-2 20 02 20 04 42 20 0-2 2 (1 1+j j)-4 4j j2 2 (1 1-j j)3 32 2-2 24 41 16 60 0-8 8-1 16 62 24 48 8 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 204W04W14W14W08W18W28W38W乘以1
21、/N=1/8得原序列43 4.2.4 FFT的计算成本的计算成本最简单的最简单的FFT 是是Cooley-Tukey法法因为因为NMNM2log2N点点DFT有有M级蝶形运算，每一级都有级蝶形运算，每一级都有N/2个个蝶形运算结构；蝶形运算结构；每个蝶形运算结构都有每个蝶形运算结构都有1次复数乘和次复数乘和2次复数次复数加；加；8442222111111118-DFT8-DFT4-DFT4-DFT2-DFT2-DFT1-DFT1-DFTstage 1stage 1stage 2stage 2stage 3stage 345所以，所以，M级蝶形运算有复数乘级蝶形运算有复数乘M级蝶形运算有复数加级
22、蝶形运算有复数加NNMNCM2log22)2(NNMNCA2log)2(46直接直接DFT的复数乘和复数加均近似为的复数乘和复数加均近似为 。当。当N1时，时，所以所以DIT-FFT大大减少了运算次数。大大减少了运算次数。NNN22log)2(2N47 Example:for N=8 Solution:M=3 stages（三级）（三级）for FFT,the total cost is（FFT的总计算成本是）的总计算成本是）MN/2=3 8/2=12 complex multiplications （复数乘）（复数乘）8log248 for directly DFT,the total co
23、st is（直接计算（直接计算DFT的成本是）的成本是）N=8=64(complex multiplications)The ratio is:（比值是）（比值是）实际上，当实际上，当N=2048时，直接运算与时，直接运算与 FFT算法的计算量之比为算法的计算量之比为372.433.5126449 4.2.5 DIT-FFT的运算规律及编程思想的运算规律及编程思想1、原位运算：、原位运算：利用同一单元存储蝶形计算的输入、输出利用同一单元存储蝶形计算的输入、输出数据。数据。每个蝶形的输入和输出均为相同位数。每个蝶形的输入和输出均为相同位数。原位运算可节省大量内存，因而硬件成本低。原位运算可节省大
24、量内存，因而硬件成本低。2、旋转因子的变化规律：、旋转因子的变化规律：N点点DFT，共有，共有M级蝶形运算，每级有级蝶形运算，每级有N/2个蝶形。个蝶形。50 称为旋转因子，称为旋转因子，p称为旋转因子的指数。称为旋转因子的指数。为了编写程序，应找出旋转因子为了编写程序，应找出旋转因子 与运算与运算 级数的关系。级数的关系。设设L=1，2，M，表示从左到右的运算，表示从左到右的运算 级数，第级数，第L级有级有 个不同的旋转因子，个不同的旋转因子，pNW12LpNW51对于对于 ，各级的旋转因子表示为，各级的旋转因子表示为L=1时，时，L=2时，时，L=3时，时，823N0,24/JWWWJJN
25、pNL1,0,22/JWWWJJNpNL3,2,1,0,2JWWWJJNpNL52对于对于 的一般情况，第的一般情况，第L级的旋转因子为级的旋转因子为MN2JpNLWW212,.,2,1,01LJ，由于由于MLMLMLN222253所以所以通过上式，可以计算第通过上式，可以计算第L级运算的旋转因子。级运算的旋转因子。LMMLJNJNpNWWW2212,.,1,01LJLMJp254 3、蝶形运算规律、蝶形运算规律设序列设序列x（n）经过时域倒序存入数组）经过时域倒序存入数组X如蝶形运算的两个输入数据相距如蝶形运算的两个输入数据相距B个点，个点，应用原位运算，可得：应用原位运算，可得：pNLLL
26、pNLLLWBJXJXBJXWBJXJXJX)()()()()()(1111MLJJpLLM,.,1;12,.,1,0;2155 4、编程思想及程序框图、编程思想及程序框图其它运算规律：其它运算规律：第第L级中，每个蝶形的两个输入数据相距级中，每个蝶形的两个输入数据相距 个点；个点；同一旋转因子对应着间隔为同一旋转因子对应着间隔为 点的点的 个蝶形（对应第几组蝶形）个蝶形（对应第几组蝶形）12LBLM2L256 8点点DIT-FFT运算流图运算流图x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)X X3 3(0)(
27、0)X X3 3(1)(1)X X4 4(0)(0)X X4 4(1)(1)X X5 5(0)(0)X X5 5(1)(1)X X6 6(0)(0)X X6 6(1)(1)X X1 1(0)(0)X X1 1(1)(1)X X1 1(2)(2)X X1 1(3)(3)X X2 2(0)(0)X X2 2(1)(1)X X2 2(2)(2)X X2 2(3)(3)X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)02W02W02W02W04W04W14W14W08W28W18W38W57编程的运算方法：编程的运算方法：
28、从输入端（第从输入端（第1级）开始，逐级进行，共级）开始，逐级进行，共进行进行M级运算。级运算。在进行第在进行第L级运算时，依次求出级运算时，依次求出 个不同的个不同的旋转因子，然后计算其对应相同的旋转因子的旋转因子，然后计算其对应相同的旋转因子的 个蝶形。个蝶形。可用三重循环程序实现可用三重循环程序实现DIT-FFT运算。运算。12LLM258(3)595、序列的倒置（、序列的倒置（bit reversal）倒序是有规律的。倒序是有规律的。由于由于 ，所以顺序数可用，所以顺序数可用M位二进制位二进制数（数（）表示。）表示。MN20121.nnnnMM60用硬件电路和汇编语言程序产生倒序很容易
29、，用硬件电路和汇编语言程序产生倒序很容易，用高级语言倒序的规律为：用高级语言倒序的规律为：倒序数是在倒序数是在M位二进制数最高位加位二进制数最高位加1，逢逢2向右进位。向右进位。614.2.6 频率抽取法频率抽取法FFT(DIF-FFT)设序列设序列x（n）长度为）长度为 ，将其前后，将其前后对半分开，得：对半分开，得：MN212/12/010)()()()()(NNnknNNnknNNnknNWnxWnxWnxnxDFTkX62式中12/0)2/(12/0)2()(NnNnkNNnknNWNnxWnx奇数，偶数kkWkkNN1,1)1(2/knNkNNNnWNnxWnx)2()(2/12/0
30、63再将再将X（k）分解成偶数组和奇数组）分解成偶数组和奇数组 k为偶数时：为偶数时：/2 120/2 1/20(2)()()2()()2NrnNnNrnNnNXrx nx nWNx nx nW64k为奇数时：为奇数时：rnNnNNnnrNNnWWNnxnxWNnxnxrX2/12/0)12(12/0)2()()2()()12(65令令12,.,1,0,)2()()()2()()(21NnWNnxnxnxNnxnxnxnN66得得12/02/212/02/1)()12()()2(NnrnNNnrnNWnxrXWnxrX12,.,1,0Nr67 DIF-FFT蝶形运算流图蝶形运算流图x x(n)
31、(n)x x(n+N/2)(n+N/2)x x1 1(n)=x(n)+x(N/2)(n)=x(n)+x(N/2)x x2 2(n)=x(n)-x(n+N/2)(n)=x(n)-x(n+N/2)nNWnNW-+68N=8时，时，DIF-FFT蝶形运算流图蝶形运算流图x(0)x(0)x(4)x(4)x(2)x(2)x(6)x(6)x(1)x(1)x(5)x(5)x(3)x(3)x(7)x(7)X(0)X(0)X(1)X(1)X(2)X(2)X(3)X(3)X(4)X(4)X(5)X(5)X(6)X(6)X(7)X(7)0NW1NW2NW2NW2NW3NW0NW0NW0NW0NW0NW0NW69注意
32、：注意：DIT-FFT和和DIF-FFT的算法流图不的算法流图不 是唯一的。是唯一的。其变形运算流图见其变形运算流图见P108。70 4.3 进一步减少运算量的措施进一步减少运算量的措施以程序的复杂度换取计算量的进一步提高。以程序的复杂度换取计算量的进一步提高。4.3.1 多类蝶形单元运算多类蝶形单元运算第一级旋转因子可简化：第一级旋转因子可简化：第二级旋转因子可简化：第二级旋转因子可简化：称为无关紧要的旋转因子称为无关紧要的旋转因子。10NWjWWNNN4/01和71其复数乘法次数可减少为：其复数乘法次数可减少为：CM（2）=（M-2）*N/2当当L=3时，第三级蝶形有两个无关紧要旋转时，第
33、三级蝶形有两个无关紧要旋转因子，同一因子对应因子，同一因子对应2(M-L)=N/2L级蝴级蝴蝶结蝶结,所以第三级共有所以第三级共有(书书110页页)32*24NN72依次类推，从依次类推，从L=3到到L=M共可减少复数乘共可减少复数乘法次数为：法次数为：DIT-FFT的复数乘法次数为：的复数乘法次数为：22231MLLNN2)3(2)22()2(2)2(MNNMNCM73另外，另外，也可用实数乘法减少计算量。也可用实数乘法减少计算量。包含所有旋转因子称为一类蝶形单元运算；包含所有旋转因子称为一类蝶形单元运算；去掉去掉 为二类；为二类；去掉去掉 为三类；为三类；依次类推，称为多类蝶形单元运算。依
34、次类推，称为多类蝶形单元运算。N=4096时，三类与一类比，仅时，三类与一类比，仅75%。2/2)1(8/jWNN1rNWjWrN744.3.2 旋转因子的生成旋转因子的生成 直接查表，提高速度，多占内存。直接查表，提高速度，多占内存。4.3.3 实序列的实序列的FFT算法算法 两个实序列，构造序列两个实序列，构造序列y（n）1212()()()()()()()epopy nx njx nDFT x nYkDFT jx nYk75由于由于x（n）为实序列，所以）为实序列，所以X（k）具有共轭）具有共轭对称性：对称性：)()()(21kXWkXkXkN12,.,1,0),()(*NkkXkNX76 4.4 分裂基分裂基FFT算法算法任意基：任意基：基较大，则程序或硬件将很复杂，基大基较大，则程序或硬件将很复杂，基大于于8无意义。无意义。分裂基：分裂基：采用基采用基2和基和基4的混合算法，可有效提高的混合算法，可有效提高运算速度。运算速度。77第四章作业n1.课本127页第1题,第4题n2.请用学号的后8位做N=8点的DIT-FFT,DIF-FFT,要求画出求解的各步骤.
