河海大學-通信工程-DSP實驗2

1、實驗二用 FFT 對信號進行頻譜分析班級：姓名：學號：一、實驗目的1、加深對離散信號DFT的理解；2、掌握 FFT 算法的流程及其MATLAB實現；3、利用 FFT 對典型信號進行頻譜分析；4、結合理論知識，對頻譜分析中出現的有關現象進行理論分析。二、實驗原理采樣序列的頻譜是被采樣模擬信號頻譜的周期延拓，當采樣頻率不滿足奈奎斯特采樣定理時， 就會發(fā)生頻譜的混疊。解決混疊問題的唯一辦法是保證采樣頻率足夠高，以防止頻譜混疊。對一個時間無限的信號，雖然頻帶有限，但在實際DFT 運算中，時間長度總是取有限值，在將信號截短的過程中，出現了分散的擴展譜線的現象，稱為頻譜泄漏或功率泄漏。泄漏也會引起混疊。從

2、某種意義上來講，用DFT 來觀察頻譜就如同通過一個柵欄來觀看景象一樣，只能在離散點上看到真實的頻譜，這樣一些頻譜的峰點或谷點就可能被“尖樁的柵欄”擋住，也就是正好落在兩個離散采樣點之間，不能被觀察到。 減小柵欄效應的一個方法是在原序列的末端填補一些零值，從而變動DFT的點數。三、實驗內容1、對高斯序列 xa(n)，令 p=8， q 分別等于 2、 4 和 8，觀察 q 值的改變對高斯序列時域特性和幅頻特性的影響，并給出理論解釋：(np )2xa (n)eq0 n150 其他2、對衰減正弦序列xb(n)， a=0.1 ， f=0.0625 ，觀察其時域和幅頻特性，檢查譜峰位置是否正確；改變 f，

3、使 f 分別等于 0.4375和 0.5625 ，比較這兩種情況下頻譜形狀和譜峰位置，并給出理論解釋：e an sin( 2 fn)0 n 15xb (n)其他03、觀察三角波序列xc(n)和反三角波序列xd(n)的時域特性和幅頻特性，用8點 FFT分析幅頻特性， 觀察兩者的序列形狀和頻譜曲線有什么異同？為什么？在xc(n)和x (n)d的末尾補零，用32 點 FFT分析幅頻特性，觀察幅頻特性發(fā)生了什么變化？這些變化說明了什么？n0n3xc (n)8n4n70其他4n0n3xd (n )n44n70其他4、對連續(xù)單頻周期信號 sin( 2fat ) ，按采樣頻率 f s8 fa 進行采樣，截取

