實驗三快速傅里葉變換及其應用

1、DSP 試驗040080121DSP 試驗04008012實驗三快速傅里葉變換及其應用:實驗目的(1) 加深對FFT的理解，熟悉 matlab中的有關(guān)函數(shù)。(2) 應用FFT對典型信號進行頻譜分析。(3) 了解應用FFT進行信號頻譜分析過程中可能出現(xiàn)的問題，以便在實際中正確應用FFT.(4) 應用FFT實現(xiàn)序列的線性卷積和相關(guān)。Fouier 變換(DF。這x(n)的長度為N時,:實驗原理：在各種信號序列中，有限長序列信號處理占有很重要地位，對有限長序列，我們可以使用離散一變換不但可以很好的反映序列的頻譜特性，而且易于用快速算法在計算機上實現(xiàn)，當序列它的DFT定義為:反變換為:有限長序列的DFT

2、是其Z變換在單位圓上的等距采樣，或者說是序列Fourier變換的等距采樣，因此可以用于序列的譜分析。FFT并不是與DFT不同的另一種變換，而是為了減少DFT運算次數(shù)的一種快速算法。它是對變換式進行一次次分解， 使其成為若干小點數(shù)的組合，從而減少運算量。常用的FFT是以2為基數(shù)的，其長度二一。它的效率高，程序簡單，使用非常方便，當要變換的序列長度不等于 2的整數(shù)次方時，為了使用以2為基數(shù)的FFT,可以用末位補零的方法， 使其長度延長至2的整數(shù)次方。(一)在運用DFT進行頻譜分析的過程中可能的產(chǎn)生三種誤差(1) 混疊序列的頻譜是被采樣信號的周期延拓，當采樣速率不滿足 Nyquist定理時，就會發(fā)生

3、頻譜混疊， 使得采樣后的信號序列頻譜不能真實的反映原信號的頻譜。避免混疊現(xiàn)象的唯一方法是保證采樣速率足夠高，使頻譜混疊現(xiàn)象不致出現(xiàn)，即在確定采樣頻率之前，必須對頻譜 的性質(zhì)有所了解，在一般情況下，為了保證高于折疊頻率的分量不會出現(xiàn)，在采樣前，先用低通模擬濾波器對信號進 行濾波。(2) 泄漏實際中我們往往用截短的序列來近似很長的甚至是無限長的序列，這樣可以使用較短的 DFT來對信號進行頻譜分析，這種截短等價于給原信號序列乘以一個矩形窗函數(shù)，也相當于在頻域?qū)⑿盘柕念l譜和矩形窗函數(shù)的頻譜卷積，所得的 頻譜是原序列頻譜的擴展。泄漏不能與混疊完全分開，因為泄漏導致頻譜的擴展，從而造成混疊。為了減少泄漏的

4、影響，可以選擇適當?shù)拇昂?數(shù)使頻譜的擴散減至最小。(3) 柵欄效應DFT是對單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個連續(xù)函數(shù)，就一定意義上看，用DFT來觀察頻譜就好像通過一個柵欄來觀看一個圖景一樣，只能在離散點上看到真實的頻譜，這樣就有可能發(fā)生一些頻譜的峰點或谷 點被“尖樁的柵欄”所攔住，不能別我們觀察到。減小柵欄效應的一個方法就是借助于在原序列的末端填補一些零值，從而變動DFT的點數(shù)，這一方法實際上是人為地改變了對真實頻譜采樣的點數(shù)和位置，相當于搬動了每一根“尖樁柵欄”的位置，從而使得頻譜的峰點或谷點暴露 出來。(二 )用FFT計算線性卷積用FFT可以實現(xiàn)兩個序列的圓周卷積。在一

5、定的條件下，可以使圓周卷積等于線性卷積。一般情況，設(shè)兩個序列的長度分別為 N1和N2,要使圓周卷積等于線性卷積的充要條件是FFT的長度NA N1 + N2對于長度不足N的兩個序列，分別將他們補零延長到N。當兩個序列中有一個序列比較長的時候，我們可以采用分段卷積的方法。有兩種方法：(1) 重疊相加法。將長序列分成與短序列相仿的片段，分別用FFT對它們作線性卷積，再將分段卷積各段重疊的部分 相加構(gòu)成總的卷積輸出。(2) 重疊保留法。這種方法在長序列分段時，段與段之間保留有互相重疊的部分，在構(gòu)成總的卷積輸出時只需將各段線性卷積部分直接連接起來，省掉了輸出段的直接相加。(三)用FFT計算相關(guān)函數(shù)兩個長

