第4章快速傅里葉變換(FFT)1

1、第第4章章 快速傅里葉變換快速傅里葉變換(FFT) n4.1 引言引言n4.2 基基2FFT算法算法n4.3 進一步減少運算量的措施進一步減少運算量的措施n4.4 其他快速算法其他快速算法4.1 引言引言n DFT是信號分析與處理中的一種是信號分析與處理中的一種重要變換。因直接計算重要變換。因直接計算DFT的計算量與的計算量與變換區(qū)間長度變換區(qū)間長度N的平方成正比，當?shù)钠椒匠烧?，當N較大較大時，計算量太大，所以在快速傅里葉變時，計算量太大，所以在快速傅里葉變換換(簡稱簡稱FFT)出現(xiàn)以前，直接用出現(xiàn)以前，直接用DFT算法算法進行譜分析和信號的實時處理是不切實進行譜分析和信號的實時處理是不切實

2、際的。直到際的。直到1965年發(fā)現(xiàn)了年發(fā)現(xiàn)了DFT的一種快的一種快速算法以后，情況才發(fā)生了根本的變化。速算法以后，情況才發(fā)生了根本的變化。 4.2 基基2FFT算法算法n 4.2.1 直接計算直接計算DFT的特點及減少運算量的的特點及減少運算量的基本途徑基本途徑 n 長度為長度為N的有限長序列的有限長序列x(n)的的DFT為為n 考慮考慮x(n)為復數(shù)序列的一般情況，對為復數(shù)序列的一般情況，對某一個某一個k值，直接計算值，直接計算X(k)值需要值需要N次復數(shù)乘法、次復數(shù)乘法、(N-1)次復數(shù)加法。次復數(shù)加法。 10( )( ),0,1,1NknNnX kx n WkNn 如前所述，如前所述，N

3、點點DFT的復乘次數(shù)的復乘次數(shù)等于等于N2。顯然，把。顯然，把N點點DFT分解為幾個較分解為幾個較短的短的DFT，可使乘法次數(shù)大大減少。另，可使乘法次數(shù)大大減少。另外，旋轉(zhuǎn)因子外，旋轉(zhuǎn)因子WmN具有明顯的周期性和具有明顯的周期性和對稱性。其周期性表現(xiàn)為對稱性。其周期性表現(xiàn)為22()jm lNjmm lNmNNNNWeeW其對稱性表現(xiàn)為2mN mN mmNNNNNmmNNWWWWWW 4.2.2 時域抽取法基時域抽取法基2FFT基本原理基本原理 FFT算法基本上分為兩大類：時域算法基本上分為兩大類：時域抽取法抽取法FFT(Decimation In Time FFT,簡簡稱稱 D I T - F

4、 F T ) 和 頻 域 抽 取 法和 頻 域 抽 取 法FFT(Decimation In Frequency FFT,簡簡稱稱DIFFFT)。下面先介紹。下面先介紹DIFFFT算算法。法。設序列設序列x(n)的長度為的長度為N，且滿足，且滿足2 ,MNM為自然數(shù) 按按n的奇偶把的奇偶把x(n)分解為兩個分解為兩個N/2點的子序列點的子序列12( )(2 ),0,1,12( )(21),0,1,12Nx rxrrNx rxrrn 則x(n)的DFT為/2 1/2 12(21)00/2 1/2 121200( )( )( )(2 )(21)( )( )knknNNnnNNkrkrNNrrNNk

5、krNNrrX kx n Wx n Wxr WxrWx rWx r W由于222222/2jkrNjkrkrkrNNNWeeW所以 /2 1/2 11/22/21200( )( )( )( )( )NNkrkkrkNNNNrrX kx r WWx r WX kW Xkn 其中X1(k)和X2(k)分別為x1(r)和x2(r)的N/2點DFT，即 /2 111/210/2 122/220( )( )( )( )( )( )NkrNrNkrNrX kx r WDFT x rXkx r WDFT x r 由于X1(k)和X2(k)均以N/2為周期，且 ，所以X(k)又可表示為2NkkNNWW 121

6、2( )( )( )0,1,12()( )( )0,1,122kNkNNX kX kW XkkNNX kX kW Xkk圖4.2.1 蝶形運算符號 CABA BCA BC圖4.2.2 N點DFT的一次時域抽取分解圖(N=8) 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)X(7)n 與第一次分解相同，將x1(r)按奇偶分解成兩個N/4長的子序列x3(l)和x4(l)，即3241( )(2 )

