1、2.3 快速傅里葉變換 （FFT） 快速傅里葉變換（FFT）是計算DFT的一種快速有效方法。從前面的討論中看到，有限長序列在數字技術中占有很重要的地位。有限長序列的一個重要特點是其頻域也可以離散化，即離散傅里葉變換（DFT）。 第1頁/共91頁 雖然頻譜分析和DFT運算很重要，但在很長一段時間里，由于DFT運算復雜，并沒有得到真正的運用，而頻譜分析仍大多采用模擬信號濾波的方法解決，直到1965年首次提出DFT運算的一種快速算法以后，情況才發生了根本變化，人們開始認識到DFT運算的一些內在規律，從而很快地發展和完善了一套高速有效的運算方法快速付里變換（FFT）算法。FFT的出現，使DFT的運算大

2、大簡化，運算時間縮短一二個數量級，使DFT的運算在實際中得到廣泛應用。第2頁/共91頁1、DFT運算的特點： 首先分析有限長序列 x（n）進行一次DFT運算所需的運算量。 一般，x（n）和wnkN都是復數，因此，每計算一個X（k）值，要進行N次復數相乘，和N-1次復數相加，X（k）一共有N個點，故完成全部DFT運算，需要N2次復數相乘和N（N-1）次復數相加，在這些運算中，乘法比加法復雜，需要的運算時間多，尤其是復數相乘，每個復數相乘包括4個實數相乘和2個實數相加，例 又每個復數相加包括2個實數相加，所以，每計算一個 X（k）要進行4N次實數相乘和2N+2（N-1）=2（2N-1）次實數相加，

