第4章快速傅里叶变换FFT

1、2021-11-12124.1 4.1 引言引言4.24.2 基基2 FFT2 FFT算法算法4.3 4.3 进一步进一步减少运算量的措施减少运算量的措施3一、一、DFT的运算量估计的运算量估计1, 1 , 0,)()(10NkWnxkXNnnkN 101, 1 , 0,)(1)(NknkNNnWkXNnx 两者的差别仅在指数的符号和因子两者的差别仅在指数的符号和因子1/N。 4.1 引言引言41210) 1() 2() 1 () 0() 1 (NNNNNWNxWxWxWxX一个一个X(k)X(k)的值的工作量，如的值的工作量，如X(1) X(1) ：则计算完全部的则计算完全部的X(k)X(k

2、)，k=0,1,2,N-1k=0,1,2,N-1，需要，需要共共需需N N次复数乘法次复数乘法N-1N-1次复数加法次复数加法N NN=N2N=N2次复数乘法次复数乘法N N(N-1)(N-1) N2N2次复数加法次复数加法nkNW需将每一个需将每一个x(n)x(n)乘以相应的乘以相应的 ，再加起来。，再加起来。5如果换算成实数运算：如果换算成实数运算：一次复乘：一次复乘：(a+jb)(c+jd)=(ac-bd)+j(ad+bc)需需4次实乘，次实乘，2次实加次实加一次复加一次复加：(a+jb)+(c+jd)=(a+c)+j(b+d)需需2次实加次实加所以，直接计算全部所以，直接计算全部X(k

3、)共需共需4N2次实乘次实乘2N2+2N2=4N2次实加次实加6例例 N=210=1024，总共需，总共需实时信号处理对实时信号处理对计算速度计算速度有十分苛刻的要求。有十分苛刻的要求。如何能减少如何能减少DFT运算量？运算量？ 400万万次乘法次乘法 400万万次加法次加法7 2. 2. 对称对称性：性：)(*)(kNnNnkNnkNWWW的特征如下：的特征如下：nkNW 1. 周期周期性性：)()(NknNkNnNnkNWWW3. 3. 可约可约性：性：kNkNjkNjkNWeeW224. 4. 正交正交性：性：rNmnrNmnWNNkkmnN,0,1110)(二、提高二、提高DFTDFT

4、运算效率的途径运算效率的途径1428WW如85. 特殊的特殊的:jWWWNNNNN420, 1, 1 如果充分利用复函数如果充分利用复函数 的这些特征，就可以改善的这些特征，就可以改善DFT的运算的运算效率。为了能利用上述特性，需将效率。为了能利用上述特性，需将x(n)或或X(k)这些长序列按一定这些长序列按一定的规律的规律分解成一些短序列分解成一些短序列，即重排，以减少运算次数。，即重排，以减少运算次数。nkNW二、提高二、提高DFT运算效率的途径运算效率的途径的特征如下：的特征如下：nkNW91965年，年，库利库利(cooley)和和图基图基(Tukey)首先提出首先提出FFT算算法。对

5、于法。对于N点点DFT，仅，仅需需(N/2)log2N 次复数乘法运算。次复数乘法运算。例例 N=210= 1024 时：时： 需要的复乘次数需要的复乘次数 (1024/2)log2210=51210=5120次，次， 直接计算需要的复乘次数为直接计算需要的复乘次数为1048576次，次， 5120/1048576=4.88 , 所以利用所以利用FFT算法，速度提高算法，速度提高200多倍。多倍。 1969年，用一组2048点的傅里叶变换分析地震信号时，DFT需要13.5小时；而FFT仅用了2.4秒。10一、基一、基(Radix)(Radix)基基2 FFT2 FFT算法的思想就是将长序列划为

6、短序列，短到什么程度算法的思想就是将长序列划为短序列，短到什么程度？用？用“基（基（RadixRadix）”来表示，即来表示，即基是基是FFTFFT中用到的最小中用到的最小DFTDFT单单元。元。 Radix-2Radix-2（基（基2 2），代表），代表2 2点点DFTDFT。2 2点点DFTDFT实际上只是加减运算。实际上只是加减运算。4.24.2 基基2 FFT2 FFT算法算法112点点DFT02021002) 1 ()0()()0(WxWxWnxXn1 , 0,)()(102kWnxkXnnk 1202102) 1 ()0()() 1 (WxWxWnxXnn即：即：12画出画出流图流

7、图：) 0 ( x) 1 ( x) 0 (X) 1 (X1) 1 () 0 () 1 () 1 () 0 () 0 () 1 () 0 () 1 () 0 (12020202xxXxxXxxWWWWXX2点点DFT基基2FFT算法分为两类：算法分为两类： 时域时域抽取法抽取法FFT（DIT-FFT,Decimation-In-Time FFT)，它是它是库利库利-图基算法图基算法； 频域频域抽取法抽取法FFT（DIF-FFT,Decimation-In-Frequency FFT)。二、二、DIT-FFT算法原理算法原理14 设序列设序列x(n)的长度为：的长度为：N=2M（0,1,2,)，按

