數(shù)字信號處理實(shí)驗(yàn)例題.doc

《數(shù)字信號處理實(shí)驗(yàn)例題.doc》由會員分享，可在線閱讀，更多相關(guān)《數(shù)字信號處理實(shí)驗(yàn)例題.doc（27頁珍藏版）》請在裝配圖網(wǎng)上搜索。

實(shí)驗(yàn)一 離散時間信號與系統(tǒng) 例1-1 用MATLAB計算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的離散卷積。 解 MATLAB程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel(n); ylabel(幅度); 圖1.1給出了卷積結(jié)果的圖形，求得的結(jié)果存放在數(shù)組c中為：{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB計算差分方程 當(dāng)輸入序列為 時的輸出結(jié)果 。 解 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel(n);ylabel(幅度) 圖 1.2 給出了該差分方程的前41個樣點(diǎn)的輸出，即該系統(tǒng)的單位脈沖響應(yīng)。 　 　 例1-3 用MATLAB計算例1-2差分方程 所對應(yīng)的系統(tǒng)函數(shù)的DTFT。 解 例1-2差分方程所對應(yīng)的系統(tǒng)函數(shù)為： 其DTFT為 用MATLAB計算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title(實(shí)部) xlabel(\omega/\pi);ylabel(幅度) subplot(2,2,2); plot(w/pi,imag(h));grid title(虛部) xlabel(\omega/\pi);ylabel(Amplitude) subplot(2,2,3); plot(w/pi,abs(h));grid title(幅度譜) xlabel(\omega/\pi);ylabel(幅值) subplot(2,2,4); plot(w/pi,angle(h));grid title(相位譜) xlabel(\omega/\pi);ylabel(弧度) 　 　 實(shí)驗(yàn)二 離散傅里葉變換及其快速算法 例2-1 對連續(xù)的單一頻率周期信號 按采樣頻率 采樣，截取長度N分別選N =20和N =16，觀察其DFT結(jié)果的幅度譜。 解 此時離散序列 ，即k=8。用MATLAB計算并作圖，函數(shù)fft用于計算離散傅里葉變換DFT，程序如下： k=8; n1=[0:1:19]; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel(t/T);ylabel(x(n)); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel(k);ylabel(X(k)); n2=[0:1:15]; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel(t/T);ylabel(x(n)); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel(k);ylabel(X(k)); 計算結(jié)果示于圖2.1，(a)和(b)分別是N=20時的截取信號和DFT結(jié)果，由于截取了兩個半周期，頻譜出現(xiàn)泄漏；(c) 和(d) 分別是N=16時的截取信號和DFT結(jié)果，由于截取了兩個整周期，得到單一譜線的頻譜。上述頻譜的誤差主要是由于時域中對信號的非整周期截斷產(chǎn)生的頻譜泄漏。 例2-2 用FFT計算兩個序列 的互相關(guān)函數(shù) 。 解 用MATLAB計算程序如下： 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(幅度); 其計算結(jié)果如圖2.2所示。 　 　 例2-3計算兩個序列的的互相關(guān)函數(shù)，其中 x(n)={2 3 5 2 1 –1 0 0 12 3 5 3 0 –1 –2 0 1 2} y(n)=x(n-4)+e(n), e(n)為一隨機(jī)噪聲，在MATLAB中可以用隨機(jī)函數(shù)rand產(chǎn)生 解 用MATLAB計算程序如下： x=[2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2]; y=[0 0 0 0 2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2]; k=length(y); e=rand(1,k)-0.5; y=y+e; 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(幅度); 計算結(jié)果如圖2.3(a)，我們看到最大值出現(xiàn)在m=4處，正好是y(n)對于x(n)的延遲。2. 3(b)是x(n)的自相關(guān)函數(shù)，他和y(n)的區(qū)別除時間位置外，形狀也略不同，這是由于y(n)受到噪聲的干擾。 實(shí)驗(yàn)三 無限長單位脈沖響應(yīng)(IIR)濾波器的設(shè)計方法 例3-1 設(shè)采樣周期T=250μs（采樣頻率fs =4kHz），用脈沖響應(yīng)不變法和雙線性變換法設(shè)計一個三階巴特沃茲濾波器，其3dB邊界頻率為fc =1kHz。 [B,A]=butter(3,2*pi*1000,s); [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); [B,A]=butter(3,2/0.00025,s); [num2,den2]=bilinear(B,A,4000); [h2,w]=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),-.,f,abs(h2),-); grid; xlabel(頻率/Hz ) ylabel(幅值/dB) 程序中第一個butter的邊界頻率2π1000，為脈沖響應(yīng)不變法原型低通濾波器的邊界頻率；第二個butter的邊界頻率2/T=2/0.00025，為雙線性變換法原型低通濾波器的邊界頻率.圖3.1給出了這兩種設(shè)計方法所得到的頻響，虛線為脈沖響應(yīng)不變法的結(jié)果；實(shí)線為雙線性變換法的結(jié)果。脈沖響應(yīng)不變法由于混疊效應(yīng)，使得過渡帶和阻帶的衰減特性變差，并且不存在傳輸零點(diǎn)。同時，也看到雙線性變換法，在z=-1即ω=π或f=2000Hz處有一個三階傳輸零點(diǎn)，這個三階零點(diǎn)正是模擬濾波器在Ω=∞處的三階傳輸零點(diǎn)通過映射形成的。 例3-2 設(shè)計一數(shù)字高通濾波器，它的通帶為400～500Hz，通帶內(nèi)容許有0.5dB的波動，阻帶內(nèi)衰減在小于317Hz的頻帶內(nèi)至少為19dB，采樣頻率為1,000Hz。 wc=2*1000*tan(2*pi*400/(2*1000)); wt=2*1000*tan(2*pi*317/(2*1000)); [N,wn]=cheb1ord(wc,wt,0.5,19,s); [B,A]=cheby1(N,0.5,wn,high,s); [num,den]=bilinear(B,A,1000); [h,w]=freqz(num,den); f=w/pi*500; plot(f,20*log10(abs(h))); axis([0,500,-80,10]); grid; xlabel() ylabel(幅度/dB) 圖3.2給出了MATLAB計算的結(jié)果，可以看到模擬濾波器在Ω=∞處的三階零點(diǎn)通過高通變換后出現(xiàn)在ω=0（z=1）處，這正是高通濾波器所希望得到的。 例3-3 設(shè)計一巴特沃茲帶通濾波器，其３dB邊界頻率分別為f2=110kHz和f1=90kHz，在阻帶f3 = 120kHz處的最小衰減大于１０dB，采樣頻率fs=400kHz。 w1=2*400*tan(2*pi*90/(2*400)); w2=2*400*tan(2*pi*110/(2*400)); wr=2*400*tan(2*pi*120/(2*400)); [N,wn]=buttord([w1 w2],[0 wr],3,10,s); [B,A]=butter(N,wn,s); [num,den]=bilinear(B,A,400); [h,w]=freqz(num,den); f=w/pi*200; plot(f,20*log10(abs(h))); axis([40,160,-30,10]); grid; xlabel(頻率/kHz) ylabel(幅度/dB) 圖3.3給出了MATLAB計算的結(jié)果，可以看出數(shù)字濾波器將無窮遠(yuǎn)點(diǎn)的二階零點(diǎn)映射為z=1的二階零點(diǎn)，數(shù)字帶通濾波器的極點(diǎn)數(shù)是模擬低通濾波器的極點(diǎn)數(shù)的兩倍。 例3-4 一數(shù)字濾波器采樣頻率fs = 1kHz，要求濾除100Hz的干擾，其３dB的邊界頻率為95Hz和105Hz，原型歸一化低通濾波器為 w1=95/500; w2=105/500; [B,A]=butter(1,[w1, w2],stop); [h,w]=freqz(B,A); f=w/pi*500; plot(f,20*log10(abs(h))); axis([50,150,-30,10]); grid; xlabel(頻率/Hz) ylabel(幅度/dB) 圖3.4為MATLAB的計算結(jié)果 　 實(shí)驗(yàn)四 有限長單位脈沖響應(yīng)(FIR)濾波器的設(shè)計方法 例2 用凱塞窗設(shè)計一FIR低通濾波器，低通邊界頻率 ，阻帶邊界頻率 ，阻帶衰減 不小于50dB。 解 首先由過渡帶寬和阻帶衰減 來決定凱塞窗的N和 圖4.1給出了以上設(shè)計的頻率特性，(a) 為N=30直接截取的頻率特性(b)為凱塞窗設(shè)計的頻率特性。凱塞窗設(shè)計對應(yīng)的MATLAB程序?yàn)椋? wn=kaiser(30,4.55); nn=[0:1:29]; alfa=(30-1)/2; hd=sin(0.4*pi*(nn-alfa))./(pi*(nn-alfa)); h=hd.*wn; [h1,w1]=freqz(h,1); plot(w1/pi,20*log10(abs(h1))); axis([0,1,-80,10]); grid; xlabel(歸一化頻率/p) ylabel(幅度/dB) 例4-2 利用雷米茲交替算法，設(shè)計一個線性相位低通FIR數(shù)字濾波器，其指標(biāo)為：通帶邊界頻率fc=800Hz，阻帶邊界fr=1000Hz，通帶波動 阻帶最小衰減At=40dB，采樣頻率fs=4000Hz。 解 在MATLAB中可以用remezord 和remez兩個函數(shù)設(shè)計，其結(jié)果如圖4.2，MATLAB程序如下： fedge=[800 1000]; mval=[1 0]; dev=[0.0559 0.01]; fs=4000; [N,fpts,mag,wt]=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); [h,w]=freqz(b,1,256); plot(w*2000/pi,20*log10(abs(h))); grid; xlabel(頻率/Hz) ylabel(幅度/dB) 函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率，數(shù)組mval是兩個邊界處的幅值，而數(shù)組dev是通帶和阻帶的波動，fs是采樣頻率單位為Hz。 第5章 數(shù)字信號處理系統(tǒng)的實(shí)現(xiàn) 例5-1求下列直接型系統(tǒng)函數(shù)的零、極點(diǎn)，并將它轉(zhuǎn)換成二階節(jié)形式 解 用MATLAB計算程序如下： num=[1 -0.1 -0.3 -0.3 -0.2]; den=[1 0.1 0.2 0.2 0.5]; [z,p,k]=tf2zp(num,den); m=abs(p); disp(零點(diǎn));disp(z); disp(極點(diǎn));disp(p); disp(增益系數(shù));disp(k); sos=zp2sos(z,p,k); disp(二階節(jié));disp(real(sos)); zplane(num,den) 輸入到“num”和“den”的分別為分子和分母多項(xiàng)式的系數(shù)。計算求得零、極點(diǎn)增益系數(shù)和二階節(jié)的系數(shù)： 零點(diǎn) 0.9615 -0.5730 -0.1443 + 0.5850i -0.1443 - 0.5850i 極點(diǎn) 0.5276 + 0.6997i 0.5276 - 0.6997i -0.5776 + 0.5635i -0.5776 - 0.5635i 增益系數(shù) 1 二階節(jié) 1.0000 -0.3885 -0.5509 1.0000 1.1552 0.6511 1.0000 0.2885 0.3630 1.0000 -1.0552 0.7679 系統(tǒng)函數(shù)的二階節(jié)形式為： 極點(diǎn)圖見圖5.1。 　 　 　 　 例5-2 分析五階橢圓低通濾波器的量化效應(yīng)，其截止頻率為0.4 ，通帶紋波為0.4dB，最小的阻帶衰減為50dB。對濾波器進(jìn)行截尾處理時，使用函數(shù)a2dT.m.。 解 用以下MATLAB程序分析量化效應(yīng) clf; [b,a]=ellip(5,0.4,50,0.4); [h,w]=freqz(b,a,512); g=20*log10(abs(h)); bq=a2dT(b,5); aq=a2dT(a,5); [hq,w]=freqz(bq,aq,512); gq=20*log10(abs(hq)); plot(w/pi,g,b,w/pi,gq,r:); grid; axis([0 1 -80 5]); xlabel(\omega/\pi); ylabel(Gain, dB); legend(量化前,量化后); figure [z1,p1,k1] = tf2zp(b,a); [z2,p2,k2] = tf2zp(bq,aq); zplaneplot([z1,z2],[p1,p2],{o,x,*,+}); legend(量化前的零點(diǎn),量化后的零點(diǎn),量化前的極點(diǎn),量化后的極點(diǎn)); 圖5.1（a）表示系數(shù)是無限精度的理想濾波器的頻率響應(yīng)（以實(shí)線表示）以及當(dāng)濾波器系數(shù)截尾到5位時的頻率響應(yīng)（以短線表示）。由圖可知，系數(shù)量化對頻帶的邊緣影響較大，經(jīng)系數(shù)量化后，增加了通帶的波紋幅度，減小了過渡帶寬，并且減小了最小的阻帶衰減。 圖5. 1（b）給出了系數(shù)量化以前和系數(shù)量化以后的橢圓低通濾波器的零極點(diǎn)位置。由圖可知，系數(shù)的量化會使零極點(diǎn)的位置與它們的理想的標(biāo)稱位置相比發(fā)生顯著的改變。在這個例子中，靠近虛軸的零點(diǎn)的位置變動最大，并且移向靠它最近的極點(diǎn)的位置。只要對程序稍作改變就可以分析舍入量化的影響。 　 為了研究二進(jìn)制數(shù)量化效應(yīng)對數(shù)字濾波器的影響，首先需要將十進(jìn)制表示的濾波器系數(shù)轉(zhuǎn)換成二進(jìn)制數(shù)并進(jìn)行量化，二進(jìn)制數(shù)的量化既可以通過截尾法也可以通過舍入法實(shí)現(xiàn)。我們提供了如下的兩個MATLAB程序：a2dT.m和a2dR.m，這兩段程序分別將向量d中的每一個數(shù)按二進(jìn)制數(shù)進(jìn)行截尾或舍入量化，量化的精度是小數(shù)點(diǎn)以后保留b位，量化后返回的向量為beq。 function beq = a2dT(d,b) % beq = a2dT(d,b) 將十進(jìn)制數(shù)利用截尾法得到b位的二進(jìn)制數(shù)， %然后將該二進(jìn)制數(shù)再轉(zhuǎn)換為十進(jìn)制數(shù) m=1; d1=abs(d); while fix(d1)>0 d1=abs(d)/(2^m); m=m+1; end beq=fix(d1*2^b); beq=sign(d).*beq.*2^(m-b-1); function beq=a2dR(d,b) % beq=a2dR(d,b)將十進(jìn)制數(shù)利用舍入法得到b位的二進(jìn)制數(shù) %然后將該二進(jìn)制數(shù)再轉(zhuǎn)換為十進(jìn)制數(shù) m=1; d1=abs(d); while fix(d1)>0 d1=abs(d)/(2^m); m=m+1; end beq=fix(d1*2^b+.5); beq=sign(d).*beq.*2^(m-b-1); 第7章 多采樣率信號處理 例7-1在時域上顯示一個 ,信號頻率為0.042 的正弦信號，然后以抽取因子3降采樣率，并在時域上顯示相應(yīng)的結(jié)果，比較兩者在時域上的特點(diǎn)。 解 用MATLAB計算程序如下： M=3; %down-sampling factor=3; fo=0.042;%signal frequency=0.042; %generate the input sinusoidal sequence n=0:N-1; m=0:N*M-1; x=sin(2*pi*fo*m); %generate the down-sampling squence y=x([1:M:length(x)]); subplot(2,1,1) stem(n,x(1:N)); title(輸入序列); xlabel(時間/n); ylabel(幅度); subplot(2,1,2) stem(n,y); title([輸出序列,抽取因子為,num2str(M)]); xlabel(時間/n); ylabel(幅度); 圖7.1 ,信號頻率為0.042 的正弦信號和降采樣率后的情況 例7-2 用漢明窗設(shè)計一長度為32的線性相位QMF濾波器組。 解 采用MATLAB設(shè)計，調(diào)用fir2函數(shù)設(shè)計公共低通濾波器，參數(shù)缺省即為漢明窗，程序如下： b1=fir2(31,[0,0.4,0.5,0.55,0.6,1],[1,1,1,0.06,0,0]); for k=1:32 b2(k)=((-1)^(k-1))*b1(k); end [H1z,w]=freqz(b1,1,256); h1=abs(H1z); g1=20*log10(h1); [H2z,w]=freqz(b2,1,256); h2=abs(H2z); g2=20*log10(h2); figure(1); plot(w/pi,g1,-,w/pi,g2,--); axis([0,1,-100,10]); grid xlabel(\omega/\pi);ylabel(幅度,dB); sum=h1.*h1+h2.*h2; d=10*log10(sum); figure(2) plot(w/pi,d);grid; xlabel(\omega/\pi);ylabel(誤差,dB); axis([0,1,-0.3,0.3]); 圖7.2(a)是一個N=32的漢明窗設(shè)計結(jié)果，圖中實(shí)線表示 = 的低通頻響，虛線表示它的鏡像 = 。圖7.2 (b)是基于這種設(shè)計方法的分析/綜合濾波器組的整個頻響 + 。從這個圖可見，重建誤差小于0.05dB。由于漢明窗設(shè) 計的頻 率響應(yīng)在通帶中近乎是平坦的，因此最大重建誤差發(fā)生在這個濾波器的通帶邊界和過渡帶內(nèi)。