數(shù)字信號處理實驗指導(dǎo)書

1、數(shù)字信號處理實驗指導(dǎo)書2實驗一 離散時間與系統(tǒng)的傅立葉分析一、實驗?zāi)康挠酶盗⑷~變換對信號和系統(tǒng)進行頻域分析。二、實驗原理對信號進行頻域分析就是對信號進行傅立葉變換。對系統(tǒng)進行頻域分析即對它的單位脈沖響應(yīng)進行傅立葉變換，得到系統(tǒng)的傳輸函數(shù)。也可以由差分方程經(jīng)過；傅立葉變換直接 求它的傳輸函數(shù)。傳輸函數(shù)代表的就是系統(tǒng)的頻率響應(yīng)特性。但傳輸函數(shù)是的連續(xù)函數(shù)，計算機只能計算出有限個離散頻率點的傳輸函數(shù)值，因此得到傳輸函數(shù)以后，應(yīng)該在 02 之間取許多點，計算這些點的傳輸函數(shù)的值，并取它們的包絡(luò)，該包絡(luò)才是需要的頻率特性。當然，點數(shù)取得多一些，該包絡(luò)才能接近真正的頻率特性。 注意：非周期信號的頻率特性是

2、的連續(xù)函數(shù)，而周期信號的頻率特性是離散譜，它們的計算公式不一樣，響應(yīng)的波形也不一樣。 三、實驗內(nèi)容 1已知系統(tǒng)用下面差分方程描述： y(n)=x(n)十a(chǎn)y(n一1) 試在a095和a一0.5兩種情況下用傅立葉變換分析系統(tǒng)的頻率特性。要求寫出系統(tǒng)的傳輸函數(shù)和幅度響應(yīng)，并打印|H(ejw)|曲線。 2已知兩系統(tǒng)分別用下面差分方程描述： y1(n)x(n)+x(n一1) y2(n)x(n)一x(n一1) 試分別寫出它們的傳輸函數(shù)和幅度響應(yīng)，并分別打印|H(ejw)|曲線。 3已知信號x(n)R3(n)，試分析它的頻域特性，要求打印|H(ejw)|曲線。 4假設(shè)x(n)a(n)，將x(n)以2為周期

3、進行周期延拓，得到x(n)，試分析它的頻率特性，并畫出它的幅頻特性。 四、實驗用MATLAB函數(shù)介紹 1abs 功能：求絕對值(復(fù)數(shù)的模)。 yabs(x)：計算實數(shù)x的絕對值。當x為復(fù)數(shù)時得到x的模(幅度值)。 當x為向量時，計算其每個元素的模，返回模向量y。 2angle 功能：求相角。 Ph=angle(x)：計算復(fù)向量x的相角(rad)。Ph值介于 - 和+之間. 3freqz： 計算數(shù)字濾波器H(z)的頻率響應(yīng)。 Hfreqz(B，A，w)：計算由向量w指定的數(shù)字頻率點上數(shù)字濾波器H(z)的頻率響應(yīng)H(ejw)，結(jié)果存于H向量中.向量B和A分別為數(shù)字濾波器系統(tǒng)函數(shù)H(z)的分子和分母

4、多項式系數(shù)。 H，wfreqz(B，A，M，whole)：計算出M個頻率點上酌頻率響應(yīng)，存放在，H向量中，M個頻率存放在向量w中。廳reqz函數(shù)自動將這M個頻點均勻設(shè)置在頻率范圍 0,2 上。缺省whole時，M個頻點均勻設(shè)置在頻率范圍（0，)上。調(diào)用參數(shù)B和A與系統(tǒng)函數(shù)的關(guān)系由下式給出： 缺省w和M時，freqz自動選取512個頻率點計算。不帶輸出向量的freqz函數(shù)將自動繪出幅頻和相頻曲線。其它幾種調(diào)用格式可用命令help查閱。Freqz函數(shù)用于計算模擬濾波器的頻率響應(yīng)函數(shù)，詳細功能及調(diào)用格式用help查看。例如，八階梳狀濾波器系統(tǒng)函數(shù)為 H(z) =B(z)=1-z -8 用下面的簡單程

