1、3.3 快速傅里葉變換快速傅里葉變換 （FFT） 雖然頻譜分析和DFT運算很重要，但在很長一段時間里，由于DFT運算復雜，并沒有得到真正的運用。直到1965年庫利（Cooley）和圖基（Tukey）首次提出DFT運算的一種快速算法以后，情況才發(fā)生了根本變化，人們開始認識到DFT運算的一些內(nèi)在規(guī)律，從而很快地發(fā)展和完善了一套一套高速有效的運算方法高速有效的運算方法快速付里變換（快速付里變換（FFT）算法）算法。3.3.1減少運算量的思路設設有限長序列 x（n）長度為N，則：101, 1 , 0)()()(NnnkNNkwnxnxDFTkX因x（n）和 都是復數(shù)，因此，每計算一個X（k）值，要進行

2、N次復數(shù)相乘，和N-1次復數(shù)相加，整個DFT運算需要: N*N次復數(shù)相乘 N*(N-1)次復數(shù)相加。nkNw可見，在DFT計算中，不論是乘法和加法，運不論是乘法和加法，運算量均與算量均與N2成正比成正比。因此，N較大時，運算量十分可觀。例，計算N=10點的DFT，需要100次復數(shù)相乘，而N=1024點時，需要1048576（一百多萬）次復數(shù)乘法，如果要求實時處理，則要求有很高的計算速度才能完成上述計算量。 反變換IDFT與DFT的運算結(jié)構(gòu)相同，只是多乘一個常數(shù)1/N，所以二者的計算量相同。FFT算法的基本思想：算法的基本思想：兩個特性可減少運算量：1）系數(shù) 是一個周期函數(shù)，它的周期性和對稱性可

3、用來改進運算，提高計算效率。 例 又如 因此 利用這些周期性和對稱性周期性和對稱性，使DFT運算中有些項可合并；nkNjnkNew2nkNnNkNkNnNwww)()(, 12/NNwkNNkNww)2/( 2）利用 的周期性和對稱性，把長度為N點的大點數(shù)的DFT運算依次分解為若干個小點數(shù)的DFT。 因為DFT的計算量正比于N2，N小，計算量也就小。 FFT算法正是基于這樣的基本思想發(fā)展起來的。它有多種形式，但基本上可分為兩類：時間抽取法和頻率抽取法。nkNw3.3.2基2-FFT算法假定N是2的整數(shù)次方，N=2M，M：正整數(shù) 首先將序列x（n）分解為兩組，一組為偶數(shù)項，一組為奇數(shù)項， r=0

4、,1，N/2-1 )() 12()()2(21rxrxrxrx（1）按時間抽取的按時間抽取的FFT2011)()(NnNnnkNnkNwnxwnx偶數(shù)奇數(shù)10)()()(NnnkNwnxnxDFTkx求x(n)的DFT：12/0)12(12/02) 12()2(NrkrNNrrkNwrxwrx12/0212/02) 12()2(NrrkNNrkNrkNwrxwwrxnNnNjnNjnNweew222222 /2 1/2 1002212(2 )(21)NNrkkrkNNNrrkNX kxr WWxrWXkW Xk /2 1102(2 )NrkNrXkxr W /2 1202(21)NrkNrXk

5、xrW注意到，X1（k）， X2（k）只有N/2個點，即k=0，1，N/2-1，還必須應用系數(shù) wkN 的周期性和對稱性表示 X（k）的 N/2 N-1點：rkNkNrNww22/2 12XX(),0,1,122kNNNX kkWkkkNNkNWW)2(可見，一個N點的DFT被分解為兩個N/2點的DFT，這兩個N/2點的DFT再合成為一個N點DFT 12,0,1,12kNNX kXkW Xkk 12(),0,1,122kNNNX kXkW Xkk結(jié)論：一個結(jié)論：一個N點的點的DFT的運算可以分解為兩個的運算可以分解為兩個N/2點點的的DFT運算運算 ，再經(jīng)上述算法計算得到。，再經(jīng)上述算法計算得

6、到。蝶形信號流圖蝶形信號流圖將X1（k）和X2（k）合成X(k)運算可歸結(jié)為： bWabWaN=8為例子為例子 兩個兩個4點點DFT組合成一個組合成一個8點點DFT 按照這個辦法，由于N=2M，仍然是偶數(shù)，可以被整除，因此可以對兩個N/2點的DFT再分別作進一步的分解。即對X1（k）， X2（k）的計算，又可以分別通過計算兩個長度為N/4=2點的DFT，進一步節(jié)省計算量。這樣，一個點的DFT就可以分解為四個點的DFT。 4個個2點點DFT組合成組合成2個個4點點DFT 最后剩下的是2點DFT，它可以用一個蝶形結(jié)表示：) 1 ()0() 1 ()0() 1 () 1 ()0() 1 ()0()0

