數(shù)字信號(hào)處理吳鎮(zhèn)楊實(shí)驗(yàn)三

1、實(shí)驗(yàn)三 快速傅里葉變換及其應(yīng)用04011344 王晨一、 實(shí)驗(yàn)?zāi)康?1) 在理論學(xué)習(xí)的基礎(chǔ)上，通過(guò)本實(shí)驗(yàn)，加深對(duì)FFT的理解，熟悉MATLAB中的有關(guān)函數(shù)。(2) 應(yīng)用FFT對(duì)典型信號(hào)進(jìn)行頻譜分析的方法。(3) 了解應(yīng)用FFT進(jìn)行信號(hào)頻譜分析過(guò)程中可能出現(xiàn)的問(wèn)題，以便在實(shí)際中正確應(yīng)用FFT。(4) 應(yīng)用FFT實(shí)現(xiàn)兩個(gè)序列的線性卷積和相關(guān)   二、 實(shí)驗(yàn)原理  在各種信號(hào)序列中，有限長(zhǎng)序列信號(hào)處理占有很重要地位，對(duì)有限長(zhǎng)序列，我們可以使用離散Fourier變換(DFT)。這一變換不但可以很好的反映序列的頻譜特性，而且易于用快速算法在計(jì)算機(jī)上實(shí)現(xiàn)，當(dāng)序列x(

2、n)的長(zhǎng)度為N時(shí)，它的DFT定義為：反變換為：有限長(zhǎng)序列的DFT是其Z變換在單位圓上的等距采樣，或者說(shuō)是序列Fourier變換的等距采樣，因此可以用于序列的譜分析。FFT并不是與DFT不同的另一種變換，而是為了減少DFT運(yùn)算次數(shù)的一種快速算法。它是對(duì)變換式進(jìn)行一次次分解，使其成為若干小點(diǎn)數(shù)的組合，從而減少運(yùn)算量。常用的FFT是以2為基數(shù)的，其長(zhǎng)度 。它的效率高，程序簡(jiǎn)單，使用非常方便，當(dāng)要變換的序列長(zhǎng)度不等于2的整數(shù)次方時(shí)，為了使用以2為基數(shù)的FFT，可以用末位補(bǔ)零的方法，使其長(zhǎng)度延長(zhǎng)至2的整數(shù)次方。     (一)、在運(yùn)用DFT進(jìn)行頻譜分析的過(guò)程中可能產(chǎn)生三種誤

3、差：    (1)    混疊    序列的頻譜時(shí)被采樣信號(hào)的周期延拓，當(dāng)采樣速率不滿足Nyquist定理時(shí)，就會(huì)發(fā)生頻譜混疊，使得采樣后的信號(hào)序列頻譜不能真實(shí)的反映原信號(hào)的頻譜。避免混疊現(xiàn)象的唯一方法是保證采樣速率足夠高，使頻譜混疊現(xiàn)象不致出現(xiàn)，即在確定采樣頻率之前，必須對(duì)頻譜的性質(zhì)有所了解，在一般情況下，為了保證高于折疊頻率的分量不會(huì)出現(xiàn)，在采樣前，先用低通模擬濾波器對(duì)信號(hào)進(jìn)行濾波。    (2)    泄漏   

4、; 實(shí)際中我們往往用截短的序列來(lái)近似很長(zhǎng)的甚至是無(wú)限長(zhǎng)的序列，這樣可以使用較短的DFT來(lái)對(duì)信號(hào)進(jìn)行頻譜分析，這種截短等價(jià)于給原信號(hào)序列乘以一個(gè)矩形窗函數(shù)，也相當(dāng)于在頻域?qū)⑿盘?hào)的頻譜和矩形窗函數(shù)的頻譜卷積，所得的頻譜是原序列頻譜的擴(kuò)展。泄漏不能與混疊完全分開(kāi)，因?yàn)樾孤?dǎo)致頻譜的擴(kuò)展，從而造成混疊。為了減少泄漏的影響，可以選擇適當(dāng)?shù)拇昂瘮?shù)使頻譜的擴(kuò)散減至最小。    (3)    柵欄效應(yīng)    DFT是對(duì)單位圓上Z變換的均勻采樣，所以它不可能將頻譜視為一個(gè)連續(xù)函數(shù)，就一定意義上看，用DFT來(lái)觀察頻譜就好像通

