數(shù)字信號(hào)處理4人文科技ppt課件

1、 第第4 4章章 FFT FFT 4.1 4.1 引言引言 4.1.1 4.1.1 離散傅里葉變換的矩陣表示及其運(yùn)算量離散傅里葉變換的矩陣表示及其運(yùn)算量 DFT DFT在數(shù)字信號(hào)處理中起著非常重要的作用，在數(shù)字信號(hào)處理中起著非常重要的作用， 這是與這是與DFTDFT存在著高效算法，存在著高效算法， 即快速傅里葉變換即快速傅里葉變換(FFT) (FFT) 分不開的。快速運(yùn)算的關(guān)鍵是減少運(yùn)算量。分不開的?？焖龠\(yùn)算的關(guān)鍵是減少運(yùn)算量。n離散傅里葉變換對(duì)為：n (4.1)n (4.2)n式中 。 下面要用矩陣來表示DFT關(guān)系。令： 101, 1 , 0,)()(:NnnkNNkWnxkXDFT1, 1

2、 , 0,)(1)(:10NnWkXNnxIDFTNknkNNjNeW2) 1() 1 ()0(Nxxxx) 1() 1 ()0(NXXXX并令TN、 表示兩個(gè)變換方陣，有：1NT2)1(2)1(1)1(22221211111111NNNNNNNNNNNNNNNWWWWWWWWWT2)1(2)1()1()1(2222)1(2111111111nNNNNNNNNNNNNNNWWWWWWWWWTn如果用i (0iN-1)表示這兩個(gè)N階方陣的行號(hào)，用j(0jN-1)表示這兩個(gè)N階方陣的列號(hào)，那末很容易看出，TN 方陣中i行j列的元素為 ，而方陣中i行j列的元素為 。n 于是4.1式可以寫成： n而4

3、.2式可以寫成： jiNWjiNWxTXN 11XTNxNn 一般情況下，信號(hào)序列x(n) 及其頻譜序列X(k) 都是用復(fù)數(shù)來表示的，WN當(dāng)然也是復(fù)數(shù)。因此，計(jì)算DFT的一個(gè)值X(k) 需要進(jìn)行N次復(fù)數(shù)乘法與1相乘也包括在內(nèi)和N-1次復(fù)數(shù)加法；所以，直接計(jì)算N點(diǎn)的DFT需要進(jìn)行N2 次復(fù)數(shù)乘法和N(N-1) 復(fù)數(shù)加法。 n顯然，直接計(jì)算N點(diǎn)的IDFT所需的復(fù)乘和復(fù)加的次數(shù)也是這么多。當(dāng)N足夠大時(shí)，N2 N(N-1), 因此，DFT與IDFT的運(yùn)算次數(shù)與N2 成正比，隨著N的增加，運(yùn)算量將急劇增加，而在實(shí)際問題中，N往往是較大的，因此有必要對(duì)DFT與IDFT的計(jì)算方法予以改進(jìn)。 4.1.2 因子

4、的特性因子的特性 DFT和和IDFT的快速算法的導(dǎo)出主要是根據(jù)的快速算法的導(dǎo)出主要是根據(jù) 因子的特性。因子的特性。 1周期性：周期性： 對(duì)離散變量對(duì)離散變量n有同樣的周期性。有同樣的周期性。nkNnNNnkNNknNWWWW)(nkNWnkNW 2對(duì)稱性：對(duì)稱性： 或或 3. 其它其它: knNNknNnkNnkNWWWW)()(*nkNNnkNknNnkNWWWW)()(*kNNNjkNNNkNNkNWeWWWW222)2(kNkNjkNjkNWeeW 2222224.2 4.2 基基2 2時(shí)間抽選的時(shí)間抽選的FFTFFT算法算法 4.2.1 4.2.1 算法推導(dǎo)算法推導(dǎo) 已經(jīng)知道：已經(jīng)知道

5、： 令令DFTDFT的長(zhǎng)度的長(zhǎng)度N=2MN=2M，M M為正整數(shù)。為正整數(shù)。1-N,0,1,k )()()(10nkNNnWnxnxDFTkXn令： n于是有：12N,0,1,r ) 12()()2()(rxrprxrg1-N,0,1,k )()()()() 12()2()(120rk 2/120 2/120)12(1202kPWkGWrpWWrgWrxWrxkXkNNrNkNNrrkNNrkrNNrrkNn其中， n是由x(n)的偶數(shù)抽樣點(diǎn)形成的DFT；而n kr 2/120120kr 2/)2()()(NNrNrNWrxWrgkG120kr 2/kr 2/120) 12()()(NrNNN

6、rWrxWnpkPn是由x(n)的奇數(shù)抽樣點(diǎn)形成的DFT。但是這兩個(gè)式子并不完全是N/2點(diǎn)的DFT，因?yàn)閗的范圍仍然是由0到N-1，因此，還應(yīng)該進(jìn)一步考慮k由N/2到N-1范圍的情況。n現(xiàn)在令 ，故對(duì)于后半段有：n n同理： n又知： 12, 1 , 0Nk)()()()()2(120kr 21202120)2(2/22kGWrgWWrgWrgNkGNrNNrrkNrNrNkrNNN)()2(kPNkPk 2NNkNWW 圖圖 4.1 N點(diǎn)點(diǎn)DFT分解為兩個(gè)分解為兩個(gè)N/2點(diǎn)的點(diǎn)的DFT (N=8) 圖圖 4.2 N/2 點(diǎn)的點(diǎn)的DFT 分解為兩個(gè)分解為兩個(gè)N/4點(diǎn)的點(diǎn)的DFT (N=8)n綜

7、上所述，可以得到： n其中G(k)、P(k) 分別是x(n)的偶數(shù)點(diǎn)和奇數(shù)點(diǎn)的N/2點(diǎn)DFT。 )()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn 這樣，我們就將一個(gè)N點(diǎn)的DFT分解成了兩個(gè)N/2點(diǎn)的DFT，由于DFT的運(yùn)算量與其點(diǎn)數(shù)的平方成正比，因此使運(yùn)算量減少了。但是，還應(yīng)該將每一個(gè)N/2點(diǎn)的DFT再分解為兩個(gè)N/4點(diǎn)的DFT，如此下去，直到分解為2點(diǎn)的DFT為止，總共需要進(jìn)行l(wèi)og2N-1=log2(N/2)次分解。 圖圖 4.3 2 點(diǎn)點(diǎn) DFT 信號(hào)流圖蝶形圖）信號(hào)流圖蝶形圖）n對(duì)于2點(diǎn)DFT，有： n所以2點(diǎn)DFT的運(yùn)算只需一次加法和一次減法