5、序繪出 H(z)的幅頻和相頻特性曲線如圖10.2.1 所示。 %Example for Freqz B=1 0 0 0 0 0 0 0 -1; A=1; freqz(B,A)五、實驗報告要求1對各實驗內(nèi)容進行理論分析和推導(dǎo)。2分析各實驗內(nèi)容，并和理論分析推導(dǎo)結(jié)果進行對比 。13實驗二 快速變換及其應(yīng)用一、實驗?zāi)康?1)在理論學(xué)習的基礎(chǔ)上。通過本實驗，加深對FFT的理解，熟悉FFT子程序；(2)熟悉應(yīng)用FFT對典型信號進行頻譜分析的方法；(3)了解應(yīng)用FFT進行信號頻譜分析過程中可能出現(xiàn)的問題，以便在實際中正確應(yīng)用FFT；(4)熟悉應(yīng)用FFT實現(xiàn)兩個序列的線性卷積的方法；(5)對DCT變換用作數(shù)

6、據(jù)壓縮有初步的認識。二、實驗內(nèi)容實驗中用到的信號序列a)余弦序列 b)衰減正弦序列c)三角波序列d)反三角波序列上機實驗內(nèi)容：(1) xn=1, 2, 3, 4, 5, hn=-3, -2, -1, 0, 1, 2, 3, 4, 5, 6. 借助FFT求它們的線性卷積。(2)觀察余弦序列的時域和幅頻特性，了解當長度N取不同值時，對信號序列的時域和幅頻特性的影響。（i）信號頻率f=500Hz。采樣頻率fs=1500 Hz。?。╝）N=16 時；（b）N=160時；（c）N=1600時。觀察出現(xiàn)的泄露現(xiàn)象?；殳B是否同時出現(xiàn)？繪出泄露現(xiàn)象最嚴重時時域序列和幅頻特性曲線。（ii）信號頻率f=500Hz

7、。N=1600。取采樣頻率（a）fs=500Hz時；（b）fs=1000Hz時；（c）fs=2000Hz時。觀察出現(xiàn)的混疊現(xiàn)象。繪出泄露現(xiàn)象最嚴重時時域序列和幅頻特性曲線。 (3)觀察衰減正弦序列Xb(n)的時域和幅頻特性，a=0.1,f=0.0625,檢查譜峰出現(xiàn)位置是否正確，注意頻譜的形狀，繪出時域序列和幅頻特性曲線。改變f，使f分別等于0.4375和0.5625,觀察這兩種情況下，頻譜的形狀和譜峰出現(xiàn)位置，有無混疊和泄露現(xiàn)象？說明產(chǎn)生現(xiàn)象的原因。(4)觀察三角波和反三角波序列的時域和幅頻特性，用N=8點FFT分析信號序列Xc(n)和Xd(n)的幅頻特性，觀察兩者的序列形狀和頻譜曲線有什么

8、異同？繪出兩時域序列及其幅頻特性曲線。在Xc(n)和Xd(n)末尾補零，用N=16點FFT分析這兩個信號的幅頻特性，觀察幅頻特性發(fā)生了什么變化？兩情況的FFT頻譜還有相同之處嗎？這些變化說明了什么？(5)一個連續(xù)信號含兩個頻率分量，經(jīng)采樣得已知N=16，f分別為1/16，1/64，觀察其頻譜；當N=128時，f不變，其結(jié)果有何不同，為什么？四、實驗報告要求(1)簡述實驗?zāi)康募霸怼?2)分析得到的波形，說明參數(shù)改變對時域和頻域的影響，按實驗要求附上實驗信號序列和幅頻特性曲線，(3)總結(jié)實驗中的主要結(jié)論實驗三 IIR數(shù)字濾波器的設(shè)計一、實驗?zāi)康模?）掌握雙線性變換法及脈沖響應(yīng)不變法設(shè)計IIR數(shù)字

9、濾波器的具體設(shè)計方法及其原理，熟悉用雙線性變換法及脈沖響應(yīng)設(shè)計低通、高通和帶通IIR數(shù)字濾波器的計算機編程。（2）觀察雙線性變換及脈沖響應(yīng)不變法設(shè)計的濾波器的頻率特性，了解雙線性變換法及脈沖響應(yīng)不變法的特點。（3）熟悉Butterworth濾波器、Chebyshew濾波器和橢圓濾波器的頻率特性。二、實驗內(nèi)容實驗中有關(guān)變量的定義：fp 通帶邊界頻率fr 阻帶邊界頻率Ap 通帶波動Ar 最小阻帶衰減fs 采樣頻率T 采樣周期（1）fp=0.3kHz,Ap=0.8dB,fr=0.2kHz,Ar=20dB,T=1ms。設(shè)計一Chebyshew高通濾波器；局部放大圖形，觀察其通帶損耗和阻帶衰減是否滿足要