5、過(guò)一個(gè)柵欄來(lái)觀看一個(gè)圖景一樣，只能在離散點(diǎn)上看到真實(shí)的頻譜，這樣就有可能發(fā)生一些頻譜的峰點(diǎn)或谷點(diǎn)被“尖樁的柵欄”所攔住，不能別我們觀察到。    減小柵欄效應(yīng)的一個(gè)方法就是借助于在原序列的末端填補(bǔ)一些零值，從而變動(dòng)DFT的點(diǎn)數(shù)，這一方法實(shí)際上是人為地改變了對(duì)真實(shí)頻譜采樣的點(diǎn)數(shù)和位置，相當(dāng)于搬動(dòng)了每一根“尖樁柵欄”的位置，從而使得頻譜的峰點(diǎn)或谷點(diǎn)暴露出來(lái)。 (二)、用FFT計(jì)算線性卷積    用FFT可以實(shí)現(xiàn)兩個(gè)序列的圓周卷積。在一定的條件下，可以使圓周卷積等于線性卷積。一般情況，設(shè)兩個(gè)序列的長(zhǎng)度分別為N1和N2，要使圓周卷積

6、等于線性卷積的充要條件是FFT的長(zhǎng)度:NN1N2對(duì)于長(zhǎng)度不足N的兩個(gè)序列，分別將他們補(bǔ)零延長(zhǎng)到N。    當(dāng)兩個(gè)序列中有一個(gè)序列比較長(zhǎng)的時(shí)候，我們可以采用分段卷積的方法。有兩種方法：a.重疊相加法。將長(zhǎng)序列分成與短序列相仿的片段，分別用FFT對(duì)它們作線性卷積，再將分段卷積各段重疊的部分相加構(gòu)成總的卷積輸出。   b.重疊保留法。這種方法在長(zhǎng)序列分段時(shí)，段與段之間保留有互相重疊的部分，在構(gòu)成總的卷積輸出時(shí)只需將各段線性卷積部分直接連接起來(lái)，省掉了輸出段的直接相加。 (三)、用周期圖法(平滑周期圖的平均法)對(duì)隨機(jī)信號(hào)作譜分析   

7、; 實(shí)際中許多信號(hào)往往既不具有有限能量，由非周期性的。無(wú)限能量信號(hào)的基本概念是隨機(jī)過(guò)程，也就是說(shuō)無(wú)限能量信號(hào)是一隨機(jī)信號(hào)。周期圖法是隨機(jī)信號(hào)作譜分析的一種方法，它特別適用于用FFT直接計(jì)算功率譜的估值。    將長(zhǎng)度為N的實(shí)平穩(wěn)隨機(jī)序列的樣本x(n)再次分割成K段，每段長(zhǎng)度為L(zhǎng),即L=N/K。每段序列仍可表示為：xi(n)=x(n+(i-1)L)，0nL-1，1iK但是這里在計(jì)算周期圖之前，先用窗函數(shù)w(n)給每段序列xi(n)加權(quán)，K個(gè)修正的周期圖定義為其中U表示窗口序列的能量，它等于在此情況下，功率譜估計(jì)量可表示為三、 實(shí)驗(yàn)程序代碼及結(jié)果討論  實(shí)驗(yàn)

8、中用到的信號(hào)序列：a)       Gaussian序列b)      衰減正弦序列c)      三角波序列d)      反三角波序列    上機(jī)實(shí)驗(yàn)內(nèi)容：(1)、觀察高斯序列的時(shí)域和幅頻特性，固定信號(hào)xa(n)中參數(shù)p=8，改變q的值，使q分別等于2，4，8，觀察它們的時(shí)域和幅頻特性，了解當(dāng)q取不同值時(shí)，對(duì)信號(hào)序列的時(shí)域幅頻特性的影響；固定q=8，改變