8、按n的奇偶的奇偶把把x(n)分解为两个分解为两个 N/2 的子序列。的子序列。(一）一）N/2点点DFT15),5(),3(),1 (1, 1 , 0),() 12(22xxxrrxrxN即 n为偶数为偶数时：时：),4(),2(),0(1, 1 , 0),()2(21xxxrrxrxN即 n为奇数为奇数时：时：1. N/2点分解点分解161010)()(NnNnnkNnkNWnxWnx( (n为偶数为偶数) ) ( (n n为奇数为奇数) )10)()()(NnnkNWnxnxDFTkX10)12(10222) 12()2(NNrkrNrrkNWrxWrx1022102122)()(NNrr

9、kNkNrrkNWrxWWrx)()(21kXWkXkN1, 2 , 1 , 0Nk17 (1) X1(k)，X2(k)均为均为N/2点的点的DFT。(2) 只能确定出只能确定出X(k)的的的的值，值，即即前一半前一半的结果。的结果。1, 1 , 02Nk)()()(21kXWkXkXkN10102210101122222222)12()()()2()()(NNNNNNNNrrkrrkrrkrrkWrxWrxkXWrxWrxkX其中其中2.结论结论12, 2 , 1 , 0Nk)()()(21kXWkXkXkN1, 2 , 1 , 0Nk18同理，同理，3. X(k)后一半的确定后一半的确定由

10、于由于 （周期性周期性），所以：），所以：rkkrNNNWW222)()()()()2(1101)(101122222kXWrxWrxkNXNNNNNrrkkrr)()2(22kXkNXX1(k)、X2(k)后一半后一半的值，分别等于其的值，分别等于其前一半前一半的值。的值。12, 2 , 1 , 0Nk19 可见，可见，X(k)的后一半，也完全由的后一半，也完全由X1(k)、X2(k)的前一半的前一半所确定。所确定。)2()2()2(212NkXWNkXNkXNkN1, 1 , 0),()(221NkNkkXWkX kNkNNkNWWWWNN22)(又又由于由于 ，所以：，所以： 前一半 后

11、一半)()()()()()(2121kXWkXkXkXWkXkXkNkN) 1,() 1, 1 , 0(22NkkNN20 只要求出只要求出0N/2-1区间内的各个区间内的各个k值所对应的值所对应的X1(k)、X2(k)的值，即可以求出的值，即可以求出0N-1整个区间内的全部整个区间内的全部X(k) 值值。这就是。这就是FFT能大量节省运算量的关键。能大量节省运算量的关键。N点的DFT可由两个N/2点的DFT来计算。 前一半 后一半)()()()()()(2121kXWkXkXkXWkXkXkNkN) 1,() 1, 1 , 0(22NkkNN214.4.蝶形运算蝶形运算碟形运算符号碟形运算符

12、号ABCBACBAC22实现上式运算，共需N/2个蝶形运算。 前一半 后一半)()()()()()(2121kXWkXkXkXWkXkXkNkN) 1,() 1, 1 , 0(22NkkNN 由由X1(k)、X 2(k)表示表示X(k)的运算是一种碟形运算：的运算是一种碟形运算：)()()(21kXWkXkXkN)()()2(21kXWkXkNXkN)(1kX)(2kXkNW23(1)N/2(1)N/2点的点的DFTDFT运算量运算量: :复乘次数，复加次数复乘次数，复加次数(2)(2)两个两个N/2N/2点的点的DFTDFT运算量运算量: :复乘次数，复加次数复乘次数，复加次数 (3)N/2

13、(3)N/2个蝶形运算的运算量个蝶形运算的运算量: :复乘次数，复加次数复乘次数，复加次数 5. .运算量分析运算量分析4)2(22NN)12(2NN22N)12(NN2NNN22按奇、偶分组后的计算量：按奇、偶分组后的计算量：24复乘复乘: :复加复加: :2)12(2NNNN22222NNN 但是，但是，N N点点DFTDFT的复乘为的复乘为N N2 2，复加，复加N(N-1)N(N-1)；与分解后相比，计；与分解后相比，计算工作量差不多算工作量差不多减少一半减少一半。5. 5.运算量分析运算量分析一个一个N N点点DFTDFT分解成两个分解成两个N/2N/2点点DFTDFT后的总共运算量

14、：后的总共运算量：25例：例：N=8N=8时的时的DFTDFT，可以分解为两个，可以分解为两个N/2=4N/2=4点的点的DFTDFT。具体。具体方法如下：方法如下： (1)(1)n n为偶数为偶数时时, ,即即分别记作分别记作: :3 , 2 , 1 , 0)2 ()()(30430411kWrxWrxkXrrkrrk，)6()3(),4()2(),2()1(),0()0(1111xxxxxxxx)6(),4(),2(),0(xxxx进行进行N/2=4点的点的DFT，得：，得：26(2)(2)n n为奇数为奇数时时，分别记作：，分别记作：)7()3(),5()2(),3()1(),1()0(

15、2222xxxxxxxx3 , 2 , 1 , 0) 12()()(30430422kWrxWrxkXrrkrrk，)()() 4()()()(2121kXWkXkXkXWkXkXkNkN进行进行N/2=4N/2=4点的点的DFTDFT，得：，得：3 , 2 , 1 , 0k27 x x1 1(0)=(0)=x x(0) (0) x x1 1(1)=(1)=x x(2) (2) N/2N/2点点 x x1 1(2)=(2)=x x(4) DFT (4) DFT x x1 1(3)=(3)=x x(6) (6) x x2 2(0)=(0)=x x(1) (1) x x2 2(1)=(1)=x x

16、(3) (3) N/2N/2点点 x x2 2(2)=(2)=x x(5) (5) DFT DFT x x2 2(3)=(3)=x x(7)(7) X X1(0)(0)X X1(1)(1)X X1(2)(2)X X1(3)(3)X X2(0)(0)X X2(1)(1)X X2(2)(2)X X2(3)(3)WN2WN1WN0WN3X(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)(3)(3)对对X1(k)和和X2(k)进行蝶形运算进行蝶形运算, ,前半部分为前半部分为X(0)X(0)X(3),X(3),后半部分为后半

17、部分为X(4)X(4)X(7)X(7)。 28 由于由于N=2N=2MM , ,所以所以N/2N/2仍为偶数仍为偶数, ,可以进一步把每个可以进一步把每个N/2N/2点的点的序列再按其奇偶部分分解为两个序列再按其奇偶部分分解为两个N/4N/4的子序列。的子序列。例如，例如，n n为为偶数时的偶数时的 N/2N/2点分解为点分解为: :( (二二) ) N/4点点DFT1, 1 , 0),4()2()(413Nllxlxlx1, 1 , 0),24() 12()(414Nllxlxlx进行进行N/4N/4点的点的DFTDFT，得：，得：klNllkNllkNllkNlWlxWlxkXWlxWlx

18、kXNNNN)12(2/1014/104422/1014/1033) 12()()()2()()(4444( (偶中偶偶中偶) )( (偶中奇偶中奇) )29)()()(4312kXWkXkXkN1, 1 , 04Nk从而可得到前从而可得到前N/4的的X1(k)：)()()4(4312kXWkXkNXkN后后N/4的的X1(k)：1, 1 , 04Nk( (二二) ) N/4点点DFT即对即对X3(k)、 X4(k)进行一次碟形运算。进行一次碟形运算。30( (奇中偶奇中偶) )104/64/1026104/5104/254444)()12()()()2()(NNNNllkNlkNlllkNl

19、lkNWlxWlxkXWlxWlxkX( (奇中奇奇中奇) )同样对同样对n n为奇数时，为奇数时，N/2N/2点分为两个点分为两个N/4N/4点点 的序列进行的序列进行DFTDFT，则有：，则有： (k)XW(k) X(k) X6kN/252 (k)XW(k) Xk)4N( X6kN/2521, 1 , 0),14()2()(425Nllxlxlx1, 1 , 0),34() 12()(426Nllxlxlx14, 1 , 0Nk由由X5(k)、 X6(k)进行碟形运算，得：进行碟形运算，得：31 所以，所以，N=8N=8时的时的DFTDFT可分解为四个可分解为四个N/4N/4的的DFT,D

20、FT, 具体步骤如下具体步骤如下: :1 , 0),4()2()(13llxlxlx(1) (1) 原序列原序列x(n)的的“偶中偶偶中偶”部分：部分：构成构成N/4点点DFT，从而得到，从而得到X3(0)，X3(1)。0,1k ,)()(104/33llkNWlxkX)4()0()1 ()0()1 ()4()0()1 ()0()0(031233030233xWxxWxXxWxxWxXNN)4()2()1()0()0()0(1313xxxxxx321 , 0),24() 12()(14llxlxlx(2)(2) 原序列原序列x(n)的的“偶中奇偶中奇”部分：部分：构成构成N/4点点DFT，从而

21、得到，从而得到X(0)，X(1)。0,1k ,)()(104/44llkNWlxkX)6()2() 1 ()0() 1 ()6()2() 1 ()0()0(041244040244xWxxWxXxWxxWxXNN)6()3()1()2()1()0(1414xxxxxx33(3)(3) 原序列原序列x(n)的的“奇中偶奇中偶”部分：部分：1 , 0),14()2()(25llxlxlx构成构成N/4点点DFT，从而得到，从而得到X(0)，X(1)。0,1k ,)()(104/55llkNWlxkX)5()2()1 ()1 ()0()0(2525xxxxxx)5()1 ()1 ()0()1 ()5