8、，這樣的運(yùn)算叫做蝶形運(yùn)算，這樣的信號(hào)流圖叫做蝶形圖。1111111122WTn該算法每次分解都是將時(shí)域序列按奇偶分為兩組，因此要求N等于2的正整數(shù)冪，故將這種FFT算法叫做基2時(shí)間抽選法。 4.2.2 算法特點(diǎn) 1. 倒序重排這種FFT算法的每次分解都是將輸入序列按照奇偶分為兩組，故要不斷地將每組輸入數(shù)據(jù)按奇偶重排，直到最后分解為2點(diǎn)的DFT，輸入數(shù)據(jù)才不再改變順序。這樣做的結(jié)果，使得作FFT運(yùn)算時(shí)，輸入序列的次序要按其序號(hào)的倒序進(jìn)行重新排列。 n現(xiàn)在將圖4.4中輸入序號(hào)以及重排后的序號(hào)按二進(jìn)制寫出如下注：下標(biāo)“2表示二進(jìn)制數(shù)）?？梢钥闯?，將輸入序號(hào)的二進(jìn)制表示n2n1n0位置顛倒，得到n0n

9、1n2），就是相應(yīng)的倒序的二進(jìn)制序號(hào)。因此，輸入序列按倒序重排，實(shí)際上就是將序號(hào)為n2n1n0的元素與序號(hào)為n0n1n2的元素的位置相互交換。 2. 同址計(jì)算同址計(jì)算 從圖從圖4.4可以看出，整個(gè)算法流圖可以分為四段，（可以看出，整個(gè)算法流圖可以分為四段，（0段段為倒序重排，后面為倒序重排，后面3段為段為3(log28=3)次迭代運(yùn)算：首先計(jì)算次迭代運(yùn)算：首先計(jì)算2點(diǎn)點(diǎn)DFT，然后將，然后將2點(diǎn)點(diǎn)DFT的結(jié)果組合成的結(jié)果組合成4點(diǎn)點(diǎn)DFT，最后將，最后將4點(diǎn)點(diǎn)DFT組合為組合為8點(diǎn)點(diǎn)DFT。因此，對(duì)于。因此，對(duì)于N點(diǎn)點(diǎn)FFT，只需要一列，只需要一列存儲(chǔ)存儲(chǔ)N個(gè)復(fù)數(shù)的存儲(chǔ)器。個(gè)復(fù)數(shù)的存儲(chǔ)器。 n

10、開始時(shí)，輸入序列的N個(gè)數(shù)放于此N個(gè)存儲(chǔ)器內(nèi)，倒序重排后仍存于這N個(gè)存儲(chǔ)器中，每一次迭代運(yùn)算后的結(jié)果也仍然存于這N個(gè)存儲(chǔ)器中，整個(gè)運(yùn)算完成后，這N個(gè)存儲(chǔ)器中即為所求的頻譜序列X(k) （k = 0、1、.、 N-1）。這就是所謂的同址計(jì)算，這樣可以大大節(jié)約存儲(chǔ)元件。 3. 運(yùn)算量運(yùn)算量觀察圖觀察圖4.4可知，圖可知，圖4.3所示的蝶形圖實(shí)際上代表了所示的蝶形圖實(shí)際上代表了FFT的基的基本運(yùn)算，它實(shí)際上只包含了兩次復(fù)數(shù)加法運(yùn)算。一個(gè)本運(yùn)算，它實(shí)際上只包含了兩次復(fù)數(shù)加法運(yùn)算。一個(gè)N(N=2M) 點(diǎn)的點(diǎn)的FFT，需要進(jìn)行，需要進(jìn)行M=log2N次迭代運(yùn)算。每次迭代運(yùn)算。每次迭代運(yùn)算包含了次迭代運(yùn)算包含