7、(1202xWxxWxXxWxxWxXoNoN這樣，一個8點的完整的按時間抽取運算的流圖基基2-DIT的的8點點DFT運算流圖運算流圖 由于這種方法每一步分解都是按輸入時間序列是屬于偶數(shù)還是奇數(shù)來抽取的，所以稱為“時域抽取法”或“時間抽取法時間抽取法”。時間抽取法時間抽取法FFT的運算特點：的運算特點：1）蝶形運算）蝶形運算2）同址計算）同址計算 3）序數(shù)重排）序數(shù)重排4）蝶形類型隨迭代次數(shù)成倍增加）蝶形類型隨迭代次數(shù)成倍增加 1）蝶形運算）蝶形運算對于N=2M，總是可以通過M次分解最后成為2點的DFT運算。這樣構(gòu)成從x(n)到X(k)的M級運算過程。每一級運算都由N/2個蝶形運算個蝶形運算構(gòu)

8、成。因此每一級運算都需要N/2次復乘和N次復加，經(jīng)過時間抽取后M級運算總共需要的運算： 復乘 復加 NNMN2log22NNMN2log 而直接運算時則與N2 成正比。例N=2048，N2=4194304，(N/2)log2N=11264，N2 / (N/2)log2N=392.4。FFT顯然要比直接法快得多。2）同址同址計算計算 當數(shù)據(jù)輸入到存儲器中以后，每一級運算的結(jié)果仍然儲存在同一組存儲器中，直到最后輸出，中間無需其它存儲器，這叫同址同址計算計算。每一級運算均可在同址存儲器同址存儲器進行，這種同址同址運算結(jié)構(gòu)可節(jié)省存儲單元，降低設備成本，還可節(jié)省尋址的時間。3）序數(shù)重排）序數(shù)重排 對按時

9、間抽取FFT的原位運算結(jié)構(gòu)，當運算完畢時，正好順序存放著正好順序存放著 X（0），X（1），X（2），X（7），因此可直接按順序輸出按順序輸出。但這種原位運算的輸入 x（n）卻是按x(0)，x(4)，x(2)，x(6)，x(7)的順序存入存儲單元，這種順序看起來相當雜亂，然而它也是有規(guī)律的。當用二進制表示這個順序時，它正好是“碼位倒置碼位倒置”的順序。例如，原來的自然順序應是 x(1)的地方，現(xiàn)在放著 x(4)，用二進制碼表示這一規(guī)律時， 則是在 x（0 0 1）處放著 x（1 0 0）， x（0 1 1）處放著 x（1 1 0）。對于原位運算結(jié)構(gòu)，偶、奇分解過程如下： 表 碼位倒置順序自然順

10、序二進碼表示碼位倒置碼位倒序0000000010011004201001023011110641000011510110156110010371111117 在實際運算中，一般直接將輸入數(shù)據(jù) x（n）按碼位倒置的順序排好輸入很不方便，總是先按自然順序輸入存儲單元，然后再通過變址運算將通過變址運算將自然順序的存儲轉(zhuǎn)換成碼位倒置順序的存儲自然順序的存儲轉(zhuǎn)換成碼位倒置順序的存儲，然后進行FFT的原位計算。Matlab的bitrevorder（x）函數(shù)就是實現(xiàn)自然序的數(shù)據(jù)排列調(diào)整成倒序排列。 目前有許多通用DSP芯片支持這種碼位倒置碼位倒置的尋址功能。也可以通過移動流圖中的某些整行支路 ，把輸入x(n

11、)調(diào)整成自然序，但是輸出頻譜序列則成為倒序排列。輸出頻譜是倒序排列的，可以通過倒序整理環(huán)節(jié)，成為自然序。但不利于編程和效率，意義不大。 4）蝶形類型隨迭代次數(shù)成倍增加）蝶形類型隨迭代次數(shù)成倍增加觀察8點FFT的三次迭代運算：第一級迭代，有一種類型的蝶形運算系數(shù)W08，兩個數(shù)據(jù)點間隔為1。第二級迭代，有二種類型的蝶形運算系數(shù)W08、W28，參加運算的兩個數(shù)據(jù)點間隔為2。第三級迭代，有四類蝶形運算系數(shù)W08、W18、W28、W38，參加運算的兩個數(shù)據(jù)點間隔為4。 結(jié)論：每迭代一次，蝶形類型增加一倍，數(shù)據(jù)點間隔也增大一倍。每一級的取數(shù)間隔和蝶形類型種類均為2i-1，i=1,2,M。 （2）基）基2頻

12、率抽取的頻率抽取的FFT（DIF）對于N=2M情況下，輸入序列按前后對半分開，這樣便將N點DFT寫成前后兩部分：12/012/)()()(NnNNnnkNnkNWnxWnxkX12/012/0)2()2()(NnNnkNnNnkNWNnxWnx12/0)2/()2/()(NnnkNkNNWNnxWnx把X（k）進一步分解為偶數(shù)組和奇數(shù)組：奇數(shù)為偶數(shù)1, 1) 1(, 1)2/(2/kWWkkNNNN12/0)2/()1()()(NnnkNkWNnxnxkX12/02)2/()()2(NnnrNWNnxnxrX12/02/)2/()(NnnrNWNnxnx12/0) 12()2/()(1)X(2

13、rNnrnNWNnxnx12/02/)2/()(NnnrNnNWWNnxnx令 a(n)=x(n)+x(n+N/2) b(n)=x(n)-x(n+N/2wnN這兩個序列都是N/2點的序列，將其代入上兩式，得12/02/)()2(NnnrNWnarXnrNNnWnbrX2/12/0)() 12(這正是兩個N/2點的DFT運算，即將一個N點的DFT分解為兩個N/2點的DFT，上式的運算關(guān)系可用下圖表示. x(n)x(n+N/2)a(n)=x(n)+x(n+N/2)b(n)=x(n)-x(n+N/2)WNnWNn1比較發(fā)現(xiàn)，比較發(fā)現(xiàn)，DIF和和DIT的流圖是流圖轉(zhuǎn)置關(guān)系，即把輸入名和的流圖是流圖轉(zhuǎn)置

14、關(guān)系，即把輸入名和輸出名對調(diào)，流向倒轉(zhuǎn)但增益和結(jié)構(gòu)不變。輸出名對調(diào)，流向倒轉(zhuǎn)但增益和結(jié)構(gòu)不變。 因此，兩者運算量完全一樣。因此，兩者運算量完全一樣。3.3.3 N為組合數(shù)的為組合數(shù)的FFT算法算法 以上討論的都是以2為基數(shù)的FFT算法，即N=2M。 實際應用時，有限長序列的長度N很大程度上由人為因素確定，因此多數(shù)場合可取 N=2M，從而直接使用以 2 為基數(shù)的FFT算法。 如N不能人為確定 , N的數(shù)值也不是以2為基數(shù)的整數(shù)次方，處理方法處理方法有兩種: 補零: 將x(n)補零 , 使N=2M. 例如N=30，補上x(30)=x(31)=0兩點,使N=32=25,這樣可直接采用以2為基數(shù)M=5

15、的FFT程序。 有限長度序列補零后并不影響其頻譜 X(ejw) ，只是頻譜的采樣點數(shù)增加了,上例中由30點增加到32點，所以在許多場合這種處理是可接受的。 如要求準確的N點DFT值，可采用任意數(shù)為基數(shù)的 FFT 算法 , 其 計算效率低于以2為基數(shù)FFT算法。 如 N 為復合數(shù)，可分解為兩個整數(shù)兩個整數(shù)P與與Q的乘積的乘積，像前面以2為基數(shù)時一樣，F(xiàn)FT的基本思想是將DFT的點數(shù)盡量分小，因此，在N=PQ情況下，也希望將N點的DFT分解為P個Q點DFT或Q個P點DFT，以減少計算量。步驟：kPkknQnn10,kn分別為0, 1,，Q-1; 01,kn分別為0, 1,，P-1。 N點DFT可以

16、重新寫成為 100101)(,)()(NnknNWnxkkXkPkXkX1010)(01010101QnPnnQnkPkNWnQnxQnPnnkNPnkNQnkNQnPnnkNPnkNQnkNPQnkNWWWnnxWWWWnnxkkX,考慮到 再令1010nkPQnkNWW令1010010100100110,QnnkQnkNPnnkPWWWnnxkkXPnnkPWnnxnkX,01000,0011001nkQnkNQnWWnkXkkX00001001,nkNWnkXnkXQnnkQWnkXkkX,以P=3，Q=4, N=12為例 (1) 先將 x(n)通過 x(n1Q+n0)改寫成 x(n1,

17、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)(2) 求個點的DFT （）X1(k0,n0)乘以 得到X1(k0,n0)。（）求P個點的DFT，參變量是k020301001110,nnkWnnxnkX304001102001,nnkWnkXkkX（）將X2(k0,k1)

18、通過X(k0+k1P)恢復為X(k)00nkNW運算量：（）求個P點DFT需要P2次復數(shù)乘法和P(P-1)次復數(shù)加法；（）乘個W因子需要次復數(shù)乘法；（）求P個點DFT需要PQ2 次復數(shù)乘法和P（-1)次復數(shù)加法?？偟膹蛿?shù)乘法量: QP2+N+PQ2=N(P+Q+1)；總的復數(shù)加法量: QP(P-1)+PQ(Q-1)=N(P+Q-2) N=667=23*29，乘法：667*（23+29+1）=35351 加法：667*（23+29-2）=33350N=667，乘法：667*667=444889； 12.58 加法：667*（667-1）=444222； 13.32N=1024；乘法：1024/2