9、p，使p分別等于8，13，14，觀察參數(shù)p變化對(duì)信號(hào)序列的時(shí)域及幅頻特性的影響，觀察p等于多少時(shí)，會(huì)發(fā)生明顯的泄漏現(xiàn)象，混疊是否也隨之出現(xiàn)？記錄實(shí)驗(yàn)中觀察到的現(xiàn)象，繪出相應(yīng)的時(shí)域序列和幅頻特性曲線。高斯序列定義：function Xa,Fa =gauss(p,q)n=0:15;Xa(n+1)=exp(-(n+1-p).2./q); F=fft(Xa);Fa=abs(F);endclear all; Xa1,Fa1= gauss(8,2);k=0:15;subplot(5,2,1);stem(k,Xa1);hold; plot(k,Xa1);xlabel('n');ylabel(

10、'時(shí)域特性');text(10,0.5,'p=8,q=2');subplot(5,2,2);stem(k,Fa1);hold; plot(k,Fa1);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=2');Xa2,Fa2= gauss(8,4);subplot(5,2,3);stem(k,Xa2); hold; plot(k,Xa2);xlabel('n');ylabel('時(shí)域特性');text(10,0.5,'p=8,q=4'

11、;);subplot(5,2,4);stem(k,Fa2);hold; plot(k,Fa2);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=4');Xa3,Fa3= gauss(8,8);subplot(5,2,5);stem(k,Xa3); hold; plot(k,Xa3);xlabel('n');ylabel('時(shí)域特性');text(10,0.5,'p=8,q=8');subplot(5,2,6);stem(k,Fa3); hold; plot(k,

12、Fa3);xlabel('n');ylabel('幅頻特性');text(8,3,'p=8,q=8');Xa4,Fa4= gauss(13,8);subplot(5,2,7);stem(k,Xa4);hold; plot(k,Xa4);xlabel('n');ylabel('時(shí)域特性');text(10,0.5,'p=13,q=8');subplot(5,2,8);stem(k,Fa4);hold; plot(k,Fa4);xlabel('n');ylabel('幅頻特性&#

13、39;);text(8,3,'p=13,q=8');Xa5,Fa5= gauss(14,8);subplot(5,2,9);stem(k,Xa5);hold; plot(k,Xa5);xlabel('n');ylabel('時(shí)域特性');text(10,0.5,'p=14,q=8');subplot(5,2,10);stem(k,Fa5);hold; plot(k,Fa5);xlabel('n');ylabel('幅頻特性');text(8,3,'p=14,q=8');實(shí)驗(yàn)現(xiàn)象分析：

14、由高斯序列表達(dá)式知n=p為其對(duì)稱軸； 當(dāng)p取固定值時(shí)，時(shí)域序列圖以n=p為其對(duì)稱軸，當(dāng)q取值越來(lái)越大時(shí)，波形變化更加緩和，即向?qū)ΨQ軸兩邊下降越加平緩。q值越大，高頻分量越少，混疊效應(yīng)減弱。由于原序列相當(dāng)于是對(duì)稱截取，且窗口包含序列主要部分，故泄露現(xiàn)象并不明顯。當(dāng)q值固定不變，p變化時(shí)，時(shí)域?qū)ΨQ軸右移，截取窗口不能捕捉信號(hào)主要部分，開(kāi)始無(wú)法代表一個(gè)周期，泄漏現(xiàn)象也越來(lái)越明顯，泄露也產(chǎn)生一定的頻譜混疊。p=14時(shí)的泄漏現(xiàn)象最為明顯，混疊也相應(yīng)產(chǎn)生；(2)、觀察衰減正弦序列xb(n)的時(shí)域和幅頻特性，a=0.1，f=0.0625，檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出幅頻特性曲線，