7、,0,1,1( )(21)4x lxlNlx lxl那么，X1(k)又可表示為 /4 1/4 12(21)11/21/200/4 1/4 13/4/24/4003/24( )(2 )(21)( )( )( )( ),0,1,/21NNklklNNiiNNklkklNNNiikNX kxl WxlWx l WWx l Wx kWXk kN式中 /4 133/430/4 144/440( )( )( )( )( )( )NklNiNklNix kx l WDFT x lxkxl WDFT xl 同理，由X3(k)和X4(k)的周期性和Wm N/2的對稱性 Wk+N/4 N/2=-Wk N/2 最后

8、得到： 13/2413/24( )( )( ),0,1,/41(/4)( )( )kNkNX kXkWXkkNX kNXkWXkn 用同樣的方法可計算出25/2625/26( )( )( ),0,1,/41(/4)( )kNkNXkXkWXkkNXkNX kWXk其中 /4 155/450/4 166/4605262( )( )( )( )( )( )( )(2 ),0,1,/41( )(21)NklNiNklNiXkx l WDFT x lXkx l WDFT x lx lxllNx lxl圖4.2.3 N點DFT的第二次時域抽取分解圖(N=8) N/4點DFTWN12WN12WN0WN1W

9、N2WN3X1(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 圖4.2.4 N點DITFFT運算流圖(N=8) WN0WN1WN2WN3WN0WN2WN0WN2WN0WN0WN0WN0 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

10、)A(4)A(5)A(6)A(7)A(0)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)n 4.2.3 DITFFT算法與直接計算DFT運算量的比較n 每一級運算都需要N/2次復數(shù)乘和N次復數(shù)加(每個蝶形需要兩次復數(shù)加法)。所以，M級運算總共需要的復數(shù)乘次數(shù)為22(2)log22(2)logMANNCMNCN MNN復數(shù)加次數(shù)為 例如，N=210=1024時221048576204.8(/2)log5120NNN圖4.2.5 FFT算法與直接計算DFT所需乘法次數(shù)的比較曲線 由于N=2L，因此N/2仍為偶

11、數(shù)，可以依照上面方法進一步把每個N/2點子序列，再按輸入n的奇偶分解為兩個N/4點的子序列，按這種方法不斷劃分下去，直到最后剩下的是2點DFT，兩點DFT實際上只是加減運算加減運算。蝶形結(jié)即蝶式計算結(jié)構(gòu)也即為蝶式信號流圖上面頻域頻域中前/后半部分表示式可以用蝶形信號流圖表示。X1(k)X2(k)kNW)()(21kXWkXkN)()(21kXWkXkN作圖要素：作圖要素：(1)左邊兩路為輸入左邊兩路為輸入(2)右邊兩路為輸出右邊兩路為輸出(3)中間以一個小圓表示加、中間以一個小圓表示加、減運算（右上路為相加輸減運算（右上路為相加輸出、右下路為相減輸出出、右下路為相減輸出)(4)如果在某一支路上

12、信號需要進行如果在某一支路上信號需要進行相乘運算，則在該支路上標以箭頭，相乘運算，則在該支路上標以箭頭，將相乘的系數(shù)標在箭頭旁。將相乘的系數(shù)標在箭頭旁。(5)當支路上沒有箭頭及系當支路上沒有箭頭及系數(shù)時，則該支路的傳輸比數(shù)時，則該支路的傳輸比為為1。例子：求 N=23=8點FFT變換 （1）先按）先按N=8-N/2=4，做做4點的點的DFT：)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN12/, 0Nk先將N=8 DFT分解成2個4點DFT：可知：時域上：x(0),x(2),x(4),x(6)為偶子序列 x(1),x(3),x(5),x(7)為奇子序列 頻域上：X

13、(0)X(3)，由X(k)給出 X(4)X(7),由X(k+N/2)給出(a)比較N=8點直接DFT與分解2個4點DFT的FFT運算量N=8點的直接DFT的計算量為：N2次(64次）復數(shù)相乘,N(N-1)次(8(8-1)=56次)復數(shù)相加.共計120次。（b)求 一個蝶形結(jié)需要的運算量要運算一個蝶形結(jié)，需要一次乘法 ， 兩次加法。)(2kXWkN)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN（c）分解為兩個N/2=4點的DFT的運算量分解2個N/2點（4點）的DFT：) 12()2()rxDFTrxDFTkX（偶數(shù)其復數(shù)相乘為復數(shù)相加為奇數(shù)其復數(shù)相乘為復數(shù)相加為2