11、了N/2個(gè)蝶型，因此共有個(gè)蝶型，因此共有N次復(fù)數(shù)加法；次復(fù)數(shù)加法；此外，除了第一次的此外，除了第一次的2點(diǎn)點(diǎn)DFT之外，每次迭代還包括了之外，每次迭代還包括了N/2次復(fù)數(shù)乘法即乘次復(fù)數(shù)乘法即乘WN的冪）。因此，一個(gè)的冪）。因此，一個(gè)N點(diǎn)的點(diǎn)的FFT共共有復(fù)數(shù)乘法的次數(shù)為：有復(fù)數(shù)乘法的次數(shù)為： 復(fù)數(shù)加法的次數(shù)為： 因此，F(xiàn)FT算法比直接計(jì)算DFT大大減少了運(yùn)算量，尤其是當(dāng)N較大時(shí)，計(jì)算量的減少更為顯著。比如，當(dāng)N=1024時(shí)，采用FFT算法時(shí)復(fù)數(shù)乘法的次數(shù)，不超過直接計(jì)算DFT時(shí)復(fù)乘次數(shù)的千分之五。 2log2) 1(log2222NNNNMcNNNNAc222loglog22 4.2.3 關(guān)于

12、FFT算法的計(jì)算機(jī)程序在一般情況下，進(jìn)行FFT運(yùn)算的序列至少都有幾百點(diǎn)的長(zhǎng)度，因此需要編制FFT算法程序以便能夠利用計(jì)算機(jī)來快速進(jìn)行計(jì)算?？梢詤⒄?qǐng)D4.4來編制計(jì)算N (N必須等于2的正整數(shù)冪)點(diǎn)FFT的計(jì)算機(jī)程序，整個(gè)程序可以分為兩部分：一部分是倒序重排，另一部分是用三層嵌套的循環(huán)來完成M=log2N次迭代。 n三層循環(huán)的功能是：最里的一層循環(huán)完成蝶形運(yùn)算，中間的一層循環(huán)完成因子WNk的變化，而最外的一層循環(huán)則是完成M次迭代過程。n倒序重排的程序是一段經(jīng)典程序，它以巧妙的構(gòu)思、簡(jiǎn)單的語句用高級(jí)編程語言來完成了倒序重排的功能。下面是一段用FORTRAN語言編寫的倒序重排程序。 N=2*M （表

13、示N=2M ，M是輸入的正整數(shù)） NV2=N/2 （NV2是一個(gè)整數(shù)變量） NM1=N-1 （NM1也是一個(gè)整數(shù)變量） J=1 （對(duì)變量J賦初值） 100 DO 7 I=1,NM1 (循環(huán)開始，到語句7結(jié)束； 循環(huán)變量I從1開始，到NM1結(jié)束，步長(zhǎng)為1)n IF (I.GE.J) GOTO 5 （如果IJ，就轉(zhuǎn)移到語句5）n （將輸入序列中序號(hào)互為倒序n 的兩個(gè)數(shù)值交換位置） TIXIXJXJXT)( )()()( 5 K=NV2 6 IF(K.GE.J) GOTO 7 (如果KJ，就轉(zhuǎn)移到語句 7)J=J-kK=K/2 GOTO 6 （轉(zhuǎn)移到語句6）7 J=J+K n由于是同址計(jì)算，因此程序

14、中只用了一個(gè)數(shù)組X( )，這個(gè)數(shù)組共有N個(gè)元素。X( )既表示輸入序列，也表示按倒序重排后的這N個(gè)數(shù)值。程序中的字母都是大寫的，這是FORTRAN語句的特點(diǎn)。nI既是循環(huán)變量，也代表輸入序列的正序序號(hào)，J代表倒置后的序號(hào)。由于FORTRAN語言的循環(huán)變量不能從0開始，故以8點(diǎn)FFT為例，I的范圍為18，下面是正序序號(hào)I與倒序序號(hào)J之間的對(duì)應(yīng)關(guān)系：n I: 1 2 3 4 5 6 7 8 n J: 1 5 3 7 2 6 4 8n 輸入序列按倒序重排，實(shí)際上就是將X(I)和X(J)互換位置。顯然，如果I=J，就不需要交換；而如果IJ，就表示X(I)與X(J)已經(jīng)交換過了。因此，在循環(huán)語句下面的語

15、句“IF (I.GE.J) GOTO 5就表明了交換的條件：只有當(dāng)IJ時(shí)才執(zhí)行下面三條語句，使X(I)與X(J)互換位置。n正序序號(hào)I也是循環(huán)變量，從1開始，每次循環(huán)增加1。關(guān)鍵問題是對(duì)于每一個(gè)I，所對(duì)應(yīng)的J是什么？程序中從語句5開始直到循環(huán)結(jié)束就是為下一次循環(huán)確定所對(duì)應(yīng)的J的。雖然I-1和J-1的二進(jìn)制表示是互為倒序的關(guān)系，但是，程序中并不是由I得到J，而是由上一個(gè)J來得到下一個(gè)J。思路是：I-1由0到7，用二進(jìn)制來表示就是從0002開始，每次在最低位加1，逐次增加到1112 。 n既然J-1的二進(jìn)制表示是I-1二進(jìn)制表示的倒置，那末J-1的變化就應(yīng)該是在最高位逐次加1，而最高位的二進(jìn)制數(shù)1

16、表示N/2。所以，程序中將J與N/2比較來判斷J-1的最高位是0還是1，如果是0，就執(zhí)行J+K=J+N/2，也就是將J-1的二進(jìn)制最高位由0變?yōu)?。如果J-1的最高位是1，應(yīng)該將它變?yōu)?，也就是將J減去N/2，然后考察次高位是0還是1，這應(yīng)該用N/4來考察，于是執(zhí)行“K=K/2”。如此下去，直到得到下一個(gè)J，完成此次循環(huán)。4.3 4.3 基基2 2頻率抽選的頻率抽選的FFTFFT算法算法 時(shí)間抽選法是在時(shí)域內(nèi)將輸入序列時(shí)間抽選法是在時(shí)域內(nèi)將輸入序列x(n)x(n)逐次分解為偶數(shù)逐次分解為偶數(shù)點(diǎn)子序列和奇數(shù)點(diǎn)子序列，通過求子序列的點(diǎn)子序列和奇數(shù)點(diǎn)子序列，通過求子序列的DFTDFT而實(shí)現(xiàn)整而實(shí)現(xiàn)整