15、改變f，使f分別等于0.4375和0.5625，觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)位置，有無(wú)混疊和泄漏現(xiàn)象？說(shuō)明產(chǎn)生現(xiàn)象的原因。衰減正弦序列定義：function Xb,Fb = downsin(a,f) n=0:15; Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);F = fft(Xb);Fb=abs(F);clear all;k=0:15; Xb,Fb=downsin(0.1,0.0625);subplot(3,2,1); stem(k,Xb); hold; plot(k,Xb);xlabel('n');ylabel('時(shí)域特性'

16、);text(8,0.5,'a=0.1,f=0.0625');subplot(3,2,2); stem(k,Fb); hold; plot(k,Fb);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.0625');Xb1,Fb1=downsin(0.1,0.4375);subplot(3,2,3); stem(k,Xb1); hold; plot(k,Xb1);xlabel('n');ylabel('時(shí)域特性');text(8,0.5,'a=0

17、.1,f=0.4375');subplot(3,2,4); stem(k,Fb1); hold; plot(k,Fb1);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.4375');Xb2,Fb2=downsin(0.1,0.5625);subplot(3,2,5); stem(k,Xb2); hold; plot(k,Xb2);xlabel('n');ylabel('時(shí)域特性');text(8,0.5,'a=0.1,f=0.5625');s

18、ubplot(3,2,6); stem(k,Fb2); hold; plot(k,Fb2);xlabel('n');ylabel('幅值特性');text(10,3,'a=0.1,f=0.5625');實(shí)驗(yàn)現(xiàn)象分析： 當(dāng)f=f1/fs=1/16=0.0625時(shí)，譜峰位置出現(xiàn)正確。時(shí)間長(zhǎng)度為16Ts，即采樣頻率為16f1，滿足采樣定理。由于截取長(zhǎng)度為周期的整數(shù)倍（一倍），此處不引入頻譜泄漏。但是由于截?cái)嘈?yīng)，相當(dāng)于加一矩形窗，其頻譜中旁瓣導(dǎo)致頻率分量擴(kuò)散，因此有一定程度的頻譜泄漏。增加的部分高頻分量會(huì)造成低頻成分的混疊效應(yīng)，但是由截?cái)鄮?lái)的高頻旁瓣