14、2）（N）（122NN22）（N）（122NN次。共計為次，（次，復數(shù)相加為復數(shù)相乘為5624) 123222NNN(d)用2個4點來求N=8點的FFT所需的運算量再將N/2點（4點）合成N點(8點)DFT時，需要進行N/2個蝶形運算)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN12/, 0Nk還需N/2次（4次）乘法 及 次(8次）加法運算。22N)(2kXWkN計算量。分解就可以節(jié)省約一半次。看出僅做一次共計為次復數(shù)加法次復數(shù)乘法）共需求點所以68,322) 12(,36221(22)(8222NNNNNNNNNkXFFT(e)將N=8點分解成2個4點的DFT

15、的信號流圖 4點DFTx(0)x(2)x(4)x(6)4點DFTx(1)x(3)x(5)x(7)08W18W28W38WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)X1(0)X1(1)X1(2)X1(3)X2(0)X2(1)X2(2)X2(3)偶數(shù)序列奇數(shù)序列3821282118210821) 3 () 3 (3) 2() 2(2) 1 () 1 (1) 0() 0(0WXXXWXXXWXXXWXXX）（）（）（）（如：X(4)X(7)同學們自已寫x1(r)x2(r)(2)N/2(4點)-N/4(2點)FFT(a)先將先將4點分解成點分解成2點的點的DFT：n因為4點DFT

16、還是比較麻煩，所以再繼續(xù)分解。n若將N/2(4點)子序列按奇/偶分解成兩個N/4點（2點）子序列。即對將x1(r)和x2(r)分解成奇、偶兩個N/4點（2點）點的子序列。奇序列、偶序列、) 6() 2() 4() 0(: )(1xxxxrx奇序列、偶序列、同理：)7() 3 () 5 () 1 (: )(2xxxxrx1014012()()2(4131，），在此（）奇序列（）偶序列若設：LNLLXLxLxLx1014012()()2(5252，），在此（）奇序列（）偶序列同理：LNLLXLxLxLx(b)求2點的DFT)()4/()()()()()()()4/10)()() 12()2()()

17、(62/5262/52242/3142/314/0) 12(2/114/022/1111LXWLXNkxLXWLXkxrxLXWLXNkXLLXWLXWLXWLXkXkXrxkNkNkNkNNLkLNNLLkNDFT）也可分解為：同樣，（同理：，其中（可分解為：（c)一個2點的DFT蝶形流圖2點DFT2點DFTx(0)x(4)x(2)x(6)X3(0)X3(1)X4(0)X4(1)X1(0)X1(1)X1(2)X1(3)04W14W) 1 () 1 () 3 () 0() 0() 2() 1 () 1 () 1 () 0() 0() 0(41431404314143140431XWXXXWXX

18、XWXXXWXX其中（d)另一個2點的DFT蝶形流圖2點DFT2點DFTx(1)x(5)x(3)x(7)X5(0)X5(1)X6(0)X6(1)X2(0)X2(1)X2(2)X2(3)04W14W) 1 () 1 () 3 () 0() 0() 2() 1 () 1 () 1 () 0() 0() 0(61452604526145260452XWXXXWXXXWXXXWXX其中同理：(3)將N/4（2點）DFT再分解成2個1點的DFT(a)求2個一點的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0

19、()0()()(2WWknWWWWWxWxWxWxXWxWxWxWxXWnxkXnknknkNnkNnnkN，其中，則這里用到對稱性這是一蝶形結(jié)代入上面流圖可知： 最后剩下兩點DFT,它可分解成兩個一點DFT，但一點DFT就等于輸入信號本身，所以兩點DFT可用一個蝶形結(jié)表示。取x(0)、x(4)為例。(b)2個1點的DFT蝶形流圖 1點DFTx(0)1點DFTx(4)X3(0)X3(1)02W進一步簡化為蝶形流圖：02WX3(0)X3(1)x(0)x(4)4()0()4()0() 1 ()4()0()4()0()0(023023xxxWxXxxxWxX其中：(4)一個完整N=8的按時間抽取FF

20、T的運算流圖 x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)38W28W18W08W08W08W08W08W08W08W28W28WX(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)12/0NNNWW其中旋轉(zhuǎn)因子，共有m=0m=1m=2n 4.2.4 DITFFT的運算規(guī)律及編程思想n 1.原位計算n 由圖4.2.4可以看出，DITFFT的運算過程很有規(guī)律。N=2M點的FFT共進行M級運算，每級由N/2個蝶形運算組成。同一級中，每個蝶形的兩個輸入數(shù)據(jù)只對本蝶形有用，且每個蝶形的輸入、輸出數(shù)據(jù)結(jié)點又在同一水平線上，計算完一個蝶形后，所得輸出數(shù)據(jù)可立即存入原數(shù)據(jù)占用

21、的存儲單元，這就是位計算。原位計算可大節(jié)省內(nèi)存單元，降低成本。n 原位運算（in-place)n原位運算的結(jié)構(gòu)，可以節(jié)省存儲單元，降低設備成本。n定義：當數(shù)據(jù)輸入到存儲器以后，每一組運算的結(jié)果，仍然存放在這同一組存儲器中直到最后輸出。02Wx(0)x(4)X3(0)X3(1)中放在一個暫存器中放在單元中，將放在單元例：將RWAxAx08) 1 ()4()0()0(單元送回單元送回將) 1 ()4()0()0()4()0(0808AxWxAxWxn例：N=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)