10、求。記錄fp=0.3kHz和fr=0.2kHz處的通帶波動和阻帶衰減，繪制所設(shè)計濾波器的圖形。（2）fp=0.2kHz,Ap=1dB,fr=0.3kHz,Ar=25dB,T=1ms。用雙線性變換法設(shè)計一Butterworth數(shù)字低通濾波器，觀察其通帶損耗和阻帶衰減是否滿足要求。記錄fp=0.3kHz和fr=0.2kHz處的通帶波動和阻帶衰減，繪制所設(shè)計濾波器的圖形。（3）利用雙線性變換法設(shè)計一Butterworth型數(shù)字帶通濾波器，已知fs=30kHz,其等效的模擬濾波器指標為Ap<3dB,2kHz<f3kHz; Ar 5 dB,f 6kHz; Ar 20dB,f 1.5kHz。觀

11、察其通帶損耗和阻帶衰減是否滿足要求。記錄通帶波動和阻帶衰減，繪制所設(shè)計濾波器的圖形。三、思考題（1）雙線性變換法中和之間的關(guān)系是非線性的，在實驗中你注意到這樣非線性關(guān)系了嗎？從哪幾種數(shù)字濾波器的幅頻特性曲線中可以觀察到這幾種非線性關(guān)系？四、實驗報告要求（1）簡述實驗?zāi)康募霸?。?）按實驗要求附上所設(shè)計的濾波器傳輸函數(shù)及相應(yīng)的幅頻特性曲線，定性分析他們的性能，判斷設(shè)計是否滿足要求。（3）簡要回答思考題。實驗四 FIR數(shù)字濾波器的設(shè)計一、實驗?zāi)康模?）掌握用窗函數(shù)法設(shè)計FIR濾波器的原理及方法，熟悉相應(yīng)的計算機編程。（2）熟悉線性相位FIR濾波器的幅頻特性和相頻特性。二、實驗內(nèi)容（1）給定濾波器

12、設(shè)計指標通帶邊界頻率 0.25阻帶邊界頻率 0.35通帶波動 0.03db最小阻帶衰減 50db 海明窗選擇固定窗自行設(shè)計MATLAB程序，得到滿足設(shè)計指標的濾波器。（2）給定濾波器設(shè)計指標通帶邊界頻率 0.5阻帶邊界頻率 0.7通帶波動 0.01db最小阻帶衰減 80db用Kaiser窗自行設(shè)計MATLAB程序，得到滿足設(shè)計指標的濾波器。（3）繪制所設(shè)計濾波器的幅度和相位圖形三、實驗報告要求（1）簡述實驗?zāi)康募霸怼＃?）按實驗要求附上所設(shè)計的濾波器的MATLAB程序。（3）繪出相應(yīng)的幅頻特性曲線（幅度和相位）。用海明窗設(shè)計一個FIR濾波器，其中Wp=0.2*pi，Ws=0.3*pi,通帶衰