3、因此，整個DFT運算需要4N2實數相乘和2N（2N-1）次實數相加。 101, 1 , 0)()()(NnnkNNkwnxnxDFTkX)()()()()(10nkNemnkNmeNnnkNmmnkNeewRnxIwInxRjwInxIwRnxRkX第3頁/共91頁 從上面的分析看到，在DFT計算中，不論是乘法和加法，運算量均與N2成正比。因此，N較大時，運算量十分可觀。例，計算N=10點的DFT，需要100次復數相乘，而N=1024點時，需要1048576（一百多萬）次復數乘法，如果要求實時處理，則要求有很高的計算速度才能完成上述計算量。 反變換IDFT與DFT的運算結構相同，只是多乘一個常

4、數1/N，所以二者的計算量相同。第4頁/共91頁FFT算法的基本思想：考察DFT與IDFT的運算發現，利用以下兩個特性可減少運算量：1）系數 是一個周期函數，它的周期性和對稱性可用來改進運算，提高計算效率。 例 又如 因此 利用這些周期性和對稱性，使DFT運算中有些項可合并；2）利用 的周期性和對稱性，把長度為N點的大點數的DFT運算依次分解為若干個小點數的DFT。因為DFT的計算量正比于N2，N小，計算量也就小。FFT算法正是基于這樣的基本思想發展起來的。它有多種形式，但基本上可分為兩類：時間抽取法和頻率抽取法。nkNjnkNew2nkNnNkNkNnNwww)()(, 12/NNwkNNk

5、Nww)2/(nkNw第5頁/共91頁2、按時間抽取的FFT（N點DFT運算的分解） 先從一個特殊情況開始，假定N是2的整數次方， N=2M，M：正整數 首先將序列x（n）分解為兩組，一組為偶數項，一組為奇數項， r=0,1，N/2-1 )() 12()()2(21rxrxrxrx第6頁/共91頁2011)()(NnNnnkNnkNwnxwnx偶數奇數10)()()(NnnkNwnxnxDFTkx：將DFT運算也相應分為兩組12/0)12(12/02) 12()2(NrkrNNrrkNwrxwrx12/0212/02) 12()2(NrrkNNrkNrkNwrxwwrx第7頁/共91頁因為 故

6、 其中nNnNjnNjnNweew222222 kHWkGWrxWWrxkXkNrkNNrkNNrrkN/)()( /)(NrrkNWrxkG rkNNrWrxkH/)(第8頁/共91頁 注意到，H（k），G（k）有N/2個點，即k=0，1， N/2-1，還必須應用系數 wkN 的周期性和對稱性 表示 X（k）的 N/2 N-1點： 由 得：rkNkNrNww22/2 12, 1 , 0,)2(NkkHWkGNkXkNkNNkNWW)2( 可見，一個N點的DFT被分解為兩個N/2點的DFT，這兩個N/2點的DFT再合成為一個N點DFT.第9頁/共91頁 NkkHWkGkXkN, 12, 1 ,

7、 0,)2(NkkHWkGNkXkN 依此類推，G(k)和H(k)可以繼續分下去，這種按時間抽取算法是在輸入序列分成越來越小的子序列上執行DFT運算，最后再合成為點的DFT。 第10頁/共91頁蝶形信號流圖將G(k)和H(k) 合成X(k)運算可歸結為： bWabWaWa+bWa-bW-W1a+bWa-bWWabab蝶 形 運 算的簡化(a)(b)第11頁/共91頁 圖 (a)為實現這一運算的一般方法，它需要兩次乘法、兩次加減法。考慮到-bW和bW兩個乘法僅相差一負號，可將圖 (a)簡化成圖2.7(b)，此時僅需一次乘法、兩次加減法。圖 (b)的運算結構像一蝴蝶通常稱作蝶形運算結構簡稱蝶形結，

8、采用這種表示法，就可以將以上所討論的分解過程用流圖表示。 第12頁/共91頁N=23=8 的例子。 第13頁/共91頁 按照這個辦法，繼續把N/2用除，由于N=2M，仍然是偶數，可以被整除，因此可以對兩個N/2點的DFT再分別作進一步的分解。即對G(k)和H(k)的計算，又可以分別通過計算兩個長度為N/4=2點的DFT，進一步節省計算量，見圖。這樣，一個點的DFT就可以分解為四個點的DFT。 第14頁/共91頁第15頁/共91頁最后剩下的是2點DFT，它可以用一個蝶形結表示：這樣，一個8點的完整的按時間抽取運算的流圖由于這種方法每一步分解都是按輸入時間序列是屬于偶數還是奇數來抽取的，所以稱為“

9、按時間抽取法”或“時間抽取法”。) 1 ()0() 1 ()0() 1 () 1 ()0() 1 ()0()0(1202xWxxWxXxWxxWxXoNoN第16頁/共91頁按時間抽取的按時間抽取的8點點FFT第17頁/共91頁 時間抽取法FFT的運算特點：（1）蝶形運算（2）原位計算（3）序數重排（4）蝶形類型隨迭代次數成倍增加第18頁/共91頁（1）蝶形運算 對于N=2M，總是可以通過M次分解最后成為2點的DFT運算。這樣構成從x(n)到X(k)的M級運算過程。從上面的流圖可看到，每一級運算都由N/2個蝶形運算構成。因此每一級運算都需要N/2次復乘和N次復加，這樣，經過時間抽取后M級運算總

10、共需要的運算： 復乘 復加 而直接運算時則與N2 成正比。例N=2048，N2=4194304，(N/2)log2N=11264，N2 / (N/2)log2N=392.4。FFT顯然要比直接法快得多。NNMN2log22NNMN2log第19頁/共91頁（2）原位計算 當數據輸入到存儲器中以后，每一級運算的結果仍然儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，這叫原位計算。 每一級運算均可在原位進行，這種原位運算結構可節省存儲單元，降低設備成本，還可節省尋址的時間。(考填空題)第20頁/共91頁（3）序數重排 對按時間抽取FFT的原位運算結構，當運算完畢時，正好順序存放著 X（0）

11、，X（1），X（2），X（7），因此可直接按順序輸出，但這種原位運算的輸入 x（n）卻不能按這種自然順序存入存儲單元中，而是按x(0)，x(4)，x(2)，x(6)，x(7)的順序存入存儲單元，這種順序看起來相當雜亂，然而它也是有規律的。當用二進制表示這個順序時，它正好是“碼位倒置”的順序。例如，原來的自然順序應是 x(1)的地方，現在放著 x(4)，用二進制碼表示這一規律時， 則是在 x（0 0 1）處放著 x（1 0 0）， x（0 1 1）處放著 x（1 1 0）。第21頁/共91頁 表 碼位倒置順序自然順序二進碼表示碼位倒置碼位倒置順序0000000010011004201001023

12、011110641000011510110156110010371111117 在實際運算中，一般直接將輸入數據 x（n）按碼位倒置的順序排好輸入很不方便，總是先按自然順序輸入存儲單元，然后再通過變址運算將自然順序的存儲轉換成碼位倒置順序的存儲，然后進行FFT的原位計算。目前有許多通用DSP芯片支持這種碼位倒置的尋址功能。第22頁/共91頁 第一次分偶、奇，根據最低位n0的0、1狀態來分，若n0=0，則為偶序列；n0=1則為奇序列，得到兩組序列：000 010 100 110 001 011 101 111 第二次對這兩個偶、奇序列再分一次偶、奇序列，這就要根據n1的、狀態。若n1=0，則為偶

13、序列；n1=1則為奇序列，得到四組序列：000 100 010 110 001 101 011 111 同理，再根據n2的、狀態來分偶、奇序列，直到不能再分偶、奇時為止。對于N=8, n2已是最高位，最后一次分得結果為 000 100 010 110 001 101 011 111第23頁/共91頁（4）蝶形類型隨迭代次數成倍增加 觀察8點FFT的三次迭代運算: 第一級迭代，有一種類型的蝶形運算系數W08，兩個數據點間隔為1 第二級迭代，有二種類型的蝶形運算系數W08、W28，參加 運算的兩個數據點間隔為2。 第三級迭代，有四類蝶形運算系數W08、W18、W28、W38，參加運算的兩個數據點間

14、隔為4。 結論：每迭代一次，蝶形類型增加一倍，數據點間隔也增大一倍。 每一級的取數間隔和蝶形類型種類均為2i-1，i=1,2,M。 第24頁/共91頁3、按頻率抽取的FFT（按輸出X（k）在頻域的順序上屬于偶數還是奇數分解為兩組） 對于N=2M情況下的另外一種普遍使用的FFT結構是頻率抽取法。對于頻率抽取法，輸入序列不是按偶數奇數，而是按前后對半分開，這樣便將N點DFT寫成前后兩部分： 12/012/)()()(NnNNnnkNnkNWnxWnxkX12/012/0)2()2()(NnNnkNnNnkNWNnxWnx12/0)2/()2/()(NnnkNkNNWNnxWnx第25頁/共91頁又

15、把X（k）進一步分解為偶數組和奇數組： = X(2r+1)= =奇數為偶數1, 1) 1(, 1)2/(2/kWWkkNNNN12/0)2/() 1()()(NnnkNkWNnxnxkX12/02)2/()()2(NnnrNWNnxnxrX12/02/)2/()(NnnrNWNnxnx12/0)12()2/()(NnrnNWNnxnx12/02/)2/()(NnnrNnNWWNnxnx第26頁/共91頁令 a(n)=x(n)+x(n+N/2) b(n)=x(n)-x(n+N/2wnN這兩個序列都是N/2點的序列，將其代入上兩式，得 這正是兩個N/2點的DFT運算，即將一個N點的DFT分解為兩個

16、N/2點的DFT，上式的運算關系可用下圖表示. 12/02/)()2(NnnrNWnarXnrNNnWnbrX2/12/02)() 12(x1(n)+ x2(n)WNnx1(n)x2(n)x1(n)-x2(n) WNn-1第27頁/共91頁以N=8的頻率抽取為例x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)-1-1-1-1N/2點DFTa(0)a(1)a(2)a(3)N/2點DFTb(0)b(1)b(2)b(3)WN1WN2WN3X(0)X(2)X(4)X(6)X(1)X(3)X(5)X(7)按頻率抽取將8點DFT分解成兩個4點DFT第28頁/共91頁按頻率抽取將8點DFT分

17、解成四個2點DFT第29頁/共91頁 與時間抽取法一樣，由于N=2M，N/2仍是一個偶數,這樣，一個N=2M點的DFT通過M次分解后，最后只剩下全部是2點的DFT，2點DFT實際上只有加減運算。但為了比較，也為了統一運算的結構，仍用一個系數為W0N 的蝶形運算來表示。頻率抽取法的流圖是時域抽取法流圖的左右翻轉。下圖是N=8的頻率抽取法FFT流圖。 N=8的按頻率抽取FFT運算流圖第30頁/共91頁頻率抽取法FFT的運算特點：（1）蝶形運算 對于任何一個2的整數冪N=2M，總是可以通過M次分解最后完全成為2點的DFT運算。這樣的M次分解，就構成從x(n)到X(k)的M級運算過程。將頻率抽取法的流

18、圖反轉，并將輸入變輸出，輸出變輸入，得到的便是時間抽取法流圖（反映了時域，頻域的對稱法）。頻率抽取法也共有M級運算 （N=2M），其運算量與時間抽取法相同。第31頁/共91頁（2）原位計算 類似于時間抽取法，當數據輸入到存儲器中以后，每一級運算的結果仍然儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，所以頻域抽取法也可進行原位算。第32頁/共91頁（3）序數重排 它的輸入正好是自然順序。但它的輸出卻是碼位倒置的順序。因此運算完畢后，要通過變址運算將碼位倒置的順序轉換為自然順序，然后輸出，變址方法同時間抽取法。第33頁/共91頁（4）蝶形類型隨迭代次數成倍減少(與時間抽取法相反) 第一級

19、迭代中有N/2種蝶形運算系數，參加蝶形運算的兩個數據相隔N/2，隨后每次迭代，蝶形類形比前一級減少一倍，間距也減少一倍，最后一級迭代，蝶形類形只有一種W0N，數據間隔為1。由這幾點規律可以看出，頻率抽取法與時間抽到法是兩種等價的FFT運算。第34頁/共91頁5. N為組合數的FFT（任意基數的FFT算法） 以上討論的都是以2為基數的FFT算法，即N=2M，這種情況實際上使用得最多。 優點：程序簡單，效率高，使用方便。 實際應用時，有限長序列的長度N很大程度上由人為因素確定，因此多數場合可取 N=2M，從而直接使用以 2 為基數的FFT算法。 如N不能人為確定 , N的數值也不是以2為基數的整數

20、次方,處理方法有兩種: 補零: 將x(n)補零 , 使N=2M. 例如N=30,補上x(30)=x(31)=0兩點,使N=32=25,這樣可直接采用以2為基數M=5的FFT程序。 有限長度序列補零后并不影響其頻譜 X(ejw) ，只是頻譜的采樣點數增加了,上例中由30點增加到32點，所以在許多場合這種處理是可接受的。第35頁/共91頁 如要求準確的N點DFT值,可采用任意數為基數的 FFT 算法 , 其 計算效率低于以2為基數FFT算法。 如 N 為復合數，可分解為兩個整數p與q的乘積,像前面以2為基數時一樣,FFT的基本思想是將DFT的運算盡量分小,因此,在N=pq情況下,也希望將N點的DF

21、T分解為p個q點DFT或q個p點DFT,以減少計算量。步驟：kPkknQnn10,kn 分別為0, 1,，Q-1; 01,kn分別為0, 1,，P-1。 第36頁/共91頁 N點DFT可以重新寫成為 100101)(,)()(NnknNWnxkkXkPkXkX1010)(01010101QnPnnQnkPkNWnQnxQnPnnkNPnkNQnkNQnPnnkNPnkNQnkNPQnkNWWWnnxWWWWnnxkkX,第37頁/共91頁考慮到 再令1010nkPQnkNWW令 1010010100100110,QnnkQnkNPnnkPWWWnnxkkXPnnkPWnnxnkX,01000,

22、0011001nkQnkNQnWWnkXkkX00001001,nkNWnkXnkXQnnkQWnkXkkX,第38頁/共91頁以P=3，Q=4, N=12為例 (1) 先將 x(n)通過 x(n1Q+n0)改寫成 x(n1,n0)。因為 Q=4, n1=0,1,2, n0=0,1,2,3,故輸入是按自然順序的，即x(0,0)=x(0) x(0,1)=x(1) x(0,2)=x(2) x(0,3)=x(3)x(1,0)=x(4) x(1,1)=x(5) x(1,2)=x(6) x(1,3)=x(7)x(2,0)=x(8) x(2,1)=x(9) x(2,2)=x(10) x(2,3)=x(11

23、)第39頁/共91頁(2)(2) 求個點的DFT (3)X1(k0,n0)乘以得到X1(k0,n0)。(4)求P個點的DFT，參變量是k020301001110,nnkWnnxnkX304001102001,nnkWnkXkkX(5)將X2(k0,k1)通過X(k0+k1P)恢復為X(k)第40頁/共91頁N=12為組合數時的FFT第41頁/共91頁（）求個P點DFT需要P2次復數乘法和P(P-1)次復數加 法；（）乘個W因子需要次復數乘法；（）求P個點DFT需要PQ2 次復數乘法和P（-1)次復數加法。總的復數乘法量: QP2+N+PQ2=N(P+Q+1)；總的復數加法量: QP(P-1)+