22、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(3)A(4)=x(4)A(5)=x(5)A(6)=x(6)A(7)=x(7)382818084,3,2,1WRWRWRWR暫存器R1R1R1R1R1R3R1R1R3R2R3R4看出：用原位運算結(jié)構(gòu)運算后，A(0)A(7)正好順序存放X(0)X(7),可以直接順序輸出。2.旋轉(zhuǎn)因子的變化規(guī)律 如上所述，N點DITFFT運算流圖中，每級都有N/2個蝶形。每個蝶形都要乘以因子WpN，稱其為旋轉(zhuǎn)因子，

23、p稱為旋轉(zhuǎn)因子的指數(shù)。 觀察圖4.2.4不難發(fā)現(xiàn)，第L級共有2L-1個不同的旋轉(zhuǎn)因子。N=23=8時的各級旋轉(zhuǎn)因子表示如下： L=1時，WpN=WJ N/4=WJ2L， J=0 L=2時， WpN =WJ N/2=WJ2L， J=0，1 L=3時， WpN =WJN=WJ2L， J=0，1，2，3 382818082814080402080,;2, 1 ,0,2,; 1 ,0,2;0, 1,8:WWWWJLWWWWJLWWWJLNN如LMLJNJNpNMLMLMLLppNJPJWWWNJWWLMMLL212, 2 , 1 , 0,222212, 2 , 1 , 0,12212對N=2M的一般情

24、況，第L級的旋轉(zhuǎn)因子為 3. 蝶形運算規(guī)律 設序列x(n)經(jīng)時域抽選(倒序)后，存入數(shù)組X中。如果蝶形運算的兩個輸入數(shù)據(jù)相距B個點，應用原位計算，則蝶形運算可表示成如下形式： XL (J) XL-1(J)+X L-1(J+B)WpN XL(J+B) XL-1(J)-X L-1(J+B)WpN 式中 p=J2 M-L；J=0,1,，2 L-1-1；L=1,2,，M 下標L表示第L級運算，XL(J)則表示第L級運算后數(shù)組元素X(J)的值。如果要用實數(shù)運算完成上述蝶形運算，可按下面的算法進行。 設: T=X L-1(J+B)WpN=TR+jTI X L-1(J)=XR(J)+jXI(J) 式中下標R

25、表示取實部，I表示取虛部，pNBJXpNBJXTpNBJXpNBJXTpNjpNBJjXBJXWBJXTRIIIRRIRpNL2sin)(2cos)(2sin)(2cos)(2sin2)cos()()(1IIIRRRIIRRIRIRpNLLIRLIIIRRRIIRRIRIRpNLLIRLTJXBJXTJXBJXTJXjTJXjTTJjXJXWBJXJXBJjXBJXBJXTJXJXTJXJXTJXjTJXjTTJjXJXWBJXJXJjXJXJX)()()()()()()()()()()()()()()()()()()()()()()()()()(1111則WN0WN1WN2WN3WN0WN2

26、WN0WN2WN0WN0WN0WN0 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(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)n 4. 編程思想及程序框圖圖4.2.6 DITFFT運算和程序框圖 開 始送入x(n)， MN2 M倒 序L1 , M0 , B 1P2 M LJk J , N1 , 2LpNpNWBkXkXBkXWBkXkXkX)()()

27、()()()(輸 出結(jié) 束B 2 L1第第L級中，每個蝶形的兩個輸級中，每個蝶形的兩個輸入數(shù)據(jù)相距入數(shù)據(jù)相距B=2L-1個點；同一個點；同一旋轉(zhuǎn)因子對應著間隔為旋轉(zhuǎn)因子對應著間隔為2L點點的的2M-L個蝶形。個蝶形。L=1,B=1L=2,B=2L=3,B=4n 5. 序列的倒序n DITFFT算法的輸入序列的排序看起來似乎很亂，但仔細分析就會發(fā)現(xiàn)這種倒序是很有規(guī)律的。由于N=2M，所以順序數(shù)可用M位二進制數(shù)(nM-1nM-2n1n0)表示。 圖4.2.7 形成倒序的樹狀圖(N=23) 01010101010101(n2n1n0)200004261537100010110001101011111

28、表4.2.1 順序和倒序二進制數(shù)對照表 圖4.2.8 倒序規(guī)律 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)倒序倒序子程序具體執(zhí)行時，只須將1與4對調(diào)送入，3與6對調(diào)送入，而0，2，5，7不變，僅需要一個中間存儲單元。I01234567J04261537 在實際運算時，先按自然順序?qū)⑤斎胄蛄写嫒氪鎯卧偻ㄟ^變址運算將自然順序變換成按時間抽取的FFT算法要求的順序。變址的過程可以用程