22、()1 ()1 ()0()0(051255050255xWxxWxXxWxxWxXNN34(4)(4) 原序列原序列x(n)的的“奇中奇奇中奇”部分：部分：1 , 0),34() 12()(26llxlxlx构成构成N/4点点DFT，从而得到，从而得到X(0)，X(1)。0,1k ,)()(4/1066lkNlWlxkX)7() 3() 1 ()0() 1 ()7() 3() 1 ()0()0(061266060266xWxxWxXxWxxWxXNN)7()3()1()3()1()0(2626xxxxxx35(5)(5)由由X3(0), X4(0)进行碟形运算，进行碟形运算，得到得到X1(0)

23、, X1(2)； 由由X3(1), X4(1)进行碟形运算进行碟形运算，得到得到X1(1), X1(3)。)0()0()2()0()0()0(402/31402/31XWXXXWXXNN) 1 () 1 ()3() 1 () 1 () 1 (412/31412/31XWXXXWXXNN36(6)(6)由由X5(0), X6(0)进行碟形运算，进行碟形运算，得到得到X2(0), X2(2)； 由由X5(1), X6(1)进行碟形运算进行碟形运算，得到得到X2(1), X2(3)。)0()0()2()0()0()0(602/52602/52XWXXXWXXNN) 1 () 1 ()3() 1 ()