24、PQ(Q-1)=N(P+Q-2) 第42頁/共91頁 上述分解原則可推廣至任意基數的更加復雜的情況。例如, 如果N可分解為m個質數因子p1,p2,，pm，即 N=p1p2p3pm則第一步：可把N先分解為兩個因子N=p1q1，其中q1=p2p3pm ，并用上述討論的方法將DFT分解為p1個q1點DFT；第二步：將q1分解為q1=p2q2，q2=p3p4pm，然后將每個q1 點DFT再分解為p2個q2點DFT；： 依此類推，：通過m次分解，一直分到最少點數的DFT運算，從而獲得最高的運算效率。其運算量近似為N(p1+p2+pm)次復數乘法和復數加法。但計算效率的提高是要以編程的復雜性為代價的，一般

25、較少應用。 p1=p2 = =pm =2，為基2 FFT 算法。第43頁/共91頁當組合數 N=P1P2P3Pm中所有的Pi均為4時，就是基四FFT算法 以N=43為例,第一級運算的一般形式為： ), 3(), 2(), 1 (), 0(111111111111), 3(), 2(), 1 (), 0(), 3(), 2(), 1 (), 0(212121210101010194643404644424043424140404040404011011011011nnxnnxnnxnnxjjjjnnxnnxnnxnnxWWWWWWWWWWWWWWWWnnXnnXnnXnnX30,),(),(04

26、300120101022kWnnnxnnkXknn第44頁/共91頁6、Chirp-z變換 采用FFT可以算出全部N點DFT值，即z變換X（z）在z平面單位圓上的等間隔取樣值，但要求 N 為復合數。問題的提出： 不需要計算整個單位圓上z變換的取樣，如對于窄帶信號，只需要對信號所在的一段頻帶進行分析，這時,希望頻譜的采樣集中在這一頻帶內，以獲得較高的分辨率，而頻帶以外的部分可不考慮。 對其它圍線上的z變換取樣感興趣，例如語音信號處理中，需要知道z變換的極點所在頻率，如極點位置離單位圓較遠，則其單位圓上的頻譜就很平滑，如果采樣不是沿單位圓而是沿一條接近這些極點的弧線進行，則在極點所在頻率上的將出現

27、明顯的尖峰，由此可較準確地測定極點頻率。 要求能有效地計算當N是素數時序列的DFT。 第45頁/共91頁算法原理： 螺旋線采樣是一種適合于這種需要的變換，且可以采用FFT來快速計算，這種變換也稱作Chirp-z變換。 已知x（n），0nN-1令zk=AW-k ，k=0，M-1 ，M ：采樣點數，A、W：任意復數其中：A0表示起始取樣點的半徑長度，通常A010表示起始取樣點z0的相角0表兩相鄰點之間的等分角W0螺旋線的伸展率，W01則線外伸，W01則線內縮（反時針），W0=1則表示半徑為A0的一段圓弧，若A0=1，這段圓弧則是單位圓的一部分。00jojoeWWeAA第46頁/共91頁圖 螺旋線采

28、樣 第47頁/共91頁計算z變換在采樣點 zk 的值 k=0,1, ,M-1 顯然，按照以上公式計算出全部M點采樣值需要 NM 次復乘和（N-1）M次復加，當N及M較大時，計算量迅速增加，以上運算可轉換為卷積形式，從而可采用FFT進行，這樣可大大提高計算速度。利用zk的表示式代入10)()(NnnkkznxzX nk可以用以下表示式來替換 10)()(NnnknkWAnxzX)(21222nknknk第48頁/共91頁2222)(1022)()(nkNnnnkkWWAnxWzX則222)(,)()(2nnnWnhWAnxng定義：1, 1 , 0)()()()()(210222MkkhkgWn

29、khngWzXkNnkk則第49頁/共91頁 上式說明，如對信號x（n）先進行一次加權處理，加權系數為 ，然后通過一個單位脈沖響應為h（n）的線性系統，最后，對該系統的前M點輸出再作一次 的加權 ，就可得到全部M點螺旋線采樣值。 系統的單位脈沖響應 與頻率隨時間成線性增加的線性調頻信號相似，因此稱為Chirp -z變換。 22nnWA22kW22)(nWnhx22)(nWnhxx(n)g(n)X(zk)22nnWA22kW第50頁/共91頁算法實現： 由于輸入信號 g（n）是有限長的，長為N，但序列 是無限長的，而計算 0M-1 點卷積 g（k）*h（k）所需要的 h（n）是取值在n=-（N-

30、1）M-1那一部分的值，因此，可認為h（n）是一個有限長序列，長為L=N+M-1。 所以，Chirp -z變換為兩個有限長序列的線性卷積g（k）*h（k），可用圓圈卷積通過FFT來實現。 h（n）的主值序列 可由h（n）作周期延拓后取0nL-1部分值獲得，將 與g（n）作圓周卷積后，其輸出的前M個值就是Chirp -z變換的M個值。這個圓周卷積的過程可在頻域上通過FFT實現。22)(nWnh)(nh)(nh第51頁/共91頁Chirp-z變換的計算步驟： (1) 求h（n）的主值序列 (2 ) 用FFT求 的傅里葉變換: H(k)=FFT , L點 (3) 對x（n）加權并補零1110)(2/

31、)(2/22LnNLWMnWnhnLn)(nh)(nh（4）G（k）=FFT g（n） , L點（5）Y（k）=G（k）H（k）, L點（6）y（n）=IFFTY（k） , L點（7） , 0kM-11010)()(2/2LnNNnWAnxngnn)()(22kYWzXkk第52頁/共91頁利用FFT計算Chirp-z變換第53頁/共91頁計計算量估算：乘 乘法數估計 （1）（2）兩步可以事先計算，不必實時計算。（3）（7）兩步兩次加權共計N+M次復乘,形成Y（k）需L次復乘，一個FFT與一個IFFT需Llog2L次乘，所以總乘法數：L+N+M+Llog2L，直接計算乘法數， NM N及M較大

32、時，用FFT實現Chirp-Z變換，速度上有很大的改進。 第54頁/共91頁Chirp-z變換的特點：1）輸入序列的長度N 與 輸出序列的長度 M 不需要相等；2）N 及 M 不必是高度復合數，二者均可為素數；3）相鄰采樣點 zk 之間的角間隔 0 是任意的，即頻率分辨率是任意的；4）圍線是任意的，不必是 Z 平面上的圓；5）起始點 z0 可任意選定，即可從任意頻率上開始對輸入數據進行窄帶高分辨率分析；6）若 A=1 , M=N , , 可用 Chirp-z變換變換計算DFT（即使 N 為素數）。NjeW2第55頁/共91頁2.4 FFT應用中的幾個問題1)IDFT的運算方法 所討論的FFT算