29、序安排加以實現(xiàn)-稱為“整序”或“重排”（采用碼位倒讀）且注意：(1)當I=J時，數(shù)據(jù)不必調(diào)換；(2)當IJ時，必須將原來存放數(shù)據(jù)x(n)送入暫存器R，再將x(J)送入x(I),R中內(nèi)容送x(J).進行數(shù)據(jù)對調(diào)。(3)為了避免再次考慮前面已調(diào)換過的數(shù)據(jù)，保證調(diào)換只進行一次，否則又變回原狀。IN/2=4，做做4點的點的DIF：nNWNnxnxnxNnxnxnx)2()()()2()()(21先將N=8DFT分解成2個4點DFT：可知：時域上：x(0),x(1),x(2),x(3)為偶子序列 x(4),x(5),x(6),x(7)為奇子序列 頻域上：X(0),X(2),X(4),X(6)由x1(n)

30、給出 X(1),X(3),X(5),X(7),由x2(n)給出將N=8點分解成2個4點的DIF的信號流圖 4點DFTx(0)x(1)x(2)x(3)4點DFTx(4)x(5)x(6)x(7)X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)X1(k)前半部分序列后半部分序列nNnNnNnNWxxxWxxxWxxxWxxxxxxxxxxxxxxx)7()3(3,)6()2(2,)5() 1 (1,)4()0(0)7()3(3),6()2(2),5() 1 (1),4()0(022221111）（）（）（）（）（）（）（）（如：08W18W28W38Wx1(n)x2(n)X2(k)(

31、2)N/2(4點)-N/4(2點)FFT(a)先將先將4點分解成點分解成2點的點的DIF：n因為4點DIF還是比較麻煩，所以再繼續(xù)分解。后半部分序列、前半部分序列、) 3 () 2() 1 () 0(: )(11111xxxxnx10140)2()()()2()(411311，），在此（）（若設：LNLLxWNnxnxLxNnxnxnN后半部分序列、前半部分序列、同理：) 3 () 2() 1 () 0(: )(22222xxxxnx10140)2()()()2()(622522，），在此（）（同理：LNLLxWNnxnxLxNnxnxnN（b)一個2點的DIF蝶形流圖2點DFT2點DFTx1