19、*log2（1024）=5120； 加法： 1024*log2（1024） =10240； 上述分解原則可推廣至任意基數(shù)的更加復雜的情況。例如, 如果N可分解為m個質(zhì)數(shù)因子p1,p2,，pm，即 N=p1p2p3pm則：（1）可把N先分解為兩個因子N=p1q1，其中q1=p2p3pm ，并用上述討論的方法將DFT分解為p1個q1點DFT；（2）將q1分解為q1=p2q2，q2=p3p4pm，然后將每個q1 點DFT再分解為p2個q2點DFT；： 依此類推，：通過m次分解，一直分到最少點數(shù)的DFT運算，從而獲得最高的運算效率。其運算量近似為N(p1+p2+pm)次復數(shù)乘法和復數(shù)加法。但計算效率的

20、提高是要以編程的復雜性為代價的，一般較少應用。 p1=p2 = =pm =2，為基基2 FFT 算法算法。當組合數(shù) N=P1P2P3Pm中所有的Pi均為4時，就是基基4FFT算法算法 以N=43為例,第一級運算的一般形式為： ), 3(), 2(), 1 (), 0(111111111111), 3(), 2(), 1 (), 0(), 3(), 2(), 1 (), 0(212121210101010194643404644424043424140404040404011011011011nnxnnxnnxnnxjjjjnnxnnxnnxnnxWWWWWWWWWWWWWWWWnnXnnXnn

21、XnnX30,),(),(04300120101022kWnnnxnnkXknn以N=1024=210，基2FFT 5120 復乘， 基4FFT N*（M-1）=4096 復乘3.3.4 Chirp-z變換變換 采用FFT可以算出N點DFT值，即X（z）在z平面單位圓上的N點等間隔取樣值。問題的提出： 不需要整個單位圓上的取樣，如：對于窄帶信號，只需對信號所在頻段進行分析，希望頻譜的采樣集中在這一頻帶內(nèi)，以獲得較高的分辨率。 需要對其它圍線上的z變換取樣，例如：語音信號處理中，需要知道極點所在頻率，如果極點離單位圓較遠，則其單位圓上的頻譜就很平滑，不易判斷；如果采樣可以沿一條接近這些極點的弧線

22、進行，則在極點所在頻率上的頻譜將出現(xiàn)明顯的尖峰，由此可較準確地測定極點頻率。 要求能有效地計算N是素數(shù)時，序列的DFT。 算法原理算法原理： 已知x（n），0nN-1令zk=AW-k ，k=0，M-1 ，M ：采樣點數(shù)，A、W：任意復數(shù)00jojoeWWeAA其中：A0表示起始取樣點的半徑長度，通常A010表示起始取樣點z0的相角0表兩相鄰點之間的等分角W0螺旋線的伸展率，W01則線內(nèi)縮（反時針），W0=1則表示半徑為A0的一段圓弧，若A0=1，這段圓弧則是單位圓的一部分。圖 螺旋線采樣 計算z變換在采樣點 zk 的值 10)()(NnnkkznxzX nk可以用以下表示式來替換10)()(N

23、nnknkWAnxzX)(21222nknknk顯然，按照以上公式計算出全部M點采樣值需要 NM 次復乘和（N-1）M次復加，當N及M較大時，計算量迅速增加，以上運算可轉(zhuǎn)換為卷積形式，從而可采用FFT進行，這樣可大大提高計算速度。k=0,1, ,M-12222)(1022)()(nkNnnnkkWWAnxWzX則222)(,)()(2nnnWnhWAnxng定義：1, 1 , 0)()()()()(210222MkkhkgWnkhngWzXkNnkk則)()()()(2222202202njnonjoneWeWWnhx22)(nWnhxx(n)g(n)X(zk)22nnWA22nW系統(tǒng)的單位脈

24、沖響應頻率隨時間成線性增加的線性調(diào)頻信號相似，因此稱為Chirp -z變換變換。 算法實現(xiàn)算法實現(xiàn)： 由于輸入信號 g（n）是有限長的，長為N，但序列 是無限長的，而計算 0M-1 點卷積 g（k）*h（k）所需要的 h（n）是取值在）是取值在n=-（N-1）M-1那一部分的值那一部分的值，因此，可認為h（n）是一個有限長序列，長為L=N+M-1。 所以，Chirp -z變換為兩個有限長序列的線性卷積g（k）*h（k），可用通過FFT來實現(xiàn)。 h（n）的主值序列 可由h（n）作周期延拓后取0nL-1部分值獲得，將 與g（n）作圓周卷積后，其輸出的前M個值就是Chirp -z變換的M個值。22)