33、法可用于IDFT運算簡稱為IFFT 比較IDFT的定義式： IDFT： DFT：X(k)=DFTx(n) =10)(1)()(NknkNWkXNkXIDFTnxnkNNnWnx10)(第56頁/共91頁IDFT與DFT的差別： 1）把DFT中的每一個系數 改為 ， 2）再乘以常數 1/N ， 則以上所討論的時間抽取或頻率抽取的FFT運算均可直接用于IDFT運算，當然，蝶形中的系數 應改為 。nkNWnkNWnkNWnkNW第57頁/共91頁10)(1)(NknkNWkXNnx)(1)(1)(10kXDFTNWkXNnxNknkN)(kX b）第二種方法，完全不需要改動FFT程序，而是直接利用它

34、作 IFFT。 考慮到 故 IFFT計算分三步： 將X（k）取共軛（虛部乘以-1） 對 直接作FFT 對FFT的結果取共軛并乘以1/N，得x（n）。第58頁/共91頁2.4 FFT應用中的幾個問題2）實數序列的FFT以上討論的FFT算法都是復數運算,包括序列x(n)也認為是復數,但大多數場合,信號是實數序列,任何實數都可看成虛部為零的復數,例如,求某實信號x(n)的復譜,可認為是將實信號加上數值為零的虛部變成復信號(x(n)+j0),再用FFT求其離散傅里葉變換。這種作法很不經濟,因為把實序列變成復序列,存儲器要增加一倍,且計算機運行時,即使虛部為零,也要進行涉及虛部的運算,浪費了運算量。合理

35、的解決方法是利用復數據FFT對實數據進行有效計算,下面介紹兩種方法。第59頁/共91頁（1）用一個N點FFT同時計算兩個N點實序列 的DFT 設 x (n)、y (n)是彼此獨立的兩個N點實序列,且 X (k)=DFTx (n) ， Y (k)=DFTy(n) 則X (k)、Y(k)可通過一次FFT運算同時獲得。 首先將x (n)、y(n)分別當作一復序列的實部及虛部,令 g(n)=x (n)+jy(n)第60頁/共91頁 通過FFT運算可獲得g(n)的DFT值 ) 利用離散傅里葉變換的共軛對稱性 通過g(n)的FFT運算結果G(k),由上式可得到X (k) 的值。 kjGkGkjYkjXkY