17、個(gè)序列的個(gè)序列的DFTDFT。而頻率抽選法則是在頻域內(nèi)將。而頻率抽選法則是在頻域內(nèi)將X(k)X(k)逐次分逐次分解成偶數(shù)點(diǎn)子序列和奇數(shù)點(diǎn)子序列，然后對(duì)這些分解得越解成偶數(shù)點(diǎn)子序列和奇數(shù)點(diǎn)子序列，然后對(duì)這些分解得越來越短的子序列進(jìn)行來越短的子序列進(jìn)行DFTDFT運(yùn)算，從而求得整個(gè)的運(yùn)算，從而求得整個(gè)的DFTDFT。當(dāng)然，。當(dāng)然，同樣要求同樣要求N N為為2 2的正整數(shù)冪。的正整數(shù)冪。 n 設(shè) ，則可以分別表示出k為偶數(shù)和奇數(shù)時(shí)的X(k)。n n其中， 12, 1 , 0Nr120rn 2/120rn 2/120rn 2/120120)2(22102)()2()()2()()()2(NnNNnrN

18、NNNnNNnNnrNnNnrNNnrnNWngWWNnxWnxWNnxWnxWnxrX11)2()()(2N,0,n Nnxnxng120rn 2/120rn 2/120rn 2/12022120rn 2/120120)2)(12(2)() 1()2()()2()()2()() 12 (NnNnNNnnNNNnnNNNnNNrNNnNnrNNnnNNNnNnNnrNnNnrNWWnpWWNnxWWnxWWWWNnxWWnxWNnxWWnxrXn其中， n 顯然，X(2r) 為g(n) 的N/2點(diǎn)DFT，X(2r+1) 為p(n)WNn 的N/2點(diǎn)DFT。這樣，一個(gè)N點(diǎn)的DFT分解為兩個(gè)N/2

19、點(diǎn)的DFT。將分解繼續(xù)下去，直到分解為2點(diǎn)的DFT為止。當(dāng)N=8時(shí)，基2頻率抽選的FFT算法的整個(gè)信號(hào)流圖如圖4.6所示。11)2()()(2N,0,n Nnxnxnpn將圖4.6與圖4.4比較，可知頻率抽選法的計(jì)算量與時(shí)間抽選法相同，而且都能夠同址計(jì)算。時(shí)間抽選法是輸入序列按奇偶分組，故x(n)的順序要按倒序重排，而輸出序列按前后分半，故X(k) 的順序不需要重排；頻率抽選法則是輸出序列按奇偶分組，故X(k) 的順序要按倒序重排，而輸入序列按前后分半，故x(n) 不需要重排。4.4 4.4 基基 4 4時(shí)間抽選的時(shí)間抽選的FFTFFT算法算法 上面討論的上面討論的FFTFFT算法都要求算法都

20、要求N=2M N=2M ，M M是正整數(shù)，因此是是正整數(shù)，因此是基基2 2的的FFTFFT算法。若算法。若N=4M N=4M ，則還可以采用基，則還可以采用基4 4的的FFTFFT算法。算法。與推導(dǎo)基與推導(dǎo)基2 FFT2 FFT算法類似，可以推導(dǎo)基算法類似，可以推導(dǎo)基4 FFT4 FFT算法。算法。將將x(n) x(n) 相間地分為四組，所得的相間地分為四組，所得的X(k) X(k) 按前后分為四組。第按前后分為四組。第一次分組的流圖如圖一次分組的流圖如圖4.74.7所示。其中，中間的矩陣為變換所示。其中，中間的矩陣為變換方陣方陣T4T4， j- 1- j 11- 1 1- 1j 1- j-

21、11 1 1 111111119464346444243424144WWWWWWWWWT而 Di (i = 0,1,2,3) 則為N/4 階的對(duì)角矩陣，即： 0,1,2,3i WWWDiNNiNiNi )14(21 即變換方陣，即： 4/NT2)14N(N/42)14N(N/41)14N(N/4)14N2(N/422N/412 N/414NN/42 N/41 N/44 W W W1 W W W1 W W W11 1 1 1NT 圖圖 4.7 基基 4 時(shí)間抽選的時(shí)間抽選的 FFT算法的第一次分解算法的第一次分解n為了理解基4 FFT算法，可以將基2時(shí)間抽選的FFT算法的第一次分解與基4的第一次

22、分解進(jìn)行比較。n )()()()(1- 11 1)2()(2kPWkGTkPWkGNkXkXkNkN12, 1 , 0Nkn基2的情形第一次分解為兩組，基4為四組?；?分解中的G(k) 和P(k) 分別為x(n)的偶數(shù)點(diǎn)和奇數(shù)點(diǎn)的N/2點(diǎn)DFT，而基4分解的每一項(xiàng)中的因子n ) () () (4/xxxTNn都代表一個(gè)N/4點(diǎn)的DFT，x(n) 的間隔為4。基2中的前后兩項(xiàng)用T2矩陣來連接；而在基4的情形為4項(xiàng)，是用T4矩陣來連接的?；?中的第二項(xiàng)的N/2點(diǎn)的DFT之前要乘上因子WNk；而基4的情況每個(gè)N/4點(diǎn)的DFT之前則要乘上對(duì)角矩陣Di，其對(duì)角線上的元素為WNki。因此，基4的FFT分解

23、既與基2 FFT的分解相對(duì)應(yīng)，又比之復(fù)雜。n基4的情形也要繼續(xù)分解下去，直到分解為四點(diǎn)DFT為止，共要經(jīng)過log4N-1= log4(N/4) 次分解。估算基4 FFT的計(jì)算量：一個(gè)N點(diǎn)的FFT，要經(jīng)過log4N次變換。復(fù)數(shù)乘法體現(xiàn)在與對(duì)角矩陣D1、D2、D3相乘，總的復(fù)乘次數(shù)為： ；總的復(fù)加次數(shù)為： 。 ) 1(log4344NNMcNNAc44log3n 基2 FFT算法的復(fù)乘次數(shù)為：n n復(fù)加次數(shù)為： )21(log2log222log24422NNNNNNMcNNNNAc422log2logn將Mc4與Mc2比較，Ac4和Ac2比較，可知，基4的復(fù)乘次數(shù)約減少了1/4，但復(fù)加次數(shù)卻增加

24、了。因此，在乘法運(yùn)算要轉(zhuǎn)化為加法運(yùn)算來實(shí)現(xiàn)的情況下，基4算法的計(jì)算量會(huì)比基2算法少；但是，近年來，在許多微處理器和DSP中，實(shí)現(xiàn)一次乘法運(yùn)算與實(shí)現(xiàn)一次加法運(yùn)算的時(shí)間完全相同，這時(shí)采用基4算法就不合算了，而且基4算法還比基2算法復(fù)雜。 4.5 4.5 快速傅里葉反變換快速傅里葉反變換IFFTIFFT）IFFTIFFT是是IDFTIDFT的快速算法。由于的快速算法。由于DFTDFT的正變換和反變換的表達(dá)的正變換和反變換的表達(dá)式相似，因此式相似，因此IDFTIDFT也有相似的快速算法?？梢杂萌N不同也有相似的快速算法。可以用三種不同的方法來導(dǎo)出的方法來導(dǎo)出IFFTIFFT算法。算法。方法方法1 1

25、設(shè)設(shè) ， 則則有：有： ， n=0,1,N-1 n=0,1,N-1)()(kXnxDFT 10)(1)(NknkNWkXNnxn根據(jù)基2 FFT的時(shí)間抽選法的第一次分解的結(jié)果：)()()2()()()(kPWkGNkXkPWkGkX kNkN12, 1 , 0Nkn可以得到：)2()(21)()2()(21)(NkXkXWkPNkXkXkG kN12, 1 , 0Nk 圖圖 4.8 由由 X(k)、X(k+N/2) 得到得到 G(k)、P(k)n再由N/2點(diǎn)的DFT求得N/4點(diǎn)的DFT，依次類推下去，就可推到求出x(n)的各點(diǎn)，如圖4.9所示。將此流圖與圖4.4比較，相當(dāng)于整個(gè)流向反過來，此外

26、，因子WNk成為WN-k，還增加了因子1/2。但是，圖4.9的IFFT算法不能直接利用按照?qǐng)D4.4編寫的FFT算法程序，卻可以利用圖4.6的頻率抽選FFT算法的程序，只需要將X(k)作為輸入序列，因子WNk變?yōu)閃N-k，并且將最后所得的輸出序列的每個(gè)元素都除以N。n方法方法2 2n 將將DFTDFT的正變換式：的正變換式： n與其反變換式：與其反變換式： 10)()(NnnkNWnxkX10)(1)(NknkNWkXNnxn比較，很容易知道，可以利用FFT算法的程序來計(jì)算IFFT，只需要將X(k) 作為輸入序列，x(n) 則是輸出序列，另外將因子WNk 變?yōu)閃N-k， 當(dāng)然，最后還必須將輸出序

