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

實驗一 離散時間信號與系統(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計算差分方程 當輸入序列為 時的輸出結(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個樣點的輸出，即該系統(tǒng)的單位脈沖響應。 例1-3 用MATLAB計算例1-2差分方程 所對應的系統(tǒng)函數(shù)的DTFT。 解 例1-2差分方程所對應的系統(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);gridtitle(實部)xlabel(omega/pi);ylabel(幅度)subplot(2,2,2);plot(w/pi,imag(h);gridtitle(虛部)xlabel(omega/pi);ylabel(Amplitude)subplot(2,2,3);plot(w/pi,abs(h);gridtitle(幅度譜)xlabel(omega/pi);ylabel(幅值)subplot(2,2,4);plot(w/pi,angle(h);gridtitle(相位譜)xlabel(omega/pi);ylabel(弧度) 實驗二 離散傅里葉變換及其快速算法 例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)為一隨機噪聲，在MATLAB中可以用隨機函數(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)受到噪聲的干擾。實驗三 無限長單位脈沖響應(IIR)濾波器的設計方法 例3-1 設采樣周期T=250s（采樣頻率fs =4kHz），用脈沖響應不變法和雙線性變換法設計一個三階巴特沃茲濾波器，其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的邊界頻率21000，為脈沖響應不變法原型低通濾波器的邊界頻率；第二個butter的邊界頻率2/T=2/0.00025，為雙線性變換法原型低通濾波器的邊界頻率.圖3.1給出了這兩種設計方法所得到的頻響，虛線為脈沖響應不變法的結(jié)果；實線為雙線性變換法的結(jié)果。脈沖響應不變法由于混疊效應，使得過渡帶和阻帶的衰減特性變差，并且不存在傳輸零點。同時，也看到雙線性變換法，在z=-1即=或f=2000Hz處有一個三階傳輸零點，這個三階零點正是模擬濾波器在=處的三階傳輸零點通過映射形成的。 例3-2 設計一數(shù)字高通濾波器，它的通帶為400500Hz，通帶內(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é)果，可以看到模擬濾波器在=處的三階零點通過高通變換后出現(xiàn)在=0（z=1）處，這正是高通濾波器所希望得到的。 例3-3 設計一巴特沃茲帶通濾波器，其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ù)字濾波器將無窮遠點的二階零點映射為z=1的二階零點，數(shù)字帶通濾波器的極點數(shù)是模擬低通濾波器的極點數(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é)果 實驗四 有限長單位脈沖響應(FIR)濾波器的設計方法 例2 用凱塞窗設計一FIR低通濾波器，低通邊界頻率 ，阻帶邊界頻率 ，阻帶衰減 不小于50dB。 解 首先由過渡帶寬和阻帶衰減 來決定凱塞窗的N和 圖4.1給出了以上設計的頻率特性，(a) 為N=30直接截取的頻率特性(b)為凱塞窗設計的頻率特性。凱塞窗設計對應的MATLAB程序為： 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 利用雷米茲交替算法，設計一個線性相位低通FIR數(shù)字濾波器，其指標為：通帶邊界頻率fc=800Hz，阻帶邊界fr=1000Hz，通帶波動 阻帶最小衰減At=40dB，采樣頻率fs=4000Hz。 解 在MATLAB中可以用remezord 和remez兩個函數(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)的實現(xiàn) 例5-1求下列直接型系統(tǒng)函數(shù)的零、極點，并將它轉(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(零點);disp(z);disp(極點);disp(p);disp(增益系數(shù));disp(k);sos=zp2sos(z,p,k);disp(二階節(jié));disp(real(sos);zplane(num,den)輸入到“num”和“den”的分別為分子和分母多項式的系數(shù)。計算求得零、極點增益系數(shù)和二階節(jié)的系數(shù)：零點 0.9615 -0.5730 -0.1443 + 0.5850i -0.1443 - 0.5850i 極點 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é)形式為： 極點圖見圖5.1。 例5-2 分析五階橢圓低通濾波器的量化效應，其截止頻率為0.4 ，通帶紋波為0.4dB，最小的阻帶衰減為50dB。對濾波器進行截尾處理時，使用函數(shù)a2dT.m.。解 用以下MATLAB程序分析量化效應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(量化前,量化后);figurez1,p1,k1 = tf2zp(b,a);z2,p2,k2 = tf2zp(bq,aq);zplaneplot(z1,z2,p1,p2,o,x,*,+);legend(量化前的零點,量化后的零點,量化前的極點,量化后的極點); 圖5.1（a）表示系數(shù)是無限精度的理想濾波器的頻率響應（以實線表示）以及當濾波器系數(shù)截尾到5位時的頻率響應（以短線表示）。由圖可知，系數(shù)量化對頻帶的邊緣影響較大，經(jīng)系數(shù)量化后，增加了通帶的波紋幅度，減小了過渡帶寬，并且減小了最小的阻帶衰減。 圖5. 1（b）給出了系數(shù)量化以前和系數(shù)量化以后的橢圓低通濾波器的零極點位置。由圖可知，系數(shù)的量化會使零極點的位置與它們的理想的標稱位置相比發(fā)生顯著的改變。在這個例子中，靠近虛軸的零點的位置變動最大，并且移向靠它最近的極點的位置。只要對程序稍作改變就可以分析舍入量化的影響。 為了研究二進制數(shù)量化效應對數(shù)字濾波器的影響，首先需要將十進制表示的濾波器系數(shù)轉(zhuǎn)換成二進制數(shù)并進行量化，二進制數(shù)的量化既可以通過截尾法也可以通過舍入法實現(xiàn)。我們提供了如下的兩個MATLAB程序：a2dT.m和a2dR.m，這兩段程序分別將向量d中的每一個數(shù)按二進制數(shù)進行截尾或舍入量化，量化的精度是小數(shù)點以后保留b位，量化后返回的向量為beq。 function beq = a2dT(d,b)% beq = a2dT(d,b) 將十進制數(shù)利用截尾法得到b位的二進制數(shù)，%然后將該二進制數(shù)再轉(zhuǎn)換為十進制數(shù)m=1; d1=abs(d);while fix(d1)>0 d1=abs(d)/(2m); m=m+1;endbeq=fix(d1*2b);beq=sign(d).*beq.*2(m-b-1); function beq=a2dR(d,b)% beq=a2dR(d,b)將十進制數(shù)利用舍入法得到b位的二進制數(shù)%然后將該二進制數(shù)再轉(zhuǎn)換為十進制數(shù)m=1; d1=abs(d);while fix(d1)>0 d1=abs(d)/(2m); m=m+1;endbeq=fix(d1*2b+.5);beq=sign(d).*beq.*2(m-b-1);第7章 多采樣率信號處理 例7-1在時域上顯示一個 ,信號頻率為0.042 的正弦信號，然后以抽取因子3降采樣率，并在時域上顯示相應的結(jié)果，比較兩者在時域上的特點。 解 用MATLAB計算程序如下：M=3; %down-sampling factor=3;fo=0.042;%signal frequency=0.042;%generate the input sinusoidal sequencen=0:N-1;m=0:N*M-1;x=sin(2*pi*fo*m);%generate the down-sampling squencey=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 用漢明窗設計一長度為32的線性相位QMF濾波器組。 解 采用MATLAB設計，調(diào)用fir2函數(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的漢明窗設計結(jié)果，圖中實線表示 = 的低通頻響，虛線表示它的鏡像 = 。圖7.2 (b)是基于這種設計方法的分析/綜合濾波器組的整個頻響 + 。從這個圖可見，重建誤差小于0.05dB。由于漢明窗設 計的頻率響應在通帶中近乎是平坦的，因此最大重建誤差發(fā)生在這個濾波器的通帶邊界和過渡帶內(nèi)。