36、kXkYkjYkjXkXkjYkXkGirriiririr kNGkGjkNGkGkNGkGkXngiirr212121ReDFT*第61頁/共91頁通過g(n)的FFT運算結果G(k),由上式也可得到X (k) 的值。 kNGkGjkNGkGkNGkGkjYngjiirr212121ImDFT* kNGkGjkNGkGkYrrii2121作一次點復序列的FFT，再通過加、減法運算就可以將X(k)與Y(k)分離出來。顯然，這將使運算效率提高一倍。第62頁/共91頁 （2）用一個N點的FFT運算獲得一個2N點實 序列的DFT 設x(n)是2N點的實序列,現人為地將x(n)分為偶數組x1(n) 和

37、奇數組x2(n) x1(n)=x(2n) n=0,1,N-1 x2(n)=x(2n+1) n=0,1,N-1 然后將x1(n)及x2(n)組成一個復序列： y(n)=x1(n)+jx2(n) 第63頁/共91頁通過N點FFT運算可得到： Y(k)=X1(k)+jX2(k) ，N點根據前面的討論，得到 )()(2)()()(21)(*2*1kNYkYjkXkNYkYkX第64頁/共91頁 為求 2N 點 x（n）所對應 X（k），需求出 X（k）與 X1（k）、X2（k）的關系 ： 而 10)12(21022120) 12()2()()(NnknNNnnkNNnnkWnxWnxWnxkX1021

38、0) 12()2(NnnkNkNNnnkNWnxWWnx 1011022101011) 12()()()2()()(NnnkNNnnkNNnNnnkNnkNWnxWnxkXWnxWnxkX第65頁/共91頁所以1）由x1(n)及x2(n)組成復序列，經FFT運算求得 Y(k)，2）利用共軛對稱性求出 X1(k)、X2(k)，3）最后利用上式求出 X（k）, 達到用一個N點的FFT計算一個2N點的實序列的DFT的目的。 X(k)=X1(k)+W2Nk X2(k)第66頁/共91頁3)線性卷積的FFT算法 線性卷積是求離散系統響應的主要方法之一,許多重要應用都建立在這一理論基礎上,如卷積濾波等。

