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

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

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

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

4、 F T ) 和 頻 域 抽 取 法和 頻 域 抽 取 法FFT(Decimation In Frequency FFT,簡(jiǎn)簡(jiǎn)稱稱DIFFFT)。下面先介紹。下面先介紹DIFFFT算算法。法。設(shè)序列設(shè)序列x(n)的長(zhǎng)度為的長(zhǎng)度為N，且滿足，且滿足2 ,MNM為自然數(shù) 按按n的奇偶把的奇偶把x(n)分解為兩個(gè)分解為兩個(gè)N/2點(diǎn)的子序列點(diǎn)的子序列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點(diǎn)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 蝶形運(yùn)算符號(hào) CABA BCA BC圖4.2.2 N點(diǎn)DFT的一次時(shí)域抽取分解圖(N=8) N/2點(diǎn)DFTWN0N/2點(diǎn)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)按奇偶分解成兩個(gè)N/4長(zhǎng)的子序列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的對(duì)稱性 Wk+N/4 N/2=-Wk N/2 最后

8、得到： 13/2413/24( )( )( ),0,1,/41(/4)( )( )kNkNX kXkWXkkNX kNXkWXkn 用同樣的方法可計(jì)算出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點(diǎn)DFT的第二次時(shí)域抽取分解圖(N=8) N/4點(diǎn)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點(diǎn)DFTN/4點(diǎn)DFTN/4點(diǎn)DFTWN02WN02 圖4.2.4 N點(diǎn)DITFFT運(yùn)算流圖(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算法與直接計(jì)算DFT運(yùn)算量的比較n 每一級(jí)運(yùn)算都需要N/2次復(fù)數(shù)乘和N次復(fù)數(shù)加(每個(gè)蝶形需要兩次復(fù)數(shù)加法)。所以，M級(jí)運(yùn)算總共需要的復(fù)數(shù)乘次數(shù)為22(2)log22(2)logMANNCMNCN MNN復(fù)數(shù)加次數(shù)為 例如，N=210=1024時(shí)221048576204.8(/2)log5120NNN圖4.2.5 FFT算法與直接計(jì)算DFT所需乘法次數(shù)的比較曲線 由于N=2L，因此N/2仍為偶

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

12、信號(hào)需要進(jìn)行如果在某一支路上信號(hào)需要進(jìn)行相乘運(yùn)算，則在該支路上標(biāo)以箭頭，相乘運(yùn)算，則在該支路上標(biāo)以箭頭，將相乘的系數(shù)標(biāo)在箭頭旁。將相乘的系數(shù)標(biāo)在箭頭旁。(5)當(dāng)支路上沒(méi)有箭頭及系當(dāng)支路上沒(méi)有箭頭及系數(shù)時(shí)，則該支路的傳輸比數(shù)時(shí)，則該支路的傳輸比為為1。例子：求 N=23=8點(diǎn)FFT變換 （1）先按）先按N=8-N/2=4，做做4點(diǎn)的點(diǎn)的DFT：)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN12/, 0Nk先將N=8 DFT分解成2個(gè)4點(diǎn)DFT：可知：時(shí)域上：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點(diǎn)直接DFT與分解2個(gè)4點(diǎn)DFT的FFT運(yùn)算量N=8點(diǎn)的直接DFT的計(jì)算量為：N2次(64次）復(fù)數(shù)相乘,N(N-1)次(8(8-1)=56次)復(fù)數(shù)相加.共計(jì)120次。（b)求 一個(gè)蝶形結(jié)需要的運(yùn)算量要運(yùn)算一個(gè)蝶形結(jié)，需要一次乘法 ， 兩次加法。)(2kXWkN)()() 2/()()()(2121kXWkXNkXkXWkXkXkNkN（c）分解為兩個(gè)N/2=4點(diǎn)的DFT的運(yùn)算量分解2個(gè)N/2點(diǎn)（4點(diǎn)）的DFT：) 12()2()rxDFTrxDFTkX（偶數(shù)其復(fù)數(shù)相乘為復(fù)數(shù)相加為奇數(shù)其復(fù)數(shù)相乘為復(fù)數(shù)相加為2

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