32、(0)x1(1)x1(2)x1(3)x3(0)x3(1)x4(0)x4(1)X(0)X(4)X(2)X(6)08W28W) 3 () 1 () 1 (,) 2() 0() 0(113113xxxxxx其中(c)另一個2點的DIF蝶形流圖2點DFT2點DFTx2(0)x2(1)x2(2)x2(3)x5(0)x5(1)x6(0)x6(1)X(1)X(5)X(3)X(7)08W28W(3)將N/4（2點）DFT再分解成2個1點的DFT(a)求2個一點的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0()0

33、()()(2WWknWWWWWxWxWxWxXWxWxWxWxXWnxkXnknknkNnkNnnkN，其中，則這里用到對稱性這是一蝶形結(jié)代入上面流圖可知： 最后剩下兩點DFT,它可分解成兩個一點DFT，但一點DFT就等于輸入信號本身，所以兩點DFT可用一個蝶形結(jié)表示。取x3(0)、x3(1)為例。(b)2個1點的DFT蝶形流圖 1點DFTx3(0)1點DFTx3(1)X(0)X(4)02W進一步簡化為蝶形流圖：02WX(0)X(4)x3(0)x3(1)(4)一個完整N=8的按頻率抽取FFT的運算流圖 x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)38W28W18W08W08

34、W08W08W08W08W08W28W28WX(0)X(4)X(2)X(6)X(1)X(5)X(3)X(7)12/0NNNWW其中旋轉(zhuǎn)因子，共有m=0m=1m=2圖4.2.14 DITFFT的一種變形運算流圖WN0WN0WN2WN0X(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)WN0WN2WN1WN3WN2WN0WN0WN0按時間抽取的FFT算法的若干變體 對于任何流圖，只要保持各節(jié)點所連續(xù)的支路及其傳輸系數(shù)不變，則不論節(jié)點位置怎么排列所得流圖總是等效的，最后所得結(jié)果都是x(n)的離散付里葉變換的正確結(jié)果。只是數(shù)

35、據(jù)的提取和存放的次序不同而已。圖4.2.15 DITFFT的一種變形運算流圖WN0WN0WN2WN0X(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)WN0WN2WN1WN3WN2WN0WN0WN0DIF的特點(a)輸入自然順序，輸出亂序且滿足碼位倒置規(guī)輸入自然順序，輸出亂序且滿足碼位倒置規(guī)則。則。(b)根據(jù)時域、頻域互換，可知：根據(jù)時域、頻域互換，可知：輸入亂序，輸出自然順序。輸入亂序，輸出自然順序。DIF與DIT比較n相同之處：n（1）DIF與DIT兩種算法均為原位運算。n（2）DIF與DIT運算量相同。n它們都

36、需要算法是兩種等價的與次復加次復乘FFTDITDIFNaNmNFNF22lglg2n不同之處：n（1）DIF與DIT兩種算法結(jié)構(gòu)倒過來。nDIF為輸入順序，輸出亂序。運算完畢再運行“二進制倒讀”程序。nDIT為輸入亂序，輸出順序。先運行“二進制倒讀”程序，再進行求DFT。n（2）DIF與DIT根本區(qū)別：在于蝶形結(jié)不同。nDIT的復數(shù)相乘出現(xiàn)在減法之前。nDIF的復數(shù)相乘出現(xiàn)在減法之后。n 4.2.6 IDFT的高效算法 n 上述FFT算法流圖也可以用于離散傅里葉逆變換(Inverse Discrete Fourier Transform，簡稱IDFT)。比較DFT和IDFT的運算公式： IDF

37、TFFTNWWDFTWnxnxDFTkXWkXNkXIDFTnxnkNnkNNknkNNknkN算法都可以拿來運算或頻率抽取抽取）那么以上討論的時間（將運算結(jié)果都除以改成運算中的每個系數(shù)只要把3)2() 1 ()()()()(1)()(1010圖4.2.16 DITIFFT運算流圖 WN0WN1WN2WN3WN0WN0N1x(0)x(4)x(2)x(6)x(4)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)WN2WN2N1N1N1N1N1N1N1圖4.2.17 DITIFFT運算流圖(防止溢出) WN02121x(0)x(4)x(2)x(6)x(1)x

38、(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)212121WN121WN221WN3212121WN021WN2212121WN021WN22121WN02121WN0212121WN021WN021n 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方法：n 由于 10101( )( )1( )( )NknNkNknNkx nX k WNx nXk WN對上式兩邊同時取共軛，得1011( )( )( )NknNkx nXk WDFT XkNN可知：只須將頻域成份一個求共軛變換，即（1）將X(k)的虛部乘以-1，即先取X(k)的共軛，得X*(k)。

39、(2)將X*(k)直接送入FFT程序即可得出Nx*(n)。(3)最后再對運算結(jié)果取一次共軛變換，并乘以常數(shù)1/N，即可以求出IFFT變換的x(n)的值。直接利用FFT流圖方法的注意點n（1）FFT與IFFT連接應用時，注意輸入輸出序列的排列順序，即應注意是自然順序還是倒序。n（2）FFT和IFFT共用同一個程序時，也應注意利用FFT算法輸入輸出的排列順序，即應注意自然順序還是例位序n（3）當把頻率抽取FFT流圖用于IDFT時，應改稱時間抽取IFFT流圖。n（4）當把時間抽取FFT流圖用于IDFT時，應改稱頻率抽取IFFT流圖。4.3 進一步減少運算量的措施n 4.3.1 多類蝶形單元運算n 由

40、DITFFT運算流圖已得出結(jié)論，N=2M點FFT共需要MN/2次復數(shù)乘法。n 由(4.2.12)式，當L=1時，只有一種旋轉(zhuǎn)因子W0N=1，所以，第一級不需要乘法運算。 n 綜上所述，先除去第一、二兩級后，所需復數(shù)乘法次數(shù)應是n 從L=3至L=M共減少復數(shù)乘法次數(shù)為(2)(2)2MNCM13312( )2222MMLLLLNNN因此，DITFFT的復乘次數(shù)降至 (2)(2)(2)(3)2222MNNNCMM22(1)()()222()()22()222()()22defjxjyxjyjxyxyj xyRjIRxyIxyyx n 從實數(shù)運算考慮，計算N=2M點DITFFT所需實數(shù)乘法次數(shù)為(2)