39、以前曾討論了用圓周卷積計算線性卷積的方法歸納如下: 將長為N2的序列x(n)延長到L,補L-N2個零， 將長為N1的序列h(n)延長到L,補L-N1個零， 如果LN1+N2-1,則圓周卷積與線性卷積相等,此時,可用FFT計算線性卷積，方法如下: 第67頁/共91頁 a.計算X(k)=FFTx(n) b. 求H(k)=FFTh(n) c. 求Y(k)=H(k)X(k) k=0L-1 d. 求y(n)=IFFTY(k) n=0L-1 可見，只要進行二次FFT,一次IFFT就可完成線性卷積計算。計算表明,L32時,上述計算線性卷積的方法比直接計算線卷積有明顯的優越性,因此,也稱上述循環卷積方法為快速

40、卷積法。第68頁/共91頁 上述結論適用于 x(n)、h(n) 兩序列長度比較接近或相等的情況,如果x(n)、h(n) 長度相差較多,例如, h(n) 為某濾波器的單位脈沖響應,長度有限,用來處理一個很長的輸入信號 x(n),或者處理一個連續不斷的信號，按上述方法, h(n) 要補許多零再進行計算,計算量有很大的浪費，或者根本不能實現。 為了保持快速卷積法的優越性,可將 x(n) 分為許多段,每段的長度與 h(n) 接近 , 處理方法有兩種：第69頁/共91頁(1) 重疊相加法由分段卷積的各段相加構成總的卷積輸出h(n)x(n)第70頁/共91頁 假 定 xi(n) 表示 x（n）序列的第i段

