FFT算法介绍.ppt
本章内容：介绍傅里叶变换的一些快速算法快速算法的思想 根据原是变换定义的运算规律，及其中某些算子的特殊性，找出减少乘法和加法运算次数的有效途径，实现原始变换的各种高效算法。本讲学习目标本讲学习目标了解直接计算N点DFT的运算量了解减少运算量的基本途径理解按时间抽取的基-2FFT算法的算法原理、运算流图、所需计算量和算法特点了解按时间抽取的基-2FFT算法的编程思想本章作业练习本章作业练习 P127: 4.14.24.44.5第四章第四章 快速傅里叶变换快速傅里叶变换FFT: Fast Fourier Transform1965年，Cooley, Tukey机器计算傅里叶级数的一种算法4.2 基-2FFT算法4.2.1直接计算DFT的特点及减少运算量的基本途径10: ( ) ( )( )( )NnkNNnDFTX kDFT x nx n WRk10:1 ( )( )( )( )NnkNNkIDFTx nIDFT X kX k WRnN( )Nx n点有限长序列1. DFT的运算量及特点的运算量及特点运算量运算量复数乘法复数加法一个X(k)NN 1N个X(k)(N点DFT)N 2N (N 1)实数乘法实数加法一次复乘42一次复加2一个X (k)4N2N+2 (N 1)=2 (2N 1)N个X (k)(N点DFT)4N 22N (2N 1)10( )NnkNnx n Wajbcjdacbdj adcbnkNW 的特性*()() ()nknkN n kn N kNNNNWWWW对称性()() nkN n kn N kNNNWWW周期性 nkmnkNmNWW可约性/nknk mNN mWW0/2(/2) 11Nk NkNNNNWWWW 特殊点：2jnknkNNWeNknkNNWWnNnkNNWW2jmnkmNe221NjjNee FFT算法分类:时间抽选法DIT: Decimation-In-Time频率抽选法DIF: Decimation-In-FrequencyDFFFTDFTDFTDFTT 利用系数的特性，合并运算中的某些项， 把长序列短序列，从而算法的基本思想：减少其运算量。4.2.2、时间抽取法基、时间抽取法基-2FFT算法基本思想算法基本思想 （基（基-2 Decimation-In-Time FFT）1、算法原理设序列点数 N = 2M，M 为整数。 若不满足，则补零 12221xrx rxrxr0,1,.,/2 1rN将序列x(n)按n的奇偶分成两组：N为2的整数幂的FFT算法称基-2FFT算法。n为偶数时: n为奇数时:则x(n)的DFT: 111000NNNnknknkNNNnnnX kx n Wx n Wx n Wn为偶数n为奇数/2 1/2 121200221NNrkrkNNrrxr WxrW /2 1/2 1221200NNrkrkkNNNrrx rWWxrW /2 1/2 11/22/200NNrkkrkNNNrrx r WWxr W 12kNXkW Xk,0,1,./2 1r kN 其中,2.2.两点结论: (1) (1) X X (k),X(k),X (k)(k)均为N/2N/2点的DFTDFT。 (2) X(k)=X(2) X(k)=X (k)+W(k)+W X X (k)(k)只能确定出 X(k)X(k)的k= k= 个；即前一半的结果。1 21 2kN10102210101122222222) 12()()()2()()(NNNNNNNNrrkrrkrrkrrkWrxWrxkXWrxWrxkX1, 1 , 02 N 同理, 这就是说,X1(k),X2(k)的后一半,分别 等于其前一半的值。3.X(k)3.X(k)的后一半的确定的后一半的确定rkkrNNNWW222)()()()()2(1101)(101122222kXWrxWrxkNXNNNNNrrkkrr由于 （周期性）（周期性）, ,所以:)()2(22kXkNX 可可见,X(k)X(k)的后一半，也完全由X X1 1(k)(k), X X2 2 (k)(k)的前一半所确定。 * N点的DFT可由两个N/2点的DFT来计算。kNkNNkNWWWWNN22)()2()2()2(212NkXWNkXNkXNkN1, 1 , 0 ),()(221NkNkkXWkX又由于 , ,所以实现上式运算的流图称作蝶形运算 前一半4.4.蝶形运算蝶形运算 后一半（N/2个蝶形)(前一半)(后一半)1 1 11-1-1)()()()()()(2121kXWkXkXkXWkXkXkNkN)1,()1, 1 ,0(22 N Nk kk kN NN N)()()(21kXWkXkXkN)()()2(21kXWkXkNXkN)(1kX)(2kXkNW由X1(k)、X 2(k)表示X(k)的运算是一种特殊的运算-碟形运算图4.2.1 蝶形运算符号 CABA BCA BC x1 1(0)=(0)=x(0) (0) x1(1)=(1)=x(2) (2) N/2N/2点点 x1(2)=(2)=x(4) DFT (4) DFT x1(3)=(3)=x(6) (6) x2(0)=(0)=x(1) (1) x2(1)=(1)=x(3) (3) N/2N/2点点 x2(2)=(2)=x(5) (5) DFT DFT x2(3)=(3)=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)WN2WN1WN0WN3-1-1-1-1X(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)5.5.对X X1 1 (k)(k)和X X2 2 (k)(k)进行蝶形运算,前半部为 X(0)- X(3),X(0)- X(3),后半部分为X(4)- X(7)X(4)- X(7) 整个过程如下图所示:6. N/2分解后的运算量：分解后的运算量：复数乘法复数加法一个N / 2点DFT(N / 2)2N / 2 (N / 2 1)两个N / 2点DFTN 2 / 2N (N / 2 1)一个蝶形12N / 2个蝶形N / 2N总计22/2/2/2NNN2/2 1/2N NNN运算量减少了近一半 由于N=2N=2 M , ,所以 N/2N/2仍为偶数,可以进 一步把每个N/2N/2点的序列再按其奇偶部分 分解为两个N/4N/4的子序列。例如，n为偶数时的 N/2N/2点分解为:( (二二) N/4) N/4点点DFTDFT1, 1 , 0),()2(431Nlxlx1, 1 , 0),() 12(441Nlxlx进行N/4N/4点的DFTDFT,得到klNllkNllkNllkNlWlxWlxkXWlxWlxkXNNNN) 12(2/1014/104422/1014/1033) 12()()()2()()(4444( (偶中偶) )( (偶中奇) )13/2413/24( )( )( )()( )( )4kNkNX kXkWXkNX kXkWXk0,1,.,14Nk 25/2625/26( )( )( )()( )( )4kNkNXkXkWXkNXkXkWXk0,1,.,14Nk 同理可将奇数序号组成的N/2点序列进行分解得：其中：552( )( )(2 )XkDFT x lDFT xl662( )( )(21)XkDFT x lDFT xl0,1,.,/4 1lN2/2kkNNWW统一系数：X2的偶数序列X2的偶数序列这样逐级分解，直到2点DFT当N = 8时，即分解到X3(k)，X4(k)，X5(k)，X6(k)，k = 0, 10003322301033223(0)(0)(1)(0)(4)(1)(0)(1)(0)(4)NNXxWW xxW xXxWW xxW x/4 1133/43/400( )( )( )NlklkNNllXkx l Wx l W0,1k 0004422401044224(0)(0)(1)(2)(6)(1)(0)(1)(2)(6)NNXxWW xxW xXxWW xxW x/4 1144/44/400( )( )( )NlklkNNllXkx l Wx l W0,1k 不再做DFT (0) (4) (2) (6) (1) (5) (3) (7)WN0WN0WN0W0N-1-1-1-1X (0)X (1)X (0)X (1)X (0)X (1)X (0)X (1)33445566WN0WN2WN0WN2-1-1-1-1X (0)X (1)X (2)X (3)X (0)X (1)X (2)X (3)11121222WWWWN0N1N2N3-1-1-1-1X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)xxxxxxxx因此,8点DFT的FFT的运算流图如下4.2.3 基2DIT-FFT算法与直接DFT运算量的比较当N = 2M时，共有M级蝶形，每级N / 2个蝶形，每个蝶形有1次复数乘法，2次复数加法。2(2)log22MNNCMN复数乘法：2(2)logACNMNN复数加法：222()2(2)()loglog2MMCDFTNNNCFFTNN与DFT 比较10221048576204.8lN=21og5120024NN如 4.2.4 DIF-FFT的运算规律 及编程思想1）原位计算 DIT-FFT的运算过程很有规律，共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，与其它蝶形运算无关。 这样,蝶形运算的两个输出值仍可放回蝶形运算的两个输入所在的存储器中,这种利用同一存储单元存储蝶形计算输入、输出的方法即为原位运算。每一级( (列) )有N/2N/2个蝶形运算,所以只需N N个存储单元，可以节省存储单元。运算规律运算规律 (0)=(0)=A A0 0(0)(0) A A1 1(0)(0) A A2 2(0) A(0) A3 3(0)(0)=X(0) =X(0) (4)=(4)=A A0 0(1)(1) A A1 1(1) A(1) A2 2(1) A(1) A3 3(1)(1)=X(1)=X(1) (2)=(2)=A A0 0(2)(2) A A3 3(2)(2)=X(2)=X(2) (6)=(6)=A A0 0(3)(3) A A3 3(3)(3)=X(3)=X(3) (1)=(1)=A A0 0(4)(4) A A1 1(4) A(4) A2 2(4) A(4) A3 3(4)(4)=X(4)=X(4) (5)=(5)=A A0 0(5)(5) A A3 3(5)(5)=X(5)=X(5) (3)=(3)=A A0 0(6)(6) A A3 3(6)(6)=X(6)=X(6) (7)=(7)=A A0 0(7)(7) A A1 1(7) A(7) A2 2(7) A(7) A3 3(7)(7)=X(7)=X(7)WWWWN0N0N0N0-1-1-1-1WWWWN0N2N0N2-1-1-1-1WWWWNNNN0123. 原位运算的图示输入数据、中间运算结果和最后输出均用同一存储器。xxxxxxxxL=1L=2L=3开始时，输入序列的N个数放于此N个存储器内，倒序重排后仍存于这N个存储器中，每一次迭代运算后的结果也仍然存于这N个存储器中，整个运算完成后，这N个存储器中即为所求的频谱序列X(k) （k = 0、1、.、 N-1）。这就是所谓的同址计算，这样可以大大节约存储元件。 N点DITFFT运算流图中，每级都有N/2个蝶形。每个蝶形都要乘以因子WpN，称其为旋转因子，p称为旋转因子的指数。 1210 3210310 20 1122 22L-JPNJ JNPNJJNPNJJNPN,., JWW, ,J WW WL,JWWWLJWWWLLLLL一般情况下pNW 的确定2）旋转引子的变化规律LMJNJNJMLMLMJPWWWNLMMLL2W2222 :222PNL所以：又因为3）蝶形运算规律 对N = 2M点FFT，输入倒位序，输出自然顺序，设第L级运算每个蝶形的两节点距离为 B行，则第L级运算：1111-1( )( )()()( )():201.21 1,2,. PLLLNPLLLNM LXkXJXJB WXJBXJXJB WPJJLM式中：，的意义LMJNPNWW2LM 2旋转因子的确定仅与指数P有关，当L一定时J可以确定，进而可以确定指数P。对于同一P，参与本级蝶形运算的次数为 ，每级第一次调用 的蝶形结第一个输入节点的位置为J，第二次调用 的蝶形结第一个输入节点的位置为J+2L，即以2L为步进，搜索下一蝶形运算第一个输入节点的位置。以此类推，直至做完 个蝶形运算。LM 2PNWPNW 同一级中，同一P的蝶形计算完成后，计算下一个P值对应的蝶形运算，直至完成本级所有蝶形运算。DIT-FFT原位计算步骤1.确定L；2.求出第L级的 个 ；3.对每一个 完成其所参与的 个蝶形运算，4. 重复步骤13完成M级蝶形运算。LM 2PNWPNW12L )2()20()0()0(2)W0)0(X 0J12002222 112,.0 2302230N2211XWXXXXJJJJP BJLNLML-L（可进行蝶形运算蝶形运算参数取出后，则可计算出：算的执行过程。以第二级为例，说明运)3()21 () 1 () 1 ()21 () 1 (X 132223222XWXXXWXJNN02230223X (4)(42)(4)4 (4)(42)(6)NNXWXJXXWX22232223X (5)(52)(5)5 (5)(52)(7)NNXWXJXXWX以2L为步进参与本级同一旋转因子对应的其他蝶形运算作为作为练习练习请写请写出出L=3时的时的运算运算过程过程。 4） 编程思想及程序框图编程思想及程序框图在一般情况下，进行FFT运算的序列至少都有几百点的长度，因此需要编制FFT算法程序以便能够利用计算机来快速进行计算。输入倒位序，输出顺序的DIT-FFT的编程思想， N必须等于2的正整数幂，FFT的计算程序可以分为两部分：一部分是倒序重排，另一部分是用三层嵌套的循环来完成M=log2N次迭代。 三层循环的功能是：最里的一层循环完成相同WNP 的蝶形运算，中间的一层循环完成因子WNP的变化，而最外的一层循环则是完成M次迭代过程。 开 始送入x ( n )， MN 2 M倒 序L 1 , M 0 , B 1P 2 M LJk J , N 1 , 2LpNpNWBkXkXBkXWBkXkXkX)()()()()()(输 出结 束B 2 L 1循环循环L确定输入数据所在的位确定输入数据所在的位置置,循环循环JL、P、B确定后确定后进行蝶形计算进行蝶形计算 5 5） 序列的倒序倒序 由DIT-FFT的规律可知,输出X X(k)(k)按正常顺序排列在存储单元,而输入是按顺序: 这种顺序称作倒位序,即二进制数倒位。在编程时需完成倒位序，才能执行原位计算。）；）；7 (),3 (),5 (),1 (6 (),2(),4(),0 (xxxxxxxxn =00n =10n =01n =11n =01n =1101010101 x xx xx xx xx xx xx xx x),(012nnnx(n2)x(000) 0 x(100) 4 x(010) 2 x(110) 6 x(001) 1x(101) 5 x(011) 3 x(111) 7 （偶）（奇） 倒位序由奇偶分组造成,以N=8N=8为例 说明如下: 0 00 0 0 0 0 00 0 0 0 0 0 0 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 2权重：0122 2 2X(I)X(J)最低最低位加位加1若最高位为0，二进制最高位加1，对应的十进数加N/2=4若最高位为1，最高位置0，同时J=J-N/2，若次高位为0次高位加1,J=J+N/4若次高位为1，次高位置0同时J=J-N/4，次次高位加1,J=J+N/8顺序I起始及终止序号为：16倒序J起始序号为：N/2=4 当IJ 时，A(I)和A(J)的内容调换 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)x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)形成倒序形成倒序J后，将原存储器放的输入序列重新按倒序排列后，将原存储器放的输入序列重新按倒序排列。X(I)X(J)I=J时不需要交换 图4.2.9 倒序程序框图 221NNLHJNLHI1 , N1I JTJAJXIAIXT)()()()(J KLHK KJJ2KKKJJNNY确定确定J的起始序号及的起始序号及I的终的终止序号止序号K=N/2最高位为最高位为0否？否？J=J+N/2最高位为1，将最高位置0，J=J-N/2；次高位加1，K=N/4，倒序重排的程序是一段经典程序，它以巧妙的构思、简单的语句用高级编程语言来完成了 倒 序 重 排 的 功 能 。 下 面 是 一 段 用FORTRAN语言编写的倒序重排程序。 N=2*M （表示N=2M ，M是输入的正整数） LH=N/2 （LH是一个整数变量） N1=N-2 （N1也是一个整数变量） J=1 （对变量J赋初值） 100 DO 7 I=1,N1(循环开始，到语句7结束； 循环变量I从1开始，到N1结束，步长为1) IF (I.GE.J) GOTO 5（如果IJ，就转移到语句5） （将输入序列中序号互为倒序 的两个数值交换位置） TIXIXJXJXT)( )()()(5 K=N2 6 IF(K.GE.J) GOTO 7 (如果KJ，就转移到语句7) J=J-kK=K/2 GOTO 6 （转移到语句6）7 J=J+K 附：DIT算法的其他形式流图输入倒位序输出自然序输入自然序输出倒位序输入输出均自然序相同几何形状 输入倒位序输出自然序 输入自然序输出倒位序 用Matlab方法计算信号的DFT时，是用函数fft(x,N) 和ifft(x.N)。对于这两个函数，如果N为2的正整数幂，则可以得到本章中介绍的基2 FFT快速算法；如果N既不是2的正整数幂，也不是质数，则函数将N分解成质数，得到较慢的混合基 FFT算法；如果 N 为质数，则fft函数采用原来的 DFT 算法。附：利用附：利用Matlab计算计算FFT