4、長度分別選N=16和 N 20，比較其幅頻特性，并給出理論解釋；當截取長度選N 164 時，與 N 20比較，頻譜發(fā)生了哪些變化？為什么？5、一個連續(xù)信號含兩個頻率分量，經采樣得：x( n)sin( 20.125n)cos(2(0.125f ) n) ,n0,1, N1令 N=16， f 分別為 1/16 和 1/64 ，觀察其頻譜，兩個譜峰是否可以區(qū)分開？為什么？當 N=128 時， f 不變，結果有何不同？為什么？四、 實驗結果與分析1、實驗結果如下1p=8 q=2 時 的 高 斯 序 列4p=8 q=2 時 的 幅 頻 特 性)n0.5k2(xX05101505101500t/Tk1p=

5、8 q=4 時 的 高 斯 序 列4p=8 q=4 時 的 幅 頻 特 性)n0.5k2(xX05101505101500t/Tkp=8 q=8 時 的 高 斯 序 列p=8 q=4 時 的 幅 頻 特 性110)n0.5k5(xX05101505101500t/Tk分析：由高斯序列表達式n=p 為其對稱軸，當p 取固定值時，時域圖都關于n=8 對稱截取長度為周期的整數倍，沒有發(fā)生明顯的泄露現象，但存在混疊， 當 q 由 2 增加至 8 的過程中，時域圖形變化越來越平緩，中間包絡越來越大，可能函數周期開始增加，頻率降低，漸漸小于 fs/2，混疊減弱。2、實驗結果如下：)n(x)n(x)n(x1

6、f=0.0625 時 域 圖)k0(X-1510150t/Tf=0.4375 時 域 圖1)k0(X-1051015t/Tf=0.5625 時 域 圖1)k0(X-1051015t/Tf=0.0625 幅 頻 特 性 圖420051015kf=0.4375 幅 頻 特 性 圖420051015kf=0.5625 幅 頻 特 性 圖420051015k分析：當 f= 0.0625 時，譜峰位置出現正確，存在混疊現象，時域采樣為一周期，不滿足采樣定理。當 f=0.4375 和 0.5625 時，時域圖像關于 Y 軸對稱，頻域完全相同。這是因為頻域圖是取絕對值的結果， 所以完全相同。 另外由于時域采

7、樣為 6 個半周期， 滿足采樣定理， 無混疊；但由于截取長度不是周期的整數倍，出現泄漏。3、 .實驗結果如下：N=8 時，)n(cx)n(dx時域特性圖幅頻特性圖420315)k2(10cX15024680246800nk4時域特性圖20幅頻特性圖315)k2(10dX15024680246800nk分析：由圖知， 三角波序列和反三角波序列的時域圖像成鏡像關系， 但頻域圖像完全一樣， 只是因為幅頻圖是對 X(k）的值取絕對值。N=32 時，)n(cx)n(dx時域特性圖幅頻特性圖420315)k2(10cX1501020304001020304000nk4時域特性圖20幅頻特性圖315)k2(

8、10dX1501020304001020304000nk分析：由圖知， N=32 點時， Xc(n)和 Xd(n)的幅頻特性都更加密集，更多離散點的幅值顯示， “柵欄效應”減小，分辨率提高，而對于Xd(n)來說變化更加明顯。在原序列的末端填補零值，變動了DFT 的點數，人為的改變了對真實頻譜采樣的點數和位置，相當于搬動了“尖樁柵欄”的位置， 從而使得頻譜的峰點和谷點暴露出來。N=32 時，Xc(n)和 Xd(n)的頻譜差別較大，但總體趨勢仍然都是中間的最小，兩側呈對稱。4、實驗結果如下：)n(x)n(x)n(x10-1010-1010-10N=16時域圖)k(X51015t/TN=20時域圖)

9、k(X5101520 t/TN=164 時 域 圖)k(X50100 150 200 t/TN=16幅頻特性圖1050051015kN=20幅頻特性圖105005101520kN=164幅頻特性圖100500050 100 150 200 k分析：因為所選取的采樣頻率fs=8fa，8>2，滿足奈奎斯特采樣頻率，由此可以得到單一譜線的 DFT 結果， N=16 與 N=20 之間的頻譜差距較大， 是因為 N=20 所選的時間長度不是周期的整數倍。而當 N=164 時，其圖形在形狀上與 N=20 時差不多，但是在 X（k）達最值時， k 的取值不同。5、實驗結果如下：8N=16f=1/16

10、頻 譜 圖8N=16f=1/64 頻 譜 圖66)kk(4(4XX2205101505101500kk80N=128f=1/16 頻 譜 圖80N=128f=1/64 頻 譜 圖6060)kk(40(40XX202005010015005010015000kk分析：由圖可以看出 N=16 時，當 f 由 1/16 減小為 1/64 時，頻譜圖出現失真，可能是 f 的改變引起的周期變化導致混疊。當 N 增加至 128 時，頻譜更加密集，分辨率明顯提高，混疊現象消失。五、實驗小結此次實驗主要研究利用DFT做連續(xù)頻譜的分析，并驗證了其中的混疊、泄漏、柵欄效應等現象的產生、 減小措施等，這次的實驗明顯

11、比第一次要難， 做起來也沒有那么順手， 但是經過一遍遍的看、改，也總算是勉強完成了，準確率不知道怎么樣，但收獲卻是頗豐，實驗將課本知識與實踐中的現象聯合在了一起，顯得更加形象。程序：1、n=0:1:15; p=8; q1=2; xa1=exp(-(n-p).*(n-p)/q1);subplot(3,2,1); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('x(n)');title('p=8 q=2 時的高斯序列 '); xk1=abs(fft(xa1);subplot(3,2,2); stem(

12、n,xk1); xlabel('k'); ylabel('X(k)');title('p=8 q=2 時的幅頻特性 ');n=0:1:15; p=8; q2=4;xa2=exp(-(n-p).*(n-p)/q2);subplot(3,2,3); plot(n,xa2,'-*'); xlabel('t/T'); ylabel('x(n)');title('p=8 q=4 時的高斯序列 '); xk2=abs(fft(xa2);subplot(3,2,4); stem(n,xk2); x

13、label('k'); ylabel('X(k)');title('p=8 q=4 時的幅頻特性 ');n=0:1:15; p=8; q3=8; xa3=exp(-(n-p).*(n-p)/q3);subplot(3,2,5); plot(n,xa3,'-*'); xlabel('t/T'); ylabel('x(n)');title('p=8 q=8 時的高斯序列 '); xk3=abs(fft(xa3);subplot(3,2,6); stem(n,xk3); xlabel(&#

14、39;k'); ylabel('X(k)');title('p=8 q=4 時的幅頻特性 ');2、n=0:1:15;a=0.1;f=0.0625;xb=exp(-a*n).*sin(2*pi*f*n);subplot(3,2,1);plot(n,xb,'-*');xlabel('t/T');ylabel('x(n)');title('f=0.0625 時域圖 ');xk=abs(fft(xb);subplot(3,2,2);stem(n,xk);xlabel('k');yl

15、abel('X(k)');title('f=0.0625 幅頻特性圖 ');n=0:1:15;a=0.1;f1=0.4375;xb1=exp(-a*n).*sin(2*pi*f1*n);subplot(3,2,3);plot(n,xb1,'-*');xlabel('t/T');ylabel('x(n)');title('f=0.4375 時域圖 ');xk1=abs(fft(xb1);subplot(3,2,4);stem(n,xk1);xlabel('k');ylabel('

16、;X(k)');title('f=0.4375 幅頻特性圖 ');n=0:1:15;a=0.1;f2=0.5625;xb2=exp(-a*n).*sin(2*pi*f2*n);subplot(3,2,5);plot(n,xb2,'-*');xlabel('t/T');ylabel('x(n)');title('f=0.5625 時域圖 ');xk2=abs(fft(xb2);subplot(3,2,6);stem(n,xk2);xlabel('k');ylabel('X(k)'

17、;);title('f=0.5625 幅頻特性圖 ');3、n1=0:1:7; xc1=0 1 2 3 4 3 2 1;subplot(2,2,1); plot(n1,xc1,'-*'); xlabel('n'); ylabel('xc(n)');title(' 時域特性圖 '); xk1=abs(fft(xc1);subplot(2,2,2); stem(n1,xk1); xlabel('k'); ylabel('Xc(k)');title(' 幅頻特性圖 ');n

18、2=0:1:7; xd1=4 3 2 1 0 1 2 3;subplot(2,2,3); plot(n2,xd1,'-*'); xlabel('n'); ylabel('xd(n)');title(' 時域特性圖 '); xk2=abs(fft(xd1);subplot(2,2,4); stem(n2,xk2); xlabel('k'); ylabel('Xd(k)');title(' 幅頻特性圖 ');n1=0:1:31; xc1=0 1 2 3 4 3 2 1 zeros(1,2

19、4);subplot(2,2,1); plot(n1,xc1,'-*'); xlabel('n'); ylabel('xc(n)');title(' 時域特性圖 '); xk1=abs(fft(xc1);subplot(2,2,2); stem(n1,xk1); xlabel('k'); ylabel('Xc(k)');title(' 頻域特性圖 ');n2=0:1:31; xd1=4 3 2 1 0 1 2 3 zeros(1,24);subplot(2,2,3); plot(n2

20、,xd1,'-*'); xlabel('n'); ylabel('xd(n)');title(' 時域特性圖 '); xk2=abs(fft(xd1);subplot(2,2,4); stem(n2,xk2); xlabel('k'); ylabel('Xd(k)');title(' 幅頻特性圖 ');4、n1=0:1:15;xa1=sin(2*pi*n1 /8);subplot(3,2,1); plot(n1,xa1,'-*'); xlabel('t/T&#

21、39;); ylabel('x(n)');title('N=16 時域圖 '); xk1=abs(fft(xa1);subplot(3,2,2); stem(n1,xk1); xlabel('k'); ylabel('X(k)');title('N=16 幅頻特性圖 ');n2=0:1:19;xa2=sin(2*pi*n2 /8);subplot(3,2,3); plot(n2,xa2,'-*'); xlabel('t/T'); ylabel('x(n)');titl

22、e('N=20 時域圖 '); xk2=abs(fft(xa2);subplot(3,2,4); stem(n2,xk2); xlabel('k'); ylabel('X(k)');title('N=20 幅頻特性圖 ');n3=0:1:163;xa3=sin(2*pi*n3 /8);subplot(3,2,5); plot(n3,xa3,'-*'); xlabel('t/T'); ylabel('x(n)');title('N=164 時域圖 '); xk3=abs(fft(xa3);subplot(3,2,6); stem(n3,xk3); xlabel('k'); ylabel('X(k)');title('N=164 幅頻特性圖 ');5、n1=0:1:15;deltaf1=1 /16;xa1=sin(2*pi*0.125*n1)+cos(2*pi*(0.125+deltaf1)*n1); xk