41、 ： 則 輸入序列可表為： 于是輸出可分解為： 其中 01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnhnxny)()(*)()(*)()()(*)()(nhnxnyii第71頁/共91頁 1）先對 h（n）及 xi（n）補零，補到具有N點長度， N=N1+N2-1。 一般選 N=2M。 2）用基2 FFT計算 yi（n）=xi（n）*h（n）。 3）重疊部分相加構成最后的輸出序列。 由于 yi（n）的長度為N，而xi（n）的長度為N2，因此相鄰兩段 yi（n）序列必然有N-N2=N1-1點發生重疊。 iinyny)()(第72頁/共91頁計算步驟：a.

42、事先準備好濾波器參數 H（k）=DFTh(n)，N點b. 用N點FFT計算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用N點IFFT求 yi（n）=IDFTYi(k)e. 將重疊部分相加第73頁/共91頁第74頁/共91頁(2)重疊保留 這種方法和第一種方法稍有不同，即將上面分段序列中補零的部分不是補零，而是保留原來的輸入序列值，這時，如利用DFT實現h(n)和xi(n)的循環卷積，則每段卷積結果中有N1-1個點不等于線性卷積值需舍去。 重疊保留法與重疊相加法的計算量差不多，但省去了重疊相加法最后的相加運算。 第75頁/共91頁第76頁/共91頁 4）用FFT計算相關

43、函數 相關的概念很重要，互相關運算廣泛應用于信號分析與統計分析，如通過相關函數峰值的檢測測量兩個信號的時延差等。 兩個長為N的實離散時間序列 x（n）與y（n）的互相關函數定義為 1010)()()()()(NnNnxynynxnynxr第77頁/共91頁則可以證明，rxy()的離散傅里葉變換為 Rxy（k）=X*(k)Y(k) 其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy() ， 0kN-1第78頁/共91頁互相關函數定義為 x(n)及y(n)的卷積公式 相比較，我們可以得到相關和卷積的時域關系： 1010)()()()()(NnNnxymnynxnymnxmr)()()()()(10mymxnynmxmfNn)()()()()()()(1010mymxnynmxnymnxmrNnNnxy第7