13、減不大于0.25dB,阻帶衰減不小于50dB。代碼：wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1;n=0:M-1;wc=(ws+wp)/2;%ideal LPF cutoff frequencyhd = ideal_lp(wc,M);w_ham=(hamming(M)'h=hd.*w_ham;% .*db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1:1:501) %Min stopband attenuati

14、on%Plotssubplot(2,2,1);stem(n,hd);title('Ideal Impulse Response')axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('hd(n)');%subplot(2,2,2);stem(n,w_ham);title('Hamming Window')axis(0 M-1 0 1.1);xlabel('n');ylabel('w(n)');%subplot(2,2,3);stem(n,h);title('Act

15、ual Impulse Response')axis(0 M-1 -0.1 0.3);xlabel('n');ylabel('h(n)');%實驗一離散時間與系統(tǒng)的頻譜分析%1. y(n)=x(n)+ay(n-1)%1)a=0.95B=1A=1, -0.95AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)

16、plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%a=-0.5B=1；A=1, 0.5；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|&

17、#39;)%2.y1(n)=x(n)+x(n-1), y2(n)=x(n)-x(n-1)%1)y1(n)=x(n)+x(n-1)B=1 1；A=1；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%y2(n)=x

18、(n)-x(n-1)B=1 -1；A=1；AM,w=freqz(B,A,'whole');am=abs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|')%x(n)=R3(N)B=1 1 1 1；A=1；AM,w=freqz(B,A,'whole');am=a

19、bs(AM);an=angle(AM);subplot(2,1,1)plot(w/pi,am)xlabel('omega/pi')ylabel('|H|')grid onsubplot(2,1,2)plot(w/pi,an)grid onxlabel('omega/pi')ylabel('|H|') 實驗二 快速變換及其應(yīng)用%（1）xn=1:5hn=-3:6xnd=length(xn)hnd=length(hn)ynd=xnd+hnd-1xnL=xn zeros(1,ynd-xnd)hnL=hn zeros(1,ynd-hnd)X

20、L=fft(xnL)HL=fft(hnL)YL=XL.*HLyn=ifft(YL)%conv(xn,hn);直接計算%(2)n=0:1600;omega=2*pi*500/1500x=sin(omega*n)y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%(3)n=0:15;a=input('a=');f=input('

21、f=');x=exp(-(a*n).*sin(2*pi*f*n);y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度序列')subplot(3,1,3)plot(n,angle(y);title('相位特性')%(4)%1) N=8; n=0:7; x=0,1,2,3,4,3,2,1;y=fft(x);subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,ab

22、s(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%2)N=8;n=0:7;x=4,3,2,1,0,1,2,3;y=fft(x);subplot(3,1,1)stem(n,x)subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%3) N=16; n=0:15; x=0,1,2,3,4,3,2,1,0,0,0,0,0,0,0,0;y=fft(x)

23、;subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')%4) N=16; n=0:15; x=4,3,2,1,0,1,2,3,0,0,0,0,0,0,0,0;y=fft(x);subplot(3,1,1)stem(n,x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')sub

24、plot(3,1,3)plot(n,angle(y);%(4)N=input('N=');n=0:N-1;df=input('df=');x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+df)*n);y=fft(x);subplot(3,1,1)stem(x);title('序列')subplot(3,1,2)plot(n,abs(y);title('幅度特性')subplot(3,1,3)plot(n,angle(y);title('相位特性')實驗三 IIR濾波器設(shè)計%(1)fp=0.3

25、k,Ap=0.8db;fr=0.2k,Ar=20db;T=1%ms.T=1wp=2*0.3*T;wr=2*0.2*T;Ap=0.8;Ar=20;N,wc=cheb1ord(wp,wr,Ap,Ar) b,a=cheby1(N,Ap,wc,'high') %濾波器差分方程：%y(n)=b0x(n)+ b1x(n-1) + b2x(n-2)+%+a1y(n-1)+ a2y(n-2)+sys1=tf(b,a,-1)% w=0:pi/512:pi;H, w=freqz(b,a,w);fg=20*log10(abs(H);plot(w/pi,fg); title('幅度特性'

26、;)grid%(2)fp=0.2k,Ap=1db,fr=0.3k,Ar=25db,T=1m%s%雙線性變換法wp=2*0.2*1;wr=2*0.3*1;T=0.001Ap=1;Ar=25;N,wc=buttord(wp,wr,Ap,Ar,'s');z,p,k=buttap(N)WC=(2/T)*tan(wc*pi/2)num=kn=N:-1:0den=real(poly(p)den=den./(WC.n);sys=tf(num,den)b,a=bilinear(num,den,1/T)%濾波器差分方程：%y(n)=b0x(n)+ b1x(n-1) + b2x(n-2)+%+a1y(n-1)+ a2y(n-2)+% m=linspace(0,pi,512);H,w=freqz(b,a,m);H2=20*log10(abs(H);W2=angle(H);subplot(2,1,1)hold onplot(w/pi,H2,'r');grid onsubplot(2,1,2)hold onplot(w/pi,W2,'r');grid on%(3)banpass:fs=30k;Ap<3db,2k<=f<=3k;Ar>=5db,f>=6k;Ar>=20db,f<=1.5k.wp1=2*2