27、列的每個(gè)元素除以N。n 方法方法3n 對(duì)對(duì)DFT的反變換式取共軛，有：的反變換式取共軛，有：n 1, 1 , 0, )(1)(1)(10*10*NnWkXNWkXNnxNknkNNknkNn與DFT的正變換式比較，可知完全可以利用FFT的計(jì)算程序，只需要將X*(k) 作為輸入序列，并將最后結(jié)果取共軛，再除以N就得到x(n)。4.6 4.6 線性調(diào)頻線性調(diào)頻z z變換變換CZTCZT算法算法 如果所需要的頻率抽樣點(diǎn)并不均勻地分布在單位圓上，此時(shí)，如果所需要的頻率抽樣點(diǎn)并不均勻地分布在單位圓上，此時(shí)，常常采用常常采用Chirp zChirp z變換變換CZTCZT算法算法 4.6.1 基本原理基本

28、原理用用CZT算法可以計(jì)算下列給定點(diǎn)算法可以計(jì)算下列給定點(diǎn)zk上的上的X(z)（即在這些點(diǎn)處（即在這些點(diǎn)處的復(fù)頻譜）：的復(fù)頻譜）： (4.26)式中，式中， ， ，W0 、A0 為正實(shí)數(shù)。為正實(shí)數(shù)。1M,0,1,k AWzkk,00jeWW00jeAA n這些z平面上的抽樣點(diǎn)所沿的周線是一條螺旋線，參數(shù)W0控制周線盤旋的傾斜率。若W0大于1，則隨著k的增加，周線向內(nèi)盤旋，趨向原點(diǎn)；若W0小于1，則隨著k的增加，周線向外盤旋。若W0=1，螺旋線實(shí)際上是一段圓弧，而如果又有A0=1，則這段圓弧是單位圓的一部分。A0和0分別為第一個(gè)抽樣點(diǎn)k = 0的模和幅角，其余抽樣點(diǎn)沿螺旋周線按角度間隔0分布。

29、圖圖 4.10 z 平面上一條螺旋周線上的抽樣點(diǎn)平面上一條螺旋周線上的抽樣點(diǎn)n要計(jì)算： n式中N為序列x(n) 的長(zhǎng)度。由恒等式： n并且令 ， n就可以得到 1, 0,)()()()(1010Mk WAnxznxzXnxCZTNnnknNnnkk)(21222nknknk22)()(nnWAnxny22)(nWnh)()()(1022nkhnyWzXNnkkn改變變量，將n換為r，k換為n，則有：n n=0,1,M-1n其中， 可以看作線性時(shí)不變系統(tǒng)h(n) 對(duì)輸入信號(hào)y(n) 的響應(yīng)。 )()()()()()(22102222ngWnhnyWrnhryWzXnnNrnn)()()(nhny

30、ng 圖 4.11 CZT 算法的計(jì)算過程 4.6.2 算法要點(diǎn)算法要點(diǎn)此算法主要是要計(jì)算此算法主要是要計(jì)算y(n)與與h(n) 的線性卷積的線性卷積g(n)。由于序列。由于序列x(n) 長(zhǎng)度為長(zhǎng)度為N，因此，因此y(n)的長(zhǎng)度也為的長(zhǎng)度也為N0nN-1）。雖然）。雖然 可以是無限長(zhǎng)的，但是由于可以是無限長(zhǎng)的，但是由于z 平面上的抽樣點(diǎn)只有平面上的抽樣點(diǎn)只有M 個(gè)，個(gè)，因此因此h(n)中只有有限個(gè)點(diǎn)值是所需要的。中只有有限個(gè)點(diǎn)值是所需要的。 2/2)(nWnhn由于只需要計(jì)算M個(gè)X(zn) 之值，那么也只需要M個(gè)g(n) 之值，也即對(duì)于線性卷積n n只需要求n = 0,1, M-1時(shí)的值。如圖

31、4.12所示，由于y(r) 的范圍是r由0到N-1，因此當(dāng)求g(n) 的第一個(gè)值即n = 0時(shí)，也需要h(-r) 中r由0到N-1這N個(gè)值。 10)()()()()(Nrrnhrynhnyngn為了計(jì)算g(1)到g(M-1), 還應(yīng)該將h(-r)向右移M-1次，也就是說，h(-r)中由r =-1到r =-(M-1)這M-1個(gè)值也是所需要的。這樣，我們需要并且也只需要h(-r)中的由r =-(M-1)到r =N-1這L個(gè)值，而n L=(N-1)-(M-1)+1 = N+M-1n這也就是說，對(duì)于序列h(r)(或h(n)，我們只需要其中r或n由 -(N-1) 到M-1這一范圍的L個(gè)值。M-1r 0h

32、(r)-(N-1)(b)-(M-1)r0h(-r)N-1(c) 圖圖 4.12 序列序列 y(n) 和和 h(n) 的長(zhǎng)度和所取范圍的長(zhǎng)度和所取范圍n如果要用L點(diǎn)循環(huán)卷積來計(jì)算線性卷積g(n)，則序列y(n) 后面應(yīng)補(bǔ)M-1個(gè)0，使其長(zhǎng)度為L(zhǎng)；而用作循環(huán)卷積的有限長(zhǎng)序列h(n)為：n (4.32) 1 )(10 )()( 2/)(2/22LnMWLnhMnWnhnhLnnn然后計(jì)算L點(diǎn)循環(huán)卷積：n n其結(jié)果的L個(gè)值中，gL(0)到gL(M-1)這M個(gè)值正好與所需要的線性卷積g(0) 到g(M-1) 對(duì)應(yīng)相等，這是因?yàn)?，圖4.13(c)中的序列 與圖4.12(c中的序列h(-r)在r=-(M-1

33、)到r=N-1區(qū)間完全相同。而gL(n)的后N-1個(gè)值則是不需要的。)( )( )()(10nRrnhryngLLrL)( rh n若計(jì)算循環(huán)卷積時(shí)需要利用FFT來進(jìn)行快速計(jì)算，則要求L = 2s，s為正整數(shù)。此時(shí)，若M+N-1 2s，則要在y(n)后面補(bǔ)上M-1+K個(gè)零，使L=M+N-1+K=2s。當(dāng)然，(4.32)式所示的序列h(n) 中也要補(bǔ)K個(gè)零，這K個(gè)零應(yīng)該插在圖4.13(b)中虛線所示的地方，即在序列h(r)中r=M-1之后插入K個(gè)零實(shí)際上，可以是K個(gè)任意的數(shù)），這樣所得到的y(n)與h(n)的循環(huán)卷積仍然只取前面的M個(gè)值，即為所要求的線性卷積的M個(gè)值。 4.5.4 CZT算法的特

34、點(diǎn)算法的特點(diǎn) 1計(jì)算量計(jì)算量 按照?qǐng)D按照?qǐng)D4.11所示的計(jì)算過程，可以大致統(tǒng)計(jì)出用所示的計(jì)算過程，可以大致統(tǒng)計(jì)出用CZT算法算法計(jì)算計(jì)算N點(diǎn)長(zhǎng)的時(shí)域序列點(diǎn)長(zhǎng)的時(shí)域序列x(n) 在在z平面的一段螺旋狀周線上的平面的一段螺旋狀周線上的M個(gè)點(diǎn)處的復(fù)頻譜所需要的乘法次數(shù)。個(gè)點(diǎn)處的復(fù)頻譜所需要的乘法次數(shù)。(1) 計(jì)算 需要N 次復(fù)乘。(2) y(n)的L點(diǎn)FFT運(yùn)算需要 次復(fù)乘。(3) 若保持M和N不變，則DFT H(k) 可以事先算好存起來，而不需要現(xiàn)場(chǎng)運(yùn)算操作。(4) 求G(k) = Y(k)H(k) 需要L次復(fù)乘。2/2)()(/)()(nnnWAnxnhAnxny)2/(log)2/(2LL(5

35、) 對(duì)G(k) 進(jìn)行L點(diǎn)IFFT運(yùn)算以得到g(n)，需求 次復(fù)乘。(6) 將g(n) 乘上 求得M點(diǎn)輸出X(zn)，需要M次復(fù)乘。 這些復(fù)數(shù)乘法中主要是兩次L點(diǎn)FFT(IFFT)的運(yùn)算量，因此，CZT算法的運(yùn)算量大致與 成正比，這里L(fēng)N+M-1。)2/(log)2/(2LL2/2nW2log22LLn如果用z變換直接計(jì)算復(fù)頻譜，即：n n則需要MN次復(fù)數(shù)乘法。實(shí)際上，當(dāng)M、N較小時(shí)，直接用z變換式來計(jì)算會(huì)比較方便；乘積MN越大，CZT算法的運(yùn)算量的減少越顯著。1M,0,1,k znxzXNnnkk10)()( 2 2與與FFTFFT算法比較算法比較與標(biāo)準(zhǔn)的與標(biāo)準(zhǔn)的FFTFFT算法相比較，算法相

36、比較，CZTCZT算法有以下特點(diǎn)：算法有以下特點(diǎn)：(1) (1) 輸入和輸出序列長(zhǎng)度不需要相等，即不需要輸入和輸出序列長(zhǎng)度不需要相等，即不需要M = NM = N。(2) N(2) N與與M M都不需要是都不需要是2 2的正整數(shù)冪。的正整數(shù)冪。(3) zk (3) zk 間的間隔可以任意選擇，這樣就可得到任意的分辨間的間隔可以任意選擇，這樣就可得到任意的分辨率；而且當(dāng)所需間隔不同時(shí)，可以分段計(jì)算。率；而且當(dāng)所需間隔不同時(shí)，可以分段計(jì)算。 (4) zk 點(diǎn)所沿周線不必須是圓弧。(5) 起始點(diǎn)的選擇可以是任意的，因此便于從任意的頻率上開始對(duì)輸入數(shù)據(jù)進(jìn)行頻譜分析。(6) 當(dāng)(4.26)式中A=1，W

37、=e-j2/N ，并且N=M時(shí)，X(zk)即為x(n)的DFT，因此利用CZT算法也能快速計(jì)算x(n)的DFT，而且不要求N為2的正整數(shù)冪。 4.7 4.7 實(shí)序列的實(shí)序列的FFTFFT的高效算法的高效算法 上面討論的上面討論的FFTFFT算法是將算法是將x(n)x(n)作為復(fù)數(shù)來對(duì)待的，如果作為復(fù)數(shù)來對(duì)待的，如果x(n)x(n)是實(shí)序列，計(jì)算量還能夠進(jìn)一步減少。是實(shí)序列，計(jì)算量還能夠進(jìn)一步減少。4.7.1 兩個(gè)長(zhǎng)度相同的實(shí)序列兩個(gè)長(zhǎng)度相同的實(shí)序列 可以將兩個(gè)長(zhǎng)度相同的實(shí)序列組合成一個(gè)復(fù)序列來進(jìn)行可以將兩個(gè)長(zhǎng)度相同的實(shí)序列組合成一個(gè)復(fù)序列來進(jìn)行FFT運(yùn)算，從而一次完成這兩個(gè)實(shí)序列的運(yùn)算，從而一

38、次完成這兩個(gè)實(shí)序列的FFT，減少了總，減少了總的計(jì)算量。的計(jì)算量。n 設(shè)p(n) 和g(n) 是兩個(gè)長(zhǎng)度均為N的實(shí)序列，并設(shè)：n n又設(shè) ， ， n則由DFT的線性有： (4.36)()()(njgnpny)()(kPnpDFT )()(kGngDFT )()(kYnyDFT )()()(kjGkPkYnP(k) 和G(k) 這兩個(gè)復(fù)序列的實(shí)部都是周期性的偶序列，而其虛部都是周期性的奇序列。n 對(duì)復(fù)序列Y(k) 又有 (4.37)n這里下標(biāo) r、i 分別表示實(shí)部和虛部，Y(k) 與其實(shí)部、虛部的長(zhǎng)度都為 N?，F(xiàn)將 (4.37) 式中各序列作周期延拓，有 (4.38)()()(kjYkYkYir

39、)()()(kYjkYkYirn由周期性有：n (4.39)n (4.40)()(),()(kNYkY kNYkYrrrr)()(),()(kNYkY kNYkYiiiin現(xiàn)在將序列 與 作如下分解：n (4.41)n (4.42)(kYr)(kYi)()(21)()(21)(kNYkYkNYkYkYrrrrr)()(21)()(21)(kNYkYkNYkYkYiiiiin 由4.39式和4.40式，容易證明在(4.41)和(4.42)這兩個(gè)式子中，前一項(xiàng)都是偶序列，而后一項(xiàng)都是奇序列。n 將4.41式和4.42式代入4.38式，并將各項(xiàng)進(jìn)行重新組合，得到：)( )( )()(21)()(21

40、)()(21)()(21)(kGjkPkNYkYjkNYkYjkNYkYjkNYkYkYrriiiirrn令0kN-1，則上式為： n這里的P(k) 和G(k) 的實(shí)部都是周期性偶序列，而它們的虛部都是周期性奇序列，此情況與4.36式中的復(fù)序列P(k)和G(k)的情況相同。因此有：)()()( )( )(kjGkPkjGkPkY 上兩式中0kN-1。)()(21)()(21)( )(NiiNrrkNYkYjkNYkYkPkP)()(21)()(21)( )(NrrNiikNYkYjkNYkYkGkG 4.7.2 4.7.2 一個(gè)一個(gè)2N2N點(diǎn)的實(shí)序列點(diǎn)的實(shí)序列 將一個(gè)將一個(gè)2N2N點(diǎn)的實(shí)序列點(diǎn)

41、的實(shí)序列x(n)x(n)按偶數(shù)點(diǎn)和奇數(shù)點(diǎn)分組形成兩按偶數(shù)點(diǎn)和奇數(shù)點(diǎn)分組形成兩個(gè)個(gè)N N點(diǎn)實(shí)序列：點(diǎn)實(shí)序列： 1, 1 , 0 ) 12()()2()(Nnnxngnxnpn則有：n k=0,1,N-1 (4.48)()()()()()( k 2k 2kGWkPNkXkGWkPkXNNn其中：n n實(shí)序列p(n) 和g(n) 的DFT P(k) 和G(k) 可以采用4.7.1節(jié)所說的方法作一次N點(diǎn)復(fù)序列的FFT而同時(shí)得到，然后再按4.48式進(jìn)行組合便得到了2N點(diǎn)實(shí)序列x(n) 的DFT。10)()(NnnkNWnpkP1, 1 , 0)()(10NkWngkGNnnkN4.8 Matlab方法方法 4.8.1 利用利用Matlab計(jì)算計(jì)算FFT在第在第3章中介紹用章中介紹用Matlab方法計(jì)算信號(hào)的方法計(jì)算信號(hào)的DFT時(shí)，提到了函時(shí)，提到了函數(shù)數(shù)fft(x,N) 和和ifft(x.N)。對(duì)于這兩個(gè)函數(shù)，如果。對(duì)于這兩個(gè)函數(shù)，如果N為為2的正的正整數(shù)冪，則可以得到本章中介紹的基整數(shù)冪，則可以得到本章中介紹的基2 FFT快速