19、較小，故混疊效應(yīng)較弱。當(dāng)f=f2/fs=7/16=0.4375時(shí)，fs=(16/7)f2，仍滿足采樣定理，譜峰位置正確，且為整數(shù)周期采樣。同上，較弱的頻譜泄漏及混疊效應(yīng)來(lái)自于截?cái)嘈?yīng)。當(dāng)f=f2/fs=9/16=0.5625時(shí)，fs=(16/9)f2，不滿足采樣定理，故產(chǎn)生混疊，此時(shí)譜峰位置不正確?；殳B后的頻譜由圖知與第二種情況頻譜相同，原因是超過(guò)1/2fs的高頻分量疊加到低頻。此時(shí)也產(chǎn)生了一定的頻譜泄漏。 (3)、觀察三角波和反三角波序列的時(shí)域和幅頻特性，用N=8點(diǎn)FFT分析信號(hào)序列xc(n)和xd(n)的幅頻特性，觀察兩者的序列形狀和頻譜曲線有什么異同？繪出兩序列及其幅頻特性曲線。在xc(

20、n)和xd(n)末尾補(bǔ)零，用N=32點(diǎn)FFT分析這兩個(gè)信號(hào)的幅頻特性，觀察幅頻特性發(fā)生了什么變化？?jī)汕闆r的FFT頻譜還有相同之處嗎？這些變化說(shuō)明了什么？程序代碼如下：clear all;n=0:3;k=1:8;Xc(n+1) = n;Xc(n+5) =4-n;Xd(n+1) = 4-n;Xd(n+5) =n;subplot(2,2,1);stem(k-1,Xc);hold; plot(k-1,Xc);xlabel('n');ylabel('時(shí)域特性');text(1,3,'三角波');subplot(2,2,2);stem(k-1,abs(fft

21、(Xc);hold; plot(k-1,abs(fft(Xc);xlabel('k');ylabel('幅頻特性');text(4,10,'三角波');subplot(2,2,3);stem(k-1,Xd);hold; plot(k-1,Xd);xlabel('n');ylabel('時(shí)域特性');text(3,3,'反三角波');subplot(2,2,4);stem(k-1,abs(fft(Xd);hold; plot(k-1,abs(fft(Xd);xlabel('k');yl

22、abel('幅頻特性');text(4,10,'反三角波');Xc(9:32)=0;Xd(9:32)=0;k=1:32;figure;subplot(2,2,1);stem(k-1,Xc); hold; plot(k-1,Xc);xlabel('n');ylabel('時(shí)域特性');text(8,3,'三角波');subplot(2,2,2);stem(k-1,abs(fft(Xc); hold; plot(k-1,abs(fft(Xc);xlabel('k');ylabel('幅頻特性&#

23、39;);text(4,10,'三角波');subplot(2,2,3);stem(k-1,Xd); hold; plot(k-1,Xd);xlabel('n');ylabel('時(shí)域特性');text(2,3,'反三角波');subplot(2,2,4); stem(k-1,abs(fft(Xd); hold; plot(k-1,abs(fft(Xd);xlabel('k');ylabel('幅頻特性');text(4,12,'反三角波');實(shí)驗(yàn)現(xiàn)象分析： 由圖知，三角波序列和反三

24、角波序列的時(shí)域圖像關(guān)于橫軸成鏡像關(guān)系，但頻域圖像完全一樣，只是因?yàn)榉l圖是對(duì)x（k）的值取絕對(duì)值。實(shí)際頻譜相位相反，幅值對(duì)應(yīng)相等。由實(shí)驗(yàn)所得的圖形知，N=32點(diǎn)時(shí)，兩幅度譜都更加密集，更多離散點(diǎn)的幅值顯示，“柵欄效應(yīng)”減小。在原序列的末端填補(bǔ)零值，變動(dòng)了DFT的點(diǎn)數(shù)，人為的改變了對(duì)真實(shí)頻譜采樣的點(diǎn)數(shù)和位置，相當(dāng)于搬動(dòng)了“尖樁柵欄”的位置，從而使得頻譜的峰點(diǎn)和谷點(diǎn)暴露出來(lái)。N=32時(shí)，兩頻譜差別較大，但總體趨勢(shì)仍然都是中間最小，兩側(cè)呈對(duì)稱。     (4)、一個(gè)連續(xù)信號(hào)含兩個(gè)頻率分量，經(jīng)采樣得x(n)=sin2*0.125n+cos2*(0.125+f)n n=0,

25、1,N-1已知N=16，f分別為1/16和1/64，觀察其頻譜；當(dāng)N=128時(shí)，f不變，其結(jié)果有何不同，為什么？clear all;N=16;detf=1/16;n=0:N-1;x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);subplot(2,2,1);stem(n,x1); xlabel('n');ylabel('時(shí)域特性');text(6,1,'N=16,detf

26、=1/16');subplot(2,2,2);stem(n,abs(fft(x1);xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16');subplot(2,2,3);stem(n,x2);xlabel('n');ylabel('時(shí)域特性');text(6,1,'N=16,detf=1/64');subplot(2,2,4);stem(n,abs(fft(x2);xlabel('n');ylabel('幅值特性&#

27、39;);text(6,4,'N=16,detf=1/64');N=128;detf=1/16;n=0:N-1;x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);detf = 1/64;x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);figure;subplot(2,2,1);stem(n,x3);xlabel('n');ylabel('時(shí)域特性');axis(0 128 -2 2);text(6,1.5,'N=128,det

28、f=1/16');subplot(2,2,2);stem(n,abs(fft(x3);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');subplot(2,2,3);stem(n,x3);xlabel('n');ylabel('時(shí)域特性');axis(0 128 -2 2);text(6,1.5,'N=128,detf=1/16');subplot(2,2,4);stem(n,abs(fft

29、(x4);xlabel('n');ylabel('幅值特性');axis(0 128 -10 70);text(40,60,'N=128,detf=1/16');實(shí)驗(yàn)現(xiàn)象分析： 由圖可以看出N=16時(shí)，當(dāng)f由1/16減小為1/64時(shí)，頻譜圖發(fā)生較大變化。此時(shí)分辨出1/64對(duì)頻率的分辨率要求提高，需要通過(guò)增大信號(hào)序列長(zhǎng)度以提高頻率的分辨率。f由1/16減小為1/64時(shí)，此時(shí)采樣是非周期采樣，DFT后導(dǎo)致頻譜泄漏，柵欄效應(yīng)使得無(wú)法顯現(xiàn)出0.125+1/64和0.125兩個(gè)頻率譜峰。   當(dāng)N增加至128時(shí)，頻譜更加密集，分辨率提高

30、（1/128），能夠分辨出兩譜峰。(5)、用FFT分別實(shí)現(xiàn)xa(n)（p8，q2）和 xb(n)（a0.1，f0.0625）的16點(diǎn)圓周卷積和線性卷積。clear all;N=16;n=0:N-1;p=8;q=2;Xa(n+1)=exp(-(n-p).2./q);a=0.1;f=0.0625;Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n);Fa=fft(Xa); Fb=fft(Xb);Fx=Fa.*Fb;X51=ifft(Fx);stem(n,X51);Xa(N+1:2*N-1)=0;Xb(N+1:2*N-1)=0;Fa=fft(Xa); Fb=fft(Xb);Fc=Fa

31、.*Fb;X52=ifft(Fc);figure;stem(1:2*N-1,X52);（6） 產(chǎn)生一512點(diǎn)的隨機(jī)序列xe(n)，并用xc(n)和xe(n)作線性卷積，觀察卷積前后xe(n)頻譜的變化。要求將xe(n)分成8段，分別采用重疊保留法和重疊相加法。 重疊保留法：1 在長(zhǎng)序列xe(n)前補(bǔ)N1-1=7個(gè)零，由于最后輸出y(n)要多取8個(gè)，為滿足L>=N1+N2-1,在xe(n)末尾補(bǔ)8個(gè)零，構(gòu)成長(zhǎng)度為527的序列；2 分段時(shí)從xe(n)中每次取72個(gè)數(shù)，即要做72點(diǎn)的DFT，但每次取數(shù)的起始點(diǎn)為每隔65個(gè)開(kāi)始取一次，即相鄰兩段之間有7個(gè)數(shù)是相重疊的，也就是將要去掉的部分；3 將

32、短序列xc(n)作72點(diǎn)的DFT，與每段序列做圓周卷積；4 去掉所得每段卷積序列的前N1-1=7個(gè)點(diǎn)，取其后的65個(gè)點(diǎn)。將各相鄰段收尾連結(jié)起來(lái)，構(gòu)成長(zhǎng)為519 的線性卷積序列。程序代碼：n=1:512;xe1=randn(1,512);n=1:7; xe(n)=0;n=8:519; xe(n)=xe1(n-7);n=520:527; xe(n)=0;n=1:4; xc(n)=n-1;n=5:8; xc(n)=9-n;n=0:526subplot(2,1,1);stem(n,xe);ylabel('xe(k)');title('512點(diǎn)的隨機(jī)序列');subplo

33、t(2,1,2);n=0:7;stem(n,xc); y1=conv(xe1,xc);figure; subplot(2,1,1); n=0:518; stem(n,y1); ylabel('y1(n)');title('線性卷積');hxc=fft(xc,72);y=zeros(1,527);for i=0:7 n=(65*i+1):65*(i+1)+7; xi(n-65*i)=xe(n); hxi=fft(xi,72); yi=ifft(hxc.*hxi); k=(65*i+1):65*(i+1); y(k)=yi(k-65*i+7);end subplot(2,1,2);m=0:526;stem(m,y);ylabel('y(n)');t