25、(nWnh)(nh)(nhChirp-z變換的計算步驟： (1) 求h（n）的主值序列 (2) 用FFT求 的傅里葉變換: H(k)=FFT , L點 (3) 對x（n）加權(quán)并補零1110)(2/)(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利用FFT計算Chirp-z變換計算量估算：乘法數(shù)估計 （1）（2）兩步可以事先計算，不必實時計算。 （3）（7）兩步兩次加權(quán)共

26、計N+M次復乘,形成Y（k）需L次復乘，一個FFT與一個IFFT需Llog2L次乘，所以總乘法數(shù)：L+N+M+Llog2L，直接計算乘法數(shù)， NM N及M較大時，用FFT實現(xiàn)Chirp-Z變換，速度上有很大的改進。 Chirp-z變換的特點：變換的特點：1）輸入序列的長度N 與 輸出序列的長度 M 不需要相等；2）N 及 M 不必是高度復合數(shù)，二者均可為素數(shù)；3）相鄰采樣點 zk 之間的角間隔 0 是任意的，即頻率分辨率是任意的；4）圍線是任意的，不必是 Z 平面上的圓；5）起始點 z0 可任意選定，即可從任意頻率上開始對輸入數(shù)據(jù)進行窄帶高分辨率分析；6）若 A=1 , M=N , , 可用

27、Chirp-z變換變換計算DFT（即使 N 為素數(shù)）。NjeW23.4 離散傅立葉變換的實際應用問題離散傅立葉變換的實際應用問題采樣采樣截短截短DFT)(ajXa()x t( )xn)(jeX( )( )( )NNx nxnR n)(jNeX)(kXNNkjNeX/2)( 混迭：對連續(xù)信號x(t)進行數(shù)字處理前，要進行采樣nanTttxtx)()()( 采樣序列的頻譜是連續(xù)信號頻譜的周期延拓，周期為fs，如采樣率過低，不滿足采樣定理，fs32時，上述計算線性卷積的方法比直接計算線卷積有明顯的優(yōu)越性,因此,也稱上述循環(huán)卷積方法為快速卷積法。 上述結(jié)論適用于 x(n)、h(n) 兩序列長度比較接近

28、或相等的情況，如果x(n)、h(n) 長度相差較多，例如， h(n) 為某濾波器的單位脈沖響應，長度有限，用來處理一個很長的輸入信號 x(n)，或者處理一個連續(xù)不斷的信號，按上述方法， h(n) 要補許多零再進行計算，計算量有很大的浪費，或者根本不能實現(xiàn)。 為了保持快速卷積法的優(yōu)越性，可將 x(n) 分為許多段，每段的長度與 h(n) 接近 ，處理方法有兩種：(1) 重疊相加法由分段卷積的各段相加構(gòu)成總的卷積輸出 假定 xi(n) 表示 x（n）序列的第i段 ： 則輸入序列可表為：于是輸出可分解為： 其中 01) 1()()(22NiniNnxnxiiinxnx)()(iiiinynhnxnh

29、nxny)()(*)()(*)()()(*)()(nhnxnyii 1）先對 h（n）及 xi（n）補零，補到具有L點長度，L=N1+N2-1。 一般選 L=2M。 2）用基2 FFT計算 yi（n）=xi（n）*h（n）。 3）重疊部分相加構(gòu)成最后的輸出序列。 由于 yi（n）的長度為L，而xi（n）的長度為N2，因此相鄰兩段 yi（n）序列必然有L-N2=N1-1點發(fā)生重疊。 iinyny)()(計算步驟：a. 事先準備好濾波器參數(shù) H（k）=DFTh(n)，L點b. 用N點FFT計算Xi（k）=DFTxi(n)c. Yi(k)=Xi(k)H(k)d. 用N點IFFT求 yi（n）=IDF

30、TYi(k)e. 將重疊部分相加(2)重疊保留 這種方法和第一種方法稍有不同，即將上面分段序列中補零的部分不是補零，而是保留原來的輸入序列保留原來的輸入序列值值，這時，如利用DFT實現(xiàn)h(n)和xi(n)的循環(huán)卷積，則每段卷積結(jié)果中有N1-1個點不等于線性卷積值需舍去。 重疊保留法與重疊相加法的計算量差不多，但省去了重疊相加法最后的相加運算。 相關(guān)的概念很重要，互相關(guān)運算廣泛應用于信號分析與統(tǒng)計分析，如通過相關(guān)函數(shù)峰值的檢測測量兩個信號的時延差等。 兩個長為N的實離散時間序列 x（n）與y（n）的互相關(guān)函數(shù)定義為 1010)()()()()(NnNnxymnynxnymnxmr3.5.4相關(guān)函

31、數(shù)的FFT計算互相關(guān)不具有交換特性 110010()() ( )( ) ()() ( )()NNyxnnNxynrmy nm x ny n x nmx nmy nrm 則可以證明，rxy()的離散傅里葉變換為 Rxy（k）=X*(k)Y(k) 其中 X(k)=DFTx(n), Y(k)=DFTy(n), Rxy(k)=DFTrxy() ， 0kN-1證明：互相關(guān)函數(shù)定義為 x(n)及y(n)的卷積公式 相比較，我們可以得到相關(guān)和卷積的時域關(guān)系： 1010)()()()()(NnNnxymnynxnymnxmr)()()()()(10mymxnynmxmfNn)()()()()()()(1010

32、mymxnynmxnymnxmrNnNnxy)()()()(1)()()(10102kYkXIDFTekYkXNnynxrNnNkkNjxy得： Rxy（k）=X*(k)Y(k) )(*)()(DFTkXnRnxNN因 當x（n）=y（n）時，得到x（n）的自相關(guān)函數(shù)自相關(guān)函數(shù)為： 維納辛欽定理： 自相關(guān)函數(shù)與信號功率譜互為傅立葉變換對。 10)()()(NnxxnxnxrkNjNkekXN2102| )(|1 上面的推導表明：(1)互相關(guān)和自相關(guān)函數(shù)的計算可利用FFT實現(xiàn)。(2) 由于離散傅里葉變換隱含著周期性，所以用FFT計算離散相關(guān)函數(shù)也是對周期序列而言的。直接做N點FFT相當于對兩個N

33、點序列x(n)、y(n)作周期延拓，作相關(guān)后再取主值（類似圓周卷積）。 實際一般要求的是兩個有限長序列的線性相關(guān)，為避免混淆，需采用與循環(huán)卷積求線性卷積相類似的方法，先將序列延長補0后再用上述方法。利用FFT求兩個有限長序列的線性相關(guān)： 設x(n)和y(n)的長均為N，求線性相關(guān)；（1）為了使兩個有限長序列的線性相關(guān)可用其循環(huán)相關(guān)代替而不產(chǎn)生混淆，選擇周期L2N-1，且L=2m，以使用FFT，將x(n)，y(n)補零至長為L。 （2） 用FFT計算 X(k)，Y(k)（k=0，1，L-1）（3） 計算 R（k）=X*（k）Y（k）（4） 對R（k）作 IFFT，取后N-1項，得 取前N項，得 11)(mNmrxy10)(Nmmrxy解：matlab實現(xiàn)代碼如下： x=1 3 -1 1 2 3 3 1;y=2 1 -1 1 2 0 -1 3;k=length(x);xk=fft(x,2*k);yk=fft(y,2*k);rm=real(ifft(conj(xk).*yk);rm=rm(k+2:2*k) rm(1:k);m=(-k+1):(k-1);stem(m,rm)xlabel(m); ylabel(幅度);【例3.5.1】假設x=1 3 -1 1 2 3 3 1，y=2 1 -1 1 2 0 -1 3，求它們的互相關(guān)函數(shù)?！纠?.5.2】