24、 1 () 1 (612/52612/52XWXXXWXXNN37（7 7）由）由X1(0), 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) 。)0()0()4()0()0()0(201201XWXXXWXXNN) 1 () 1 ()5() 1 () 1 () 1 (211211XWXXXWXXNN)2()2()6()2()2()2(221221XWXXXWXXNN)3()3()7()3()3()3(231231XWXXX

25、WXXNN38x3(0)=x1(0)=x(0) N/4N/4点点x3(1)=x1(2)=x(4) DFTDFTx4(0)=x1(1)=x(2) N/4N/4点点x4(1)=x1(3)=x(6) DFTDFTx5(0)=x2(0)=x(1) N/4N/4点点x5(1)=x2(2)=x(5) DFTDFT x6(0)=x2(1)=x(3) N/4N/4点点x6(1)=x2(3)=x(7) DFTDFT02NN02NNWWWW0123NNNNWWWWX (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566X (0)X (1)X (2)X (3)X (0)X (

26、1)X (2)X (3)11122221X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)xx N=8N=8时的DFTDFT可分解为四个N/4N/4点点的DFT:DFT:N点DFT的第二次时域抽取分解图(N=8) N/4点DFTWN12WN12WN0WN1WN2WN3X1(0)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)x(0)X3(0)X3(1)X4(0)X4(1)x(4)x(2)x(6)x(1)x(5)x(3)x(7)N/4点DFTN/4点DFTN/4点DFTWN02WN02)0(5X)

27、 1 (5X)0(6X) 1 (6X40 x(0) x(4)x(2) x(6)x(1)x(5)x(3)x(7)02NN02NNWWWW0123NNNNWWWWX (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566X (0)X (1)X (2)X (3)X (0)X (1)X (2)X (3)11122221X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)对于对于N=8N=8时的时的DFTDFT，N/4N/4点即为两点点即为两点DFTDFT。因此，因此，8 8点点DFTDFT的的DIT-FFTDIT-FFT的运算流图的运算流图如下：如下：

28、0NW0NW0NW0NW41 这种这种FFT算法，是在时间上对输入序算法，是在时间上对输入序列的次序是属于偶数还是属于奇数来进行分列的次序是属于偶数还是属于奇数来进行分 解的，所以称作按解的，所以称作按时域时域抽取的算法（抽取的算法（DIT）。）。（三）运算量分析（三）运算量分析 由上述分析可知，由上述分析可知，N=8N=8需三级蝶形运需三级蝶形运算算 N=2N=23 3=8,=8,由此可知，由此可知，N=2N=2M M共需共需M M级蝶形运级蝶形运算算。而且。而且每级都由每级都由N/2N/2个蝶形运算个蝶形运算组成，组成，每个蝶形运算有一次复乘，两次复加。每个蝶形运算有一次复乘，两次复加。4

29、2 因此因此,N,N点的点的DIT-FFTDIT-FFT的运算量为的运算量为 复乘复乘: : (N/2)(N/2)M=(N/2)logM=(N/2)log2 2N N 复加复加：N N=Nlog=Nlog2 2N NDIT-FFTDIT-FFT算法比直接计算算法比直接计算DFTDFT运算次数要大大减少。运算次数要大大减少。如如N=1024N=1024时：时：8.2041051210241024log222NNNFFTDITDFT：直接计算复乘即运算效率提高即运算效率提高200200多倍，原来需要多倍，原来需要3 3个多小时个多小时计算完，现在只需计算完，现在只需1 1分钟分钟即可。即可。43三

30、、三、DIT-FFTDIT-FFT算法的特点算法的特点1.1.原位运算原位运算（in-placein-place）输入数据、中间运算结果和输出数据均用同一存储器。输入数据、中间运算结果和输出数据均用同一存储器。.) 0 ( x0NW0NW0NW0NW0NW0NW2NW0NW2NW1NW2NW3NW) 4 ( x) 2 ( x) 6 ( x) 1 ( x) 5 ( x) 3 ( x) 7 ( x) 0 (1X) 1 (1X) 2 (1X) 3 (1X) 4 (1X) 5 (1X) 6 (1X) 7 (1X) 0 (2X) 1 (2X) 2 (2X) 3 (2X) 4 (2X) 5 (2X) 6