41、4(3)2(2)2213(2)102MNNRMNMn 4.3.2 旋轉(zhuǎn)因子的生成n 在FFT運算中，旋轉(zhuǎn)因子WmN=cos(2m/N)-jsin(2m/N)，求正弦和余弦函數(shù)值的計算量是很大的。 n 4.3.3 實序列的FFT算法 n 設x(n)為N點實序列，取x(n)的偶數(shù)點和奇數(shù)點分別作為新構(gòu)造序列y(n)的實部和虛部，即1212( )(2 ),( )(21),0,1,12( )( )( ),0,1,12Nx nxnx nxnnNy nx njx n n對y(n)進行N/2點FFT，輸出Y(k)，則1122( )( )( ),0,1,1( )( )( )2epopX kDFT x nYkN

42、kXkDFT x njYk 根據(jù)DITFFT的思想及式(4.2.7)和(4.2.8)，可得到 12( )( )( ),0,1,2kNNX kX kW Xk kn實信號的快速循環(huán)卷積實信號的快速循環(huán)卷積n 用DITFFT計算兩個實信號x1(n)和x2(n)的循環(huán)卷積時，最直接的方法是將信號的虛部都置為零，再按復序列用FFT計算。顯然這樣做是很浪費的?，F(xiàn)在我們可用基2DITFHT直接進行實正交變換來處理實信號的循環(huán)卷積問題。1122( )( )( )( )FHTHFHTHx nXkx nXk計算式 )()()()(21nxnxnykYIFHTH n%conv2nfunction y=circonv

43、2(x1,x2,N)nif length(x1)Nn error(N must not be less than length of x1)nendnif length(x2)Nn error(N must not be less than length of x2)nendX1k=fft(x1,N);X2k=fft(x2,N);Yk=X1k.*X2k;y=ifft(Yk);if(all(imag(x1)=0)&(all(imag(x2)=0) y=real(y);endxn=1,2,1,0;n=size(xn);N=n(2); %取得X(k)的長度subplot(2,1,1);ste

44、m(xn)Xk=fft(xn) %計算FFTX(k)magXk=abs(Xk)subplot(2,1,2)stem(real(Xk); x1=1,1,1,1,0,0,0,0;n=0:16;x2=cos(pi*n/4)+cos(pi*n/8);i=0:7;FFT的應用n凡是利用付里葉變換來進行分析、綜合、變換的地方，都可以利用FFT算法來減少其計算量。nFFT主要應用在：n（1）快速卷積n（2）快速相關n（3）頻譜分析一、快速卷積n設x(n)的長度為N點，h(n)的長度為M點，則：ny(n)的長度為N+M-1點。所以直接計算y(n)線性卷積運算量為NM。1.利用圓周卷積代替線卷積n用圓周卷積（F

45、FT）代替線卷積的必要條件：x(n)、h(n)補零加長至L=N+M-1.n然后計算圓卷積。求出y(n)代表線卷積。1010)()(1010)()(LnNNnnhnhLnNNnnxnx2、用FFT計算y(n)的步驟n(1)求H(k)=DFTh(n) (L點） n(2)求X(k)=DFTx(n) (L點）n(3)計算Y(k)=X(k)H(k) (L點）n(4)求y(n)=IDFTY(k) (L點）2、用FFT計算y(n)與直接計算y(n)的運算量比較稱分段過濾的方法。這時須采用分段卷積或的優(yōu)點。則體現(xiàn)不出其圓周卷積長度差很多時，與）當當（的優(yōu)點?？梢泽w現(xiàn)出其圓周卷積長度差不多時，與）當（計算運算量

46、用直接計算運算量求線卷積運算量用)()(2)()(1log231)1(log23)1(22nhnxnhnxMNNMFFTKLLFFTMNmL3、分段卷積的方法n（1）重疊相加法n（2）重疊保留法二、快速相關1.方法1010)()(1010)()(21)()()(:)()(10*LnMMnnynyLnNNnnxnxmLMNLFFTmnymxnrMnyNnxmLmyx為整數(shù)），令且，來代替線性相關，選擇圓周相關法來求線性相關，是用利用相關點，要求線性的長度為點，的長度為設2.步驟自相關序列。則求得若（而得到。后，取共軛再乘的可以利用即：，來求程序計算同樣可以用已有的點求求乘積點求并求點求),()(

47、),()()1)()(1)(1)()()()()()()()()()()(),()()(*)(*10*)()(nrnrnynxeNFFTkRrWkRNWkRNrrIFFTFFTkRIDFTnrIFFTNdkYkXkRcnyDFTkYFFTNbkXnxDFTkXFFTNaxxxyyxnyxNknkNyxnkNxynyxnyxyxyxyx三、頻譜分析n這里我們僅介紹FFT的儀器設備：快速付里葉分析儀。n其功能為：n（1）能分析確定性信號的頻譜。n（2）估計隨機信號的功率譜n（3）并對信號進行快速卷積濾波n（4）計算相關函數(shù)1.頻譜分析儀的框圖數(shù)據(jù)選擇器A/D保護濾波器輸入電路輸入數(shù)據(jù)存儲器運算器數(shù)

48、據(jù)選擇器控制器變址單元函數(shù)發(fā)生器輸出緩沖器D/A輸出2.部件說明n（1）保護濾波器：對輸入信號進行模擬濾波，濾掉噪聲，提取感興趣的有用信號，以便分析頻譜。n（2）運算器：完成時間抽取FFT或Chirp-Z變換所要求運算。n（3）RAM:存儲輸入數(shù)據(jù)。n（4）ROM(函數(shù)發(fā)生器）。以表格形式存放蝶形運算因子W.n(5)變址單元：根據(jù)輸入數(shù)據(jù)，進行碼位倒置排序。.2)()()2()()1.(2)()(2)(.41IDFTNnxIFFTNkXkXFFTNDFTNnxkXNnx點的實現(xiàn)求點設計用一次，若已知高效算法。完成計算點的設計用一次點的為的有限長序列，是長度為設：習題例)()(21)()()()

49、()(21)()()(1, 2 , 1 , 0),()()()()().()()()(1, 2 , 1 , 0),12()(1, 2 , 1 , 0),2()()()()(122112121212121kNYkYkYnjxDFTkjXkNYkYkYnxDFTkXNknyDFTkYnjxnxnykXkXFFTNDFTnxnxNnnxnxNnnxnxnxnxNnxFFTDITopep令：和求得點可用一次的共軛對稱性，均為實序列，根據(jù)和和點實序列得到兩個點和奇數(shù)點）在時域分別抽取偶數(shù)（解：解題思路是),()()(1, 1 , 0),()()(221221kXWkXNkXNkkXWkXkXkNkN1,

50、 2 , 1 , 0),()(1, 2 , 1 , 0),()(1, 2 , 1 , 0),12()(1, 2 , 1 , 0),2()(2221121NknxDFTkXNknxDFTkXNnnxnxNnnxnx）（kNkNkNWNkXkXkXNkXkXkXkXWkXNkXNkkXWkXkX221221221)()(21)()()(21)(),()()(1, 1 , 0),()()(1, 1 , 0)(Im)(Re)()()()()()()()()()()()()(21)()()(21)()()(2121221NnnyjnykYIFFTnykYkYkjXkXkYkYNkXkXbWNkXkXkX

51、NkXkXkXkXaopepkN點頻域序列構(gòu)成和由計算由運算過程如下：120)21()2()()()()()()()()()(21)(Im)()()()(21)(Re212121NnnnxnnxnxnxnxnxcnjxkYDFTnynynyjnxkYDFTnynynyopep，奇數(shù)，偶數(shù)，合成和由組元素中即可。的數(shù)組的偶數(shù)和奇數(shù)數(shù)依次放入的兩個數(shù)組元素分別和編程時，只要將存放)()()(21nxnxnx習題習題4.10%Xk:被變換數(shù)據(jù)被變換數(shù)據(jù)X(k)%XN:IFFTX(k)結(jié)果結(jié)果x(n)%N:x(n)/X(k)長度長度Xk=X(0) X(1) X(N-1);n=size(Xk);N=n(

52、2); %取得取得X(k)的長度的長度Xk=conj(Xk); %取復共軛取復共軛XN=fft(Xk); %計算計算FFTX(k)XN=conj(XN)/N;stem(real(XN); %繪制繪制x(n)序列波形圖序列波形圖數(shù)據(jù)向量數(shù)據(jù)向量Xk也可通過鍵盤、數(shù)據(jù)文件或函數(shù)計算等多種方法也可通過鍵盤、數(shù)據(jù)文件或函數(shù)計算等多種方法建立，本程序中直接賦值建立，本程序中直接賦值Xk向量。向量。n 如果希望直接調(diào)用FFT子程序計算IFFT，則可用下面的方法：n 由于 10101( )( )1( )( )NknNkNknNkx nX k WNx nXk WN對上式兩邊同時取共軛，得1011( )( )( )NknNkx nXk WDFT XkNN可知：只須將頻域成份一個求共軛變換，即（1）將X(k)