6、為N的實離散時間序列x(n)與y(n)的互相關(guān)函數(shù)定義為：N JN 二rxy(m)二 x(n)y(n m)二 x(n - m)y(n) =x(m) y(m)n 衛(wèi)m z0rxy(n)的離散傅里葉變換為：Rxy (k)二 X (k)Y(k),0 _ k 一 N _1N->4 N_12 j%rxx(m)x(n)x(n + m)=需送 X(k) e N當x(n) = y(n)時，得到x(n)的自相關(guān)函數(shù)為：n=0N km利用FFT求兩個有限長序列線性相關(guān)的步驟(設(shè) x(n)長Ni , y(n)長N2)：(1)為了使兩個有限長序列的線性相關(guān)可用其圓周相關(guān)代替而不產(chǎn)生混淆，選擇周期N =2l a

7、NiN2-1 ,以便使用FFT將x(n), y(n)補零至長為N。用FFT計算X(k),Y(k),0沖訃-1(3) Rxy(k) =X (k)Y(k),0 沖 沙一1(4) 對Rxy(k)作ifft；取后N -1項，得rxy(m) - N 仁m豈-1 ；取前N項，得rxy(m)。汀訃-1。三、實驗內(nèi)容及步驟 實驗中用到的信號序列:a)高斯(Gaussian)序列f -p)2xa(n)=e q 0勿 <150 其它b)衰減正弦序列en sin(2竈fn)Xb( n) =,0n <15其它2DSP 試驗040080123DSP 試驗04008012c)三角波序列d)反三角波序列#DSP

8、 試驗04008012I nxc (n) = <8 n00 _ n _34乞n空7其它4 nxd (n) = <n 400 _ n _ 34空n空7其它上機實驗內(nèi)容：(1)觀察高斯序列的時域和幅頻特性，固定信號xa(n)中參數(shù)p=8,改變q的值，使q分別等于2, 4, 8，觀察它們的時域和幅頻特性，了解當q取不同值時，對信號序列的時域幅頻特性的影響；固定q=8，改變p，使p分別等于8，13,14,觀察參數(shù)p變化對信號序列的時域及幅頻特性的影響，觀察p等于多少時，會發(fā)生明顯的泄漏現(xiàn)象，混疊是否也隨之出現(xiàn)？記錄實驗中觀察到的現(xiàn)象，繪出相應的時域序列和幅頻特性曲線。參數(shù)p=8不變主要代碼