31、(2X) 7 (2X) 0 (0X) 1 (0X) 2 (0X) 3 (0X) 4 (0X) 5 (0X) 6 (0X) 7 (0X) 0 (3X) 1 (3X) 2 (3X) 3 (3X) 4 (3X) 5 (3X) 6 (3X) 7 (3X) 0 (X) 1 (X) 2 (X) 3 (X) 4 (X) 5 (X) 6 (X) 7 (X44rNmmmrNmmmWjXkXjXWjXkXkX)()()()()()(1111),3()6(),2()2(),1 ()4(),0()0(0000XxXxXxXx).7()7(),6()3(),5()5(),4()1(0000XxXxXxXx 设用设用m(

32、m=1,2,m(m=1,2, ,L) ,L)表示第表示第m m级，用级，用k k、j j表示蝶形输入数据所在的上表示蝶形输入数据所在的上/ /下行号下行号(0,1,2,(0,1,2, ,N-1),N-1)，则任何一个蝶形运算可用下式表示：，则任何一个蝶形运算可用下式表示： 由运算流图可知：共有由运算流图可知：共有N个输入个输入/输出行输出行，共有共有log2N=L级蝶形运算级蝶形运算。 45 当当m=1时，则有（前两个蝶形）：时，则有（前两个蝶形）： 00010001)1()0()1()1()0()0(NNWXXXWXXX00010001)3()2()3()3()2()2(NNWXXXWXXX

33、46 当当m=2时时，则有（前两个蝶形）：，则有（前两个蝶形）： 当当m=3时时，则有（前两个蝶形）：，则有（前两个蝶形）： 2112211201120112)3()1()3()3()1()1()2()0()2()2()0()0(NNNNWXXXWXXXWXXXWXXX1223122302230223)5()1()5()5()1()1()4()0()4()4()0()0(NNNNWXXXWXXXWXXXWXXX47可见，由某级进行蝶形运算的可见，由某级进行蝶形运算的任意两行任意两行k和和j的的节点变量，就完全可以确定蝶形运算的结果，节点变量，就完全可以确定蝶形运算的结果，与其它行无关与其它行无

34、关。所以，蝶形运算的两个输出值。所以，蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中。仍可放回蝶形运算的两个输入所在的存储器中。这样，从第一级的输入这样，从第一级的输入x(n)开始，到最后一级开始，到最后一级输出输出X(k)，只需要，只需要N个存储单元，中间无需其个存储单元，中间无需其它存储器，这一规律就称为它存储器，这一规律就称为原位运算原位运算。原位运算大大节省了存储单元。原位运算大大节省了存储单元。 48由运算流图可知，输出由运算流图可知，输出X(k)是按自然是按自然顺序排列在存储单元，而输入顺序排列在存储单元，而输入x(n)(n)是按以是按以下顺序：下顺序： 这种顺序称作

35、这种顺序称作倒位序倒位序，即二进制数，即二进制数倒位。倒位。 倒位序是由倒位序是由奇偶分组奇偶分组造成的。造成的。2. 2. 倒位序规律倒位序规律）7(),3(),5(),1 (,6(),2(),4(),0(xxxxxxxx49n =00n =10n =01n =11n =01n =1101010101 ),(012nnnx(n2)x(000) 0 x(100) 4 x(010) 2 x(110) 6 x(001) 1x(101) 5x(011) 3x(111) 7 （偶）（偶）（奇）（奇）以以N=8N=8为例说明如下：为例说明如下：50 0 00 0 0 0 0 00 0 0 0 0 0 0

36、 0 1 01 0 0 0 1 11 1 0 0 0 4 0 4 2 02 0 1 1 0 00 0 1 1 0 2 0 2 3 03 0 1 1 1 11 1 1 1 0 6 0 6 4 14 1 0 0 0 00 0 0 0 1 1 1 1 5 15 1 0 0 1 11 1 0 0 1 5 1 5 6 16 1 1 1 0 00 0 1 1 1 3 1 3 7 17 1 1 1 1 11 1 1 1 1 7 1 7 自然顺序n n 二进制n n n n n n 倒位序二进制n n n n n n 倒位顺序 n n2 1 0 0 1 2N=8N=8时的顺序和倒序对照表时的顺序和倒序对照表

37、自然顺序加自然顺序加1 1，是在二进制数最低位加，是在二进制数最低位加1 1，并逢，并逢2 2向左进位；向左进位；倒位序加倒位序加1 1，则是在二进制数最高位加，则是在二进制数最高位加1 1，并逢，并逢2 2向右进位向右进位。51倒位序加倒位序加1 1二进制数最高位加二进制数最高位加1 1，并逢，并逢2 2向右进位向右进位nM位二进制数：从高到低，位二进制数：从高到低，各位的权值各位的权值分别为分别为2M-1、2M-2、21、20。设。设N=2M ，则，则权值分别为权值分别为N/2、 N/4、2、1。n最高位加最高位加1，相当于，相当于十进制加十进制加N/2的运算。的运算。n若最高位为若最高位

38、为0( N/2)，则直接加，则直接加N/2； 若最高位为若最高位为1 ( N/2) ，则向右进位，即最高位变为，则向右进位，即最高位变为0，然后次高位加然后次高位加N/4。52 3. 3.倒位序实现倒位序实现 u输入序列先按自然顺序存入存储单元，然后输入序列先按自然顺序存入存储单元，然后经经变址运算变址运算来实现来实现倒位序排列倒位序排列。u变址处理方法：变址处理方法：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)x(0) x(4) x(2) x(6) x(1) x(5) x(3) x

39、(7)存储单元存储单元自然顺序自然顺序变址变址倒序倒序n自然序号自然序号I与与倒序倒序号号J相等时，不交换。相等时，不交换。n为避免交换二次，只考虑为避免交换二次，只考虑I小于小于J的情况的情况。53倒序程序框图倒序程序框图544.4.蝶形运算两点的距离蝶形运算两点的距离 2 2m-1m-1 （m m表示第表示第m m级蝶形运算）级蝶形运算）例如例如N=8=2N=8=23 3, ,则第一级距离为则第一级距离为2 21-11-1=1,=1,第二级距离为第二级距离为2 22-12-1=2=2，第三级距离为，第三级距离为2 23-13-1=4=4。rNmmmrNmmmWjXkXjXWjXkXkX)(

40、)()()()()(1111即：即： j = k + 2m-1555.W5.WN Nr r 的确定的确定u考虑蝶形运算两点的距离为考虑蝶形运算两点的距离为2 2m-1m-1, ,蝶形运算表示为蝶形运算表示为 Xm(k)=Xm-1(k)+Xm-1(k+2m-1)WNr Xm(k+2m-1)=Xm-1(k)-Xm-1(k+2m-1)WNr 由于由于N N为已知为已知, ,所以确定所以确定r r的值即可。的值即可。u设共有设共有L L级蝶形运算，则级蝶形运算，则(r)(r)2 2 =(k)k)2 22 2L-mL-m 。即：即：令令k=(nk=(n2 2n n1 1n n0 0) )2 2 , ,再

41、将再将k= (nk= (n2 2n n1 1n n0 0) )2 2 左移左移(L-m)(L-m)位位, ,右右边位置补零边位置补零, ,就可得到就可得到(r)(r)2 2 的值。的值。566.6.存储单元存储单元l例：例：N=8 FFTN=8 FFT运算，输入：运算，输入：x(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)A(0)A(1)A(2)A(3)A(4)A(5)A(6)A(7)A(0)=X(0)A(1)=X(1)A(2)=X(2)A(3)=X(

42、3)A(4)=X(4)A(5)=X(5)A(6)=X(6)A(7)=X(7)012388881,23,4RWRW RWRW暂存器R1R1R1R1R1R3R1R1R3R2R3R457则一个完整的则一个完整的N点点FFT运算只需要：运算只需要：（1）N个存储单元来存放原始序列个存储单元来存放原始序列x(n)(n),n=0,1,n=0,1,N-1,N-1 ；（2）N/2个存储单元用来存放旋转因子个存储单元用来存放旋转因子WNr , , r=0,1,r=0,1, ,（N/2N/2）-1-1；共计共计(N+N/2)个存储单元。个存储单元。587.7.按时域抽取的按时域抽取的FFT算法的若干变体算法的若干

43、变体 思路思路：对于任何流图，只要保持：对于任何流图，只要保持各节点所各节点所连续的支路连续的支路及其及其传输系数不变传输系数不变，则不论节，则不论节点位置怎么排列所得流图总是等效的，点位置怎么排列所得流图总是等效的，最最后所得结果都是后所得结果都是x(n)的离散的离散傅傅里叶变换的正里叶变换的正确结果确结果。只是数据的提取和存放的次序不。只是数据的提取和存放的次序不同而已。同而已。59(1)(1)输入是自然顺序而输出是乱序输入是自然顺序而输出是乱序将原先中与将原先中与x(4)水平相邻的所有节点跟水平相邻的所有节点跟x(1)水平相邻的所有节水平相邻的所有节点位置对调；将原先中与点位置对调；将原

44、先中与x(6)水平相邻的所有节点跟水平相邻的所有节点跟x(3)水平水平相邻的所有节点位置对调，其余诸节点保持不变，则可得：相邻的所有节点位置对调，其余诸节点保持不变，则可得：x(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(3)X(5)X(7)它与输入是乱序、它与输入是乱序、输出顺序比较，看输出顺序比较，看出：出：相同点相同点：运算：运算量一样；量一样；不同点不同点：第一是数据存入方第一是数据存入方式不同；第二是取式不同；第二是取用系数的顺序不同。用系数的顺序不同。08W08W08W08W08W28W08W28W08W18W38W28W

45、60(2)(2)输入和输出都是自然顺序输入和输出都是自然顺序对输入自然顺序，输出乱序的第二级，第三级节点作调整，对输入自然顺序，输出乱序的第二级，第三级节点作调整，可得输入输出都是顺序，无需数据重排可得输入输出都是顺序，无需数据重排：x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)（1）它失掉了）它失掉了“原位运算原位运算”的的性质。（性质。（2）为）为了变换了变换N点数据，点数据，至少需要至少需要2N个个复数单元。在内复数单元。在内存比较紧张时，存比较紧张时，而对输入数据整而对输入数据整序并不困难时一序并不困难时

46、一般不用。为了争般不用。为了争取速度，才采取取速度，才采取这种变体。这种变体。08W08W08W08W08W28W08W28W08W18W28W38W61三、三、DIF-FFT算法原理算法原理频域抽取法基频域抽取法基2 2FFT算法算法（Radix-2 DIF-FFT,Decimation-In-Frequency FFT)，基本思基本思想想是在时域下对输入是在时域下对输入x(n)前后分组前后分组，输出，输出X(k)按其按其奇偶顺序奇偶顺序进行分组。进行分组。 设序列设序列x(n)的长度为：的长度为：N=2N=2M M（0,1,2,0,1,2,) )，按按n n的前后的前后把把x(n)分解为两

47、个分解为两个 N/2 N/2 的子序列。的子序列。1. N/2点分解点分解6212012)()(NnNNnnkNnkNWnxWnx( (前一半前一半) ) ( (后一半后一半) )10)()()(NnnkNWnxnxDFTkX10)2/(1022)2()(NNnkNnNnnkNWNnxWnx102/2)2()(NnnkNkNNWNnxWnx102)2() 1()(NnnkNkWNnxnx63K=2mK=2m(m=0,1,(m=0,1, ,N/2-1),N/2-1)时：时：1. N/2点分解点分解 1022)2()()2(NnnmNWNnxnxmX102/12)(NnnmNWnxK=2m+1K=

48、2m+1(m=0,1,(m=0,1, ,N/2-1),N/2-1)时：时：10)12(2)2()() 12(NnmnNWNnxnxmX102/22)(NnnmNWnx1022)2()(NnnmNnNWWNnxnx6412, 2 , 1 , 0Nn(1) x1(n)表示序列表示序列x(n)前前N/2点值与后点值与后N/2点值点值之和之和，x2(n)表示序列表示序列x(n)前前N/2点值与后点值与后N/2点值的差再与序列点值的差再与序列W WN Nn n之积之积 。(2) X(2m)为为x1(n)的的N/2点点DFT，X(2m+1)为为x2(n)的的N/2点点DFT。(3) N点的点的DFT可由两

49、个可由两个N/2点的点的DFT来计算。来计算。2.2.三三点结论点结论)2()()(1Nnxnxnx其中其中nNWNnxnxnx)2()()(265实现上式运算需实现上式运算需N/2个个DIF-FFT蝶形运算。蝶形运算。. .DIF-FFT蝶形运算蝶形运算)(nx)2(Nnx12, 2 , 1 , 0Nn)2()()(1NnxnxnxnNWNnxnxnx)2()()(2)2()()(1NnxnxnxnNWNnxnxnx)2()()(2WNn-1-1kNW66DIF-FFT碟形运算符号：碟形运算符号：ABBACBA)(C. .DIF-FFT蝶形运算蝶形运算由于由于N=2M，所以，所以N/2仍为偶

50、数。可以把仍为偶数。可以把N/2点点DFT继续奇、偶分组，这样就把一个继续奇、偶分组，这样就把一个N/2点的点的DFT分解分成两个分解分成两个N/4点的点的DFT。 如此下去，对于如此下去，对于N=2M点点DFT经过经过(M-1)次分解，最后分解成次分解，最后分解成2M-1个两点的个两点的DFT，每个两点的，每个两点的DFT就可以用一个就可以用一个DIF-FFT蝶形蝶形运算来表示，这样就形成了完整的运算来表示，这样就形成了完整的N点点DIF-FFT运运算流图。算流图。 （P120：图：图4.2.13）4.4.DIF-FFT运算流图运算流图675. .运算量分析运算量分析u当当N=2M时，经过时

51、，经过(M-1)(M-1)次分解，整个次分解，整个DIF-FFT运算有运算有M M级蝶形，每级蝶形有级蝶形，每级蝶形有N/2N/2个蝶形运算，每个蝶形运算有一次复乘和个蝶形运算，每个蝶形运算有一次复乘和两次复加。两次复加。u因此，因此，N N点的点的DIF-FFTDIF-FFT的运算量为的运算量为 复乘复乘: : (N/2)M=(N/2)log2N 复加复加：N=Nlog2NDIF-FFT算法运算量与算法运算量与DIT-FFT完全相同。完全相同。68四、四、DIT-DIT-FFT与与DIF-DIF-FFT算法的比较算法的比较u相同点相同点：（1 1）两种算法均为两种算法均为原位运算原位运算。（

52、2 2）两种算法两种算法的的运算量相同运算量相同。故。故DITDIT与与DIFDIF是两是两种等价的种等价的FFTFFT算法。算法。u不同点不同点：（1 1）两种算法两种算法结构倒过来结构倒过来：DIFDIF为输入顺序，为输入顺序，输出倒序，即输出倒序，即运算完毕再运行运算完毕再运行“倒倒位序位序”程序程序；而而DITDIT为输入倒序，输出顺序，即为输入倒序，输出顺序，即运行运行“倒倒位序位序”程序再运算程序再运算。故。故DITDIT与与DIFDIF互为转置。互为转置。（2 2）两种算法两种算法的的根本区别在于根本区别在于蝶形结不同蝶形结不同： DITDIT蝶形是先乘旋转因子，后进行加减运算；

53、而蝶形是先乘旋转因子，后进行加减运算；而DIFDIF蝶形是先进行加减运算，后乘以旋转因子。蝶形是先进行加减运算，后乘以旋转因子。69五、五、IDFTIDFT的的快速快速算法算法（IFFT）1, 1 , 0 ,)()(10NkWnxkXNnnkN101, 1 , 0 ,)(1)(NknkNNnWkXNnx比较比较IDFTIDFT与与DFTDFT的公式：的公式： 可以直接调用可以直接调用FFTFFT的程序来计算的程序来计算IDFTIDFT。70n 利用利用FFTFFT计算计算IFFTIFFT时在命名上应注意时在命名上应注意：（1 1）把把FFTFFT的时的时域域抽取法，用于抽取法，用于I IF F

54、FTFT运算时，由于运算时，由于输入变量输入变量由由时间序列时间序列x(n)x(n)改成改成频率序列频率序列X(k)X(k)，原来按原来按x(n)x(n)的奇、偶次序的奇、偶次序分组的分组的时时域域抽取法抽取法FFTFFT，现在就变成了按，现在就变成了按X(k)X(k)的奇偶次序抽的奇偶次序抽取取的的频域抽取法频域抽取法I IF FFTFT。（2 2）同样，）同样，频频域域抽取抽取法法FFTFFT运算用于运算用于I IF FFTFT运算时，也应改变运算时，也应改变为为时时域域抽取抽取法法IFFTIFFT。71IFFTIFFT的基本蝶形运算的基本蝶形运算21nNW21BWAnN21BWAnN21

55、AB21nNW21BA21nNWBA21AB(a)频频域域抽取抽取IFFT的蝶形运算的蝶形运算(b)时域抽时域抽取取IFFT的蝶形运算的蝶形运算l在在IFFTIFFT的运算中，常常把的运算中，常常把1/N1/N分解为分解为(1/2)(1/2)MM，并且在，并且在MM级运算中级运算中每一级运算都分别乘以每一级运算都分别乘以1/21/2因子，就可得到因子，就可得到IFFTIFFT的两种基本蝶的两种基本蝶形运算结构。形运算结构。72先将先将X(k)取共轭，然后直接调用取共轭，然后直接调用FFT程序，再程序，再对运算结果取共轭，并乘以常数对运算结果取共轭，并乘以常数1/N，就可以得，就可以得到时域下的

56、到时域下的x(n)。 方法一方法一：10)(1)(NknkNWkXNnx10)(1NknkNWkXN1DFT( )XkN 73先调用先调用FFT程序计算程序计算X(k)的的DFT，然后把运算，然后把运算结果翻褶后右移结果翻褶后右移N位，最后乘以常数位，最后乘以常数1/N，就可，就可以得到时域下的以得到时域下的x(n)。 方法二方法二： 令令10)()()(NknkNWkXkXDFTng10)()(1)(1)(NknNkNWkXNnNgNnx 则则10)(1NknkNWkXN74一、蝶形单元的特点一、蝶形单元的特点 )2(2MNCM4.3 4.3 进一步减少运算量的措施进一步减少运算量的措施无关

57、紧要的旋转因子无关紧要的旋转因子：值为：值为1和和j因为第一、二级蝶形运算的旋转因子都因为第一、二级蝶形运算的旋转因子都是无关紧要的旋转因子，故是无关紧要的旋转因子，故FFTFFT算法中所需的算法中所需的复乘次数可减少为：复乘次数可减少为：752) 3(2)22()2(2MNNMNCM再考虑其它级蝶形运算中的无关紧要的再考虑其它级蝶形运算中的无关紧要的旋转因子，则可减少的复乘次数为：旋转因子，则可减少的复乘次数为：22231NNMLL所以，引入无关紧要的旋转因子后，所以，引入无关紧要的旋转因子后，FFTFFT算法中的复乘次数可减少为：算法中的复乘次数可减少为：76在基在基2 FFT2 FFT算

58、法中，若包含了所有旋转因算法中，若包含了所有旋转因子，则称该算法为子，则称该算法为一类一类碟形单元运算碟形单元运算；若去；若去掉值为掉值为1的的旋转因子，则称为旋转因子，则称为二类二类碟形单元碟形单元运算运算；若再去掉值为；若再去掉值为j的的旋转因子，则称为旋转因子，则称为三类三类碟形单元运算碟形单元运算。除一类外的其它类运算称为除一类外的其它类运算称为多类多类碟形单碟形单元运算元运算。77二、旋转因子的生成二、旋转因子的生成 为了进一步减少运算量，在为了进一步减少运算量，在FFTFFT程序开程序开始前始前预先计算预先计算出旋转因子出旋转因子WNr的值，生成的值，生成旋旋转因子表转因子表，并存

59、放在，并存放在内存内存中，在程序执行过中，在程序执行过程中可以直接查表得到所需的旋转因子值。程中可以直接查表得到所需的旋转因子值。5.1 5.1 快速傅里叶变换（快速傅里叶变换（FFTFFT）的提出）的提出 DFT的快速算法； 1965年由Cooley和Tukey提出； 1969年，用一组2048点的傅里叶变换分析地震信号时，DFT需要13.5小时；而FFT仅用了2.4秒。直接计算直接计算DFTDFT的运算量的运算量复数序列的一般情况，对某一个k值，直接按上式计算X(k)值需要N次复数乘法、(N-1)次复数加法。减少运算量的途径： 将N点DFT分解为几个较短的DFT； 利用 的周期性和对称性；

60、10( )( ),0,1,1NknNnX kx n WkNnkNjnkNeW2- 周期性对称性 的周期性和对称性的周期性和对称性22()jm lNjmm lNmNNNNWeeW2mN mN mmNNNNNmmNNWWWWWW nkNW时间抽取法和频率抽取法时间抽取法和频率抽取法时间抽取法的基本原理时间抽取法的基本原理设序列x(n)的长度为N，且满足2 ,MNM为自然数 按n的奇偶把x(n)分解为两个N/2点的子序列12( )(2 ),0,1,12( )(21),0,1,12NxrxrrNxrxrr 则x(n)的DFT为由于所以 /2 1/2 11/22/21200( )( )( )( )( )