第4章快速傅里叶变换FFT

2019/12/191第4章快速傅里叶变换FFT24.1引言4.2基2FFT算法4.3进一步减少运算量的措施3一、DFT的运算量估计1,,1,0,)()(10NkWnxkXNnnkN101,,1,0,)(1)(NknkNNnWkXNnx两者的差别仅在指数的符号和因子1/N。4.1引言41210)1()2()1()0()1(NNNNNWNxWxWxWxX一个X(k)的值的工作量，如X(1)：则计算完全部的X(k)，k=0,1,2,…,N-1，需要共需N次复数乘法N-1次复数加法N×N=N2次复数乘法N×(N-1)N2次复数加法nkNW需将每一个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)共需4N2次实乘2N2+2N2=4N2次实加6[例]N=210=1024，总共需实时信号处理对计算速度有十分苛刻的要求。如何能减少DFT运算量？？400万次乘法400万次加法72.对称性：)(*)(kNnNnkNnkN的特征如下：nkNW1.周期性：)()(NknNkNnNnkN3.可约性：kNkNjkNjkNWeeW224.正交性：rNmnrNmnWNNkkmnN,0,1110)(二、提高DFT运算效率的途径1428WW如85.特殊的:j420,1,1如果充分利用复函数的这些特征，就可以改善DFT的运算效率。为了能利用上述特性，需将x(n)或X(k)这些长序列按一定的规律分解成一些短序列，即重排，以减少运算次数。nkNW二、提高DFT运算效率的途径的特征如下：nkNW91965年，库利(cooley)和图基(Tukey)首先提出FFT算法。对于N点DFT，仅需(N/2)log2N次复数乘法运算。[例]N=210=1024时：需要的复乘次数(1024/2)log2210=512×10=5120次，直接计算需要的复乘次数为1048576次，5120/1048576=4.88‰,所以利用FFT算法，速度提高200多倍。1969年，用一组2048点的傅里叶变换分析地震信号时，DFT需要13.5小时；而FFT仅用了2.4秒。10一、基(Radix)基2FFT算法的思想就是将长序列划为短序列，短到什么程度？用“基（Radix）”来表示，即基是FFT中用到的最小DFT单元。Radix-2（基2），代表2点DFT。2点DFT实际上只是加减运算。4.2基2FFT算法112点DFT02021002)1()0()()0(WxWxWnxXn1,0,)()(102kWnxkXnnk1202102)1()0()()1(WxWxWnxXnn即：12画出流图：)0(x)1(x)0(X)1(X1)1()0()1()1()0()0()1()0()1()0(12020202xxXxxXxx算法分为两类：时域抽取法FFT（DIT-FFT,Decimation-In-TimeFFT)，它是库利-图基算法；频域抽取法FFT（DIF-FFT,Decimation-In-FrequencyFFT)。二、DIT-FFT算法原理14设序列x(n)的长度为：N=2M（0,1,2,……)，按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点分解161010)()(NnNnnkNnkNWnxWnx(n为偶数)(n为奇数)10)()]([)(NnnkNWnxnxDFTkX10)12(10222)12()2(NNrkrNrrkNWrxWrx1022102122)()(NNrrkNkNrrkNWrxWWrx)()(21kXWkXkN1,,2,1,0Nk17(1)X1(k)，X2(k)均为N/2点的DFT。(2)只能确定出X(k)的的值，即前一半的结果。1,,1,02Nk)()()(21kXWkXkXkN10102210101122222222)12()()()2()()(NNNNNNNNrrkrrkrrkrrkWrxWrxkXWrxWrxkX其中2.结论12,,2,1,0Nk)()()(21kXWkXkXkN1,,2,1,0Nk18同理，3.X(k)后一半的确定由于（周期性），所以：rkkrNNNWW222)()()()()2(1101)(101122222kXWrxWrxkNXNNNNNrrkkrr)()2(22kXkNXX1(k)、X2(k)后一半的值，分别等于其前一半的值。12,,2,1,0Nk19可见，X(k)的后一半，也完全由X1(k)、X2(k)的前一半所确定。)2()2()2(212NkXWNkXNkXNkN1,,1,0),()(221NkNkkXWkXkNkNNkN22)(又由于，所以：前一半后一半)()()()()()(2121kXWkXkXkXWkXkXkNkN)1,,()1,,1,0(22NkkNN20只要求出0～N/2-1区间内的各个k值所对应的X1(k)、X2(k)的值，即可以求出0～N-1整个区间内的全部X(k)值。这就是FFT能大量节省运算量的关键。N点的DFT可由两个N/2点的DFT来计算。前一半后一半)()()()()()(2121kXWkXkXkXWkXkXkNkN)1,,()1,,1,0(22NkkNN214.蝶形运算碟形运算符号ABCBACBAC22实现上式运算，共需N/2个蝶形运算。前一半后一半)()()()()()(2121kXWkXkXkXWkXkXkNkN)1,,()1,,1,0(22NkkNN由X1(k)、X2(k)表示X(k)的运算是一种碟形运算：)()()(21kXWkXkXkN)()()2(21kXWkXkNXkN)(1kX)(2kXkNW23(1)N/2点的DFT运算量:复乘次数，复加次数(2)两个N/2点的DFT运算量:复乘次数，复加次数(3)N/2个蝶形运算的运算量:复乘次数，复加次数5.运算量分析4)2(22NN)12(2NN22N)12(NN2NNN22按奇、偶分组后的计算量：24复乘:复加:2)12(2NNNN22222NNN但是，N点DFT的复乘为N2，复加N(N-1)；与分解后相比，计算工作量差不多减少一半。5.运算量分析一个N点DFT分解成两个N/2点DFT后的总共运算量：25例：N=8时的DFT，可以分解为两个N/2=4点的DFT。具体方法如下：(1)n为偶数时,即分别记作:3,2,1,0)2()()(30430411kWrxWrxkXrrkrrk，　)6()3(),4()2(),2()1(),0()0(1111xxxxxxxx)6(),4(),2(),0(xxxx进行N/2=4点的DFT，得：26(2)n为奇数时，分别记作：)7()3(),5()2(),3()1(),1()0(2222xxxxxxxx3,2,1,0)12()()(30430422kWrxWrxkXrrkrrk，　)()()4()()()(2121kXWkXkXkXWkXkXkNkN进行N/2=4点的DFT，得：3,2,1,0k27x1(0)=x(0)x1(1)=x(2)N/2点x1(2)=x(4)DFTx1(3)=x(6)x2(0)=x(1)x2(1)=x(3)N/2点x2(2)=x(5)DFTx2(3)=x(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)WN2WN1WN0WN3X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)(3)对X1(k)和X2(k)进行蝶形运算,前半部分为X(0)～X(3),后半部分为X(4)～X(7)。28由于N=2M,所以N/2仍为偶数,可以进一步把每个N/2点的序列再按其奇偶部分分解为两个N/4的子序列。例如，n为偶数时的N/2点分解为:(二)N/4点DFT1,,1,0),4()2()(413Nllxlxlx1,,1,0),24()12()(414Nllxlxlx进行N/4点的DFT，得：klNllkNllkNllkNlWlxWlxkXWlxWlxkXNNNN)12(2/1014/104422/1014/1033)12()()()2()()(4444(偶中偶)(偶中奇)29)()()(4312kXWkXkXkN1,,1,04Nk从而可得到前N/4的X1(k)：)()()4(4312kXWkXkNXkN后N/4的X1(k)：1,,1,04Nk(二)N/4点DFT即对X3(k)、X4(k)进行一次碟形运算。30(奇中偶)104/64/1026104/5104/254444)()12()()()2()(NNNNllkNlkNlllkNllkNWlxWlxkXWlxWlxkX(奇中奇)同样对n为奇数时，N/2点分为两个N/4点的序列进行DFT，则有：(k)XW(k)X(k)X6kN/252(k)XW(k)Xk)4N(X6kN/2521,,1,0),14()2()(425Nllxlxlx1,,1,0),34()12()(426Nllxlxlx14,,1,0Nk由X5(k)、X6(k)进行碟形运算，得：31所以，N=8时的DFT可分解为四个N/4的DFT,具体步骤如下:1,0),4()2()(13llxlxlx(1)原序列x(n)的“偶中偶”部分：构成N/4点DFT，从而得到X3(0)，X3(1)。0,1k,)()(104/33llkNWlxkX)4()0()1()0()1()4()0()1()0()0(031233030233xWxxWxXxWxxWxXNN)4()2()1()0()0()0(1313xxxxxx321,0),24()12()(14llxlxlx(2)原序列x(n)的“偶中奇”部分：构成N/4点DFT，从而得到X４(0)，X４(1)。0,1k,)()(104/44llkNWlxkX)6()2()1()0()1()6()2()1()0()0(041244040244xWxxWxXxWxxWxXNN)6()3()1()2()1()0(1414xxxxxx33(3)原序列x(n)的“奇中偶”部分：1,0),14()2()(25llxlxlx构成N/4点DFT，从而得到X５(0)，X５(1)。0,1k,)()(104/55llkNWlxkX)5()2()1()1()0()0(2525xxxxxx)5()1()1()0()1()5()1()1()0()0(051255050255xWxxWxXxWxxWxXNN34(4)原序列x(n)的“奇中奇”部分：1,0),34()12()(26llxlxlx构成N/4点DFT，从而得到X６(0)，X６(1)。0,1