9、如下：i=1:15;p=8;q=2;subplot(3,2,1);x(i)=exp(-(i-p).A2)/q);stem(x);xlabel(' n');subplot(3,2,2);G=fft(x);plot(abs(G(1:15);xlabel('k');運行結(jié)果如下：4DSP 試驗04008012#DSP 試驗04008012，波形變胖，低頻分量變多，頻域信號頻譜泄露程對比三張圖可知：P不變，隨著q值的增大，時域信號幅值變化緩慢時域幅度對應變大#DSP 試驗04008012度減小。參數(shù)q=8不變代碼和之前相比只要改變相應幾個參數(shù)而已，故不列出;運行結(jié)果如下

10、：可見，當q不變，隨著p的增大，時域信號幅值不變，會在時間軸移位，對應右移，可見p決定了波形位置,當實驗中q=8,p=13時，明顯出現(xiàn)泄漏。觀察衰減正弦序列 xb(n)的時域和幅頻特性，a=0.1 , f=0.0625，檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出 幅頻特性曲線，改變f,使f分別等于0.4375和0.5625，觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)位置，有無混疊和 泄漏現(xiàn)象？說明產(chǎn)生現(xiàn)象的原因。其中一段代碼如下(其他的都是改變參數(shù)而已)：a=0.1;f=0.0625;for i=1:16x(i)=exp(-a*(i-1)*si n( 2*pi*f*(i-1);endfor i

11、=17:100x(i)=0;endn=0:15;subplot(3,2,1);plot( n(1:16),x(1:16);xlabel(' n');subplot(3,2,2);G=fft(x,16);plot( n(1:16),abs(G(1:16);xlabel('k');結(jié)果：7DSP 試驗04008012#DSP 試驗04008012#DSP 試驗04008012#DSP 試驗04008012觀察三角波和反三角波序列的時域和幅頻特性,用N=8點FFT分析信號序列xc(n)和xd(n)的幅頻特性，觀察兩者的序8DSP 試驗04008012列形狀和頻譜曲線有

12、什么異同？繪出兩序列及其幅頻特性曲線。 主要代碼如下：for i=1:4x(i)=i-1;endfor i=5:8x(i)=9-i;endn=0:7;subplot(2,2,1);plot( n,x(1:8);xlabel(' n');subplot(2,2,2);G=fft(x,8);plot( n(1:8),abs(G(1:8);xlabel('k');結(jié)果：在xc(n)和xd(n)末尾補零，用 N=32點FFT分析這兩個信號的幅頻特性，觀察幅頻特性發(fā)生了什么變化？兩情況的FFT 頻譜還有相同之處嗎？這些變化說明了什么？ 代碼只需改動幾處數(shù)據(jù)即可，其中一段如

13、下所示: for i=1:4x(i)=i-1;endfor i=5:8x(i)=9-i;endfor i=9:32x(i)=0;endn=0:31;subplot(2,1,1);G=fft(x,32);plot( n(1:32),abs(G(1:32);xlabel('k');運行結(jié)果：變化：反三角波的低頻分量增多，對信號末尾補零加長整數(shù)個周期可以對原信號達到細化頻譜的作用。一個連續(xù)信號含兩個頻率分量，經(jīng)采樣得x(n)=sin2 n *0.125n+cos2 n *(0.125+ 4 f)nF0,1-,N已知N=16, 分別為1/16和1/64，觀察其頻譜；當 N=128時，不

14、變，其結(jié)果有何不同，為什么？10DSP 試驗04008012其中一段代碼如下：實現(xiàn)方法和之前一樣N=16; f=1/16; for n=1:Nx( n)=si n( 2*pi*0.125*( n-1)+cos(2*pi*(0.125+f)*( n-1); end n=0:15;subplot(2,2,1);G=fft(x,16);plot( n(1:16),abs(G(1:16);xlabel('k');運行結(jié)果：(5)用FFT分別實現(xiàn)xa(n)( p = 8, q= 2)和xb(n) (a= 0.1, f = 0.0625)的16點循環(huán)卷積和線性卷積。x(i)=exp(-(i

15、-1-8)A2)/2);y(i)=exp(-0.1*(i-1)*si n(2*pi*0.0625*(i-1);end%for i=17:31% x(i)=0;% y(i)=0;%e nd n=0:30G1=fft(x,31);G2=fft(y,31); z=ifft(G1.*G2,31); subplot(2,1,1);plot( n(1:31),z(1:31);subplot(2,1,2);G1=fft(x,16);G2=fft(y,16);z=ifft(G1.*G2,16); n=0:15;plot( n(1:16),z(1:16);運行結(jié)果：12DSP 試驗04008012#DSP 試驗

16、04008012(6)產(chǎn)生一 512點的隨機序列xe(n),并用xc(n)和xe(n)作線性卷積，觀察卷積前后xe(n)頻譜的變化。要求將8段，分別采用重疊相加法和重疊保留法。xe( n)分成#DSP 試驗04008012xe=ra nd(1,512);n1=0:1:3;xc1= n1;n2=4:7;xc2=8-n2;xc=xc1,xc2;yn=zeros(1,519);for j=0:7xj=xe(64*j+1:64*(j+1);xak=fft(xj,71);xck=fft(xc,71);yn仁 ifft(xak.*xck);temp=zeros(1,519);temp(64*j+1:64*

17、j+71)=y n1;yn=yn+temp;en d;n=0:518;figure(1)subplot(2,1,1);plot( n,yn);xlabel(' n');ylabel('y( n)');subplot(2,1,2);plot( n,abs(fft(y n); xlabel('k');ylabel('Y(k)'); axis(0,600,0,300);xe=ra nd(1,512);k=1:7;xe1=k-k;xe_1=xe1,xe;yn_1= zeros(1,519);for j=0:7xj_1=xe_1(64*j+

18、1:64*j+71);xak_ 仁 fft(xj_1);xck_ 仁 fft(xc,71);yn1_仁 ifft(xak_1.*xck_1);temp_ 仁 zeros(1,519);temp_1(64*j+1:64*j+64)=y n1_1(8:71); yn_1=yn _1+temp_1;en d;n=0:518;figure(2) subplot(2,1,1);plot( n, yn_1);xlabel(' n');ylabel('y( n)'); subplot(2,1,2);plot( n,abs(fft(yn_1); xlabel('k

19、9;);ylabel('Y(k)'); axis(0,600,0,300);運行結(jié)果：(7)用FFT分別計算xa(n)(p=8,q=2)和xb(n) (a= 0.1, f = 0.0625)的16點循環(huán)相關(guān)和線性相關(guān)，問一共有多少種結(jié)果,它們之間有何共同點？n=0:1:15;xan=exp(-( n-8).A2/2);xbn=exp(-0.1* n).*si n(2*pi*0.0625* n);k=le ngth(xb n);xan 1=xa n zeros(1,k-1);xbn 1=xb n zeros(1,k-1);xak=fft(xa n1);xbk=fft(xb n1); rm=real(ifft(conj(xak).*xbk);rm1= rm(k+1:2*k-1) rm(1:k); m=(_k+1):(k_1);subplot(2,1,1);stem(m,rm1);xlabel(' n');subplot(2,1,2);16DSPM專 0400822xakufft(xan)八xbkufft(xbn)八 rmurea-(ifft(conj(