15、的信號(hào)流圖 4點(diǎn)DFTx(0)x(2)x(4)x(6)4點(diǎn)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)同學(xué)們自已寫x1(r)x2(r)(2)N/2(4點(diǎn))-N/4(2點(diǎn))FFT(a)先將先將4點(diǎn)分解成點(diǎn)分解成2點(diǎn)的點(diǎn)的DFT：n因?yàn)?點(diǎn)DFT

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

17、(62/5262/52242/3142/314/0) 12(2/114/022/1111LXWLXNkxLXWLXkxrxLXWLXNkXLLXWLXWLXWLXkXkXrxkNkNkNkNNLkLNNLLkNDFT）也可分解為：同樣，（同理：，其中（可分解為：（c)一個(gè)2點(diǎn)的DFT蝶形流圖2點(diǎn)DFT2點(diǎn)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)另一個(gè)2點(diǎn)的DFT蝶形流圖2點(diǎn)DFT2點(diǎn)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點(diǎn)）DFT再分解成2個(gè)1點(diǎn)的DFT(a)求2個(gè)一點(diǎn)的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0

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

20、T的運(yùn)算流圖 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的運(yùn)算規(guī)律及編程思想n 1.原位計(jì)算n 由圖4.2.4可以看出，DITFFT的運(yùn)算過(guò)程很有規(guī)律。N=2M點(diǎn)的FFT共進(jìn)行M級(jí)運(yùn)算，每級(jí)由N/2個(gè)蝶形運(yùn)算組成。同一級(jí)中，每個(gè)蝶形的兩個(gè)輸入數(shù)據(jù)只對(duì)本蝶形有用，且每個(gè)蝶形的輸入、輸出數(shù)據(jù)結(jié)點(diǎn)又在同一水平線上，計(jì)算完一個(gè)蝶形后，所得輸出數(shù)據(jù)可立即存入原數(shù)據(jù)占用

21、的存儲(chǔ)單元，這就是位計(jì)算。原位計(jì)算可大節(jié)省內(nèi)存單元，降低成本。n 原位運(yùn)算（in-place)n原位運(yùn)算的結(jié)構(gòu)，可以節(jié)省存儲(chǔ)單元，降低設(shè)備成本。n定義：當(dāng)數(shù)據(jù)輸入到存儲(chǔ)器以后，每一組運(yùn)算的結(jié)果，仍然存放在這同一組存儲(chǔ)器中直到最后輸出。02Wx(0)x(4)X3(0)X3(1)中放在一個(gè)暫存器中放在單元中，將放在單元例：將RWAxAx08) 1 ()4()0()0(單元送回單元送回將) 1 ()4()0()0()4()0(0808AxWxAxWxn例：N=8 FFT運(yùn)算,輸入: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看出：用原位運(yùn)算結(jié)構(gòu)運(yùn)算后，A(0)A(7)正好順序存放X(0)X(7),可以直接順序輸出。2.旋轉(zhuǎn)因子的變化規(guī)律 如上所述，N點(diǎn)DITFFT運(yùn)算流圖中，每級(jí)都有N/2個(gè)蝶形。每個(gè)蝶形都要乘以因子WpN，稱其為旋轉(zhuǎn)因子，

23、p稱為旋轉(zhuǎn)因子的指數(shù)。 觀察圖4.2.4不難發(fā)現(xiàn)，第L級(jí)共有2L-1個(gè)不同的旋轉(zhuǎn)因子。N=23=8時(shí)的各級(jí)旋轉(zhuǎn)因子表示如下： L=1時(shí)，WpN=WJ N/4=WJ2L， J=0 L=2時(shí)， WpN =WJ N/2=WJ2L， J=0，1 L=3時(shí)， 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對(duì)N=2M的一般情

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

25、表示取實(shí)部，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運(yùn)算和程序框圖 開(kāi) 始送入x(n)， MN2 M倒 序L1 , M0 , B 1P2 M LJk J , N1 , 2LpNpNWBkXkXBkXWBkXkXkX)()()

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

28、表4.2.1 順序和倒序二進(jìn)制數(shù)對(duì)照表 圖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í)行時(shí)，只須將1與4對(duì)調(diào)送入，3與6對(duì)調(diào)送入，而0，2，5，7不變，僅需要一個(gè)中間存儲(chǔ)單元。I01234567J04261537 在實(shí)際運(yùn)算時(shí)，先按自然順序?qū)⑤斎胄蛄写嫒氪鎯?chǔ)單元，再通過(guò)變址運(yùn)算將自然順序變換成按時(shí)間抽取的FFT算法要求的順序。變址的過(guò)程可以用程

29、序安排加以實(shí)現(xiàn)-稱為“整序”或“重排”（采用碼位倒讀）且注意：(1)當(dāng)I=J時(shí)，數(shù)據(jù)不必調(diào)換；(2)當(dāng)IJ時(shí)，必須將原來(lái)存放數(shù)據(jù)x(n)送入暫存器R，再將x(J)送入x(I),R中內(nèi)容送x(J).進(jìn)行數(shù)據(jù)對(duì)調(diào)。(3)為了避免再次考慮前面已調(diào)換過(guò)的數(shù)據(jù)，保證調(diào)換只進(jìn)行一次，否則又變回原狀。IN/2=4，做做4點(diǎn)的點(diǎn)的DIF：nNWNnxnxnxNnxnxnx)2()()()2()()(21先將N=8DFT分解成2個(gè)4點(diǎn)DFT：可知：時(shí)域上：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點(diǎn)分解成2個(gè)4點(diǎn)的DIF的信號(hào)流圖 4點(diǎn)DFTx(0)x(1)x(2)x(3)4點(diǎn)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點(diǎn))-N/4(2點(diǎn))FFT(a)先將先將4點(diǎn)分解成點(diǎn)分解成2點(diǎn)的點(diǎn)的DIF：n因?yàn)?點(diǎn)DIF還是比較麻煩，所以再繼續(xù)分解。后半部分序列、前半部分序列、) 3 () 2() 1 () 0(: )(11111xxxxnx10140)2()()()2()(411311，），在此（）（若設(shè)：LNLLxWNnxnxLxNnxnxnN后半部分序列、前半部分序列、同理：) 3 () 2() 1 () 0(: )(22222xxxxnx10140)2()()()2()(622522，），在此（）（同理：LNLLxWNnxnxLxNnxnxnN（b)一個(gè)2點(diǎn)的DIF蝶形流圖2點(diǎn)DFT2點(diǎn)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)另一個(gè)2點(diǎn)的DIF蝶形流圖2點(diǎn)DFT2點(diǎn)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點(diǎn)）DFT再分解成2個(gè)1點(diǎn)的DFT(a)求2個(gè)一點(diǎn)的DFT021022120202120230202020231021 , 0; 1 , 0)4()0()4()0() 1 ()4()0()4()0()0

33、()()(2WWknWWWWWxWxWxWxXWxWxWxWxXWnxkXnknknkNnkNnnkN，其中，則這里用到對(duì)稱性這是一蝶形結(jié)代入上面流圖可知： 最后剩下兩點(diǎn)DFT,它可分解成兩個(gè)一點(diǎn)DFT，但一點(diǎn)DFT就等于輸入信號(hào)本身，所以兩點(diǎn)DFT可用一個(gè)蝶形結(jié)表示。取x3(0)、x3(1)為例。(b)2個(gè)1點(diǎn)的DFT蝶形流圖 1點(diǎn)DFTx3(0)1點(diǎn)DFTx3(1)X(0)X(4)02W進(jìn)一步簡(jiǎn)化為蝶形流圖：02WX(0)X(4)x3(0)x3(1)(4)一個(gè)完整N=8的按頻率抽取FFT的運(yùn)算流圖 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的一種變形運(yùn)算流圖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按時(shí)間抽取的FFT算法的若干變體 對(duì)于任何流圖，只要保持各節(jié)點(diǎn)所連續(xù)的支路及其傳輸系數(shù)不變，則不論節(jié)點(diǎn)位置怎么排列所得流圖總是等效的，最后所得結(jié)果都是x(n)的離散付里葉變換的正確結(jié)果。只是數(shù)

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

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

37、TFFTNWWDFTWnxnxDFTkXWkXNkXIDFTnxnkNnkNNknkNNknkN算法都可以拿來(lái)運(yùn)算或頻率抽取抽取）那么以上討論的時(shí)間（將運(yùn)算結(jié)果都除以改成運(yùn)算中的每個(gè)系數(shù)只要把3)2() 1 ()()()()(1)()(1010圖4.2.16 DITIFFT運(yùn)算流圖 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運(yùn)算流圖(防止溢出) 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子程序計(jì)算IFFT，則可用下面的方法：n 由于 10101( )( )1( )( )NknNkNknNkx nX k WNx nXk WN對(duì)上式兩邊同時(shí)取共軛，得1011( )( )( )NknNkx nXk WDFT XkNN可知：只須將頻域成份一個(gè)求共軛變換，即（1）將X(k)的虛部乘以-1，即先取X(k)的共軛，得X*(k)。

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

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

41、4(3)2(2)2213(2)102MNNRMNMn 4.3.2 旋轉(zhuǎn)因子的生成n 在FFT運(yùn)算中，旋轉(zhuǎn)因子WmN=cos(2m/N)-jsin(2m/N)，求正弦和余弦函數(shù)值的計(jì)算量是很大的。 n 4.3.3 實(shí)序列的FFT算法 n 設(shè)x(n)為N點(diǎn)實(shí)序列，取x(n)的偶數(shù)點(diǎn)和奇數(shù)點(diǎn)分別作為新構(gòu)造序列y(n)的實(shí)部和虛部，即1212( )(2 ),( )(21),0,1,12( )( )( ),0,1,12Nx nxnx nxnnNy nx njx n n對(duì)y(n)進(jìn)行N/2點(diǎn)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實(shí)信號(hào)的快速循環(huán)卷積實(shí)信號(hào)的快速循環(huán)卷積n 用DITFFT計(jì)算兩個(gè)實(shí)信號(hào)x1(n)和x2(n)的循環(huán)卷積時(shí)，最直接的方法是將信號(hào)的虛部都置為零，再按復(fù)序列用FFT計(jì)算。顯然這樣做是很浪費(fèi)的?，F(xiàn)在我們可用基2DITFHT直接進(jìn)行實(shí)正交變換來(lái)處理實(shí)信號(hào)的循環(huán)卷積問(wèn)題。1122( )( )( )( )FHTHFHTHx nXkx nXk計(jì)算式 )()()()(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)的長(zhǎng)度subplot(2,1,1);ste

44、m(xn)Xk=fft(xn) %計(jì)算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的應(yīng)用n凡是利用付里葉變換來(lái)進(jìn)行分析、綜合、變換的地方，都可以利用FFT算法來(lái)減少其計(jì)算量。nFFT主要應(yīng)用在：n（1）快速卷積n（2）快速相關(guān)n（3）頻譜分析一、快速卷積n設(shè)x(n)的長(zhǎng)度為N點(diǎn)，h(n)的長(zhǎng)度為M點(diǎn)，則：ny(n)的長(zhǎng)度為N+M-1點(diǎn)。所以直接計(jì)算y(n)線性卷積運(yùn)算量為NM。1.利用圓周卷積代替線卷積n用圓周卷積（F

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

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

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

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

49、()(21)()()(1, 2 , 1 , 0),()()()()().()()()(1, 2 , 1 , 0),12()(1, 2 , 1 , 0),2()()()()(122112121212121kNYkYkYnjxDFTkjXkNYkYkYnxDFTkXNknyDFTkYnjxnxnykXkXFFTNDFTnxnxNnnxnxNnnxnxnxnxNnxFFTDITopep令：和求得點(diǎn)可用一次的共軛對(duì)稱性，均為實(shí)序列，根據(jù)和和點(diǎn)實(shí)序列得到兩個(gè)點(diǎn)和奇數(shù)點(diǎn)）在時(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點(diǎn)頻域序列構(gòu)成和由計(jì)算由運(yùn)算過(guò)程如下：120)21()2()()()()()()()()()(21)(Im)()()()(21)(Re212121NnnnxnnxnxnxnxnxcnjxkYDFTnynynyjnxkYDFTnynynyopep，奇數(shù)，偶數(shù)，合成和由組元素中即可。的數(shù)組的偶數(shù)和奇數(shù)數(shù)依次放入的兩個(gè)數(shù)組元素分別和編程時(shí)，只要將存放)()()(21nxnxnx習(xí)題習(xí)題4.10%Xk:被變換數(shù)據(jù)被變換數(shù)據(jù)X(k)%XN:IFFTX(k)結(jié)果結(jié)果x(n)%N:x(n)/X(k)長(zhǎng)度長(zhǎng)度Xk=X(0) X(1) X(N-1);n=size(Xk);N=n(

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