數(shù)字信號處理Matlab實現(xiàn)實例(推薦給學(xué)生)

1、數(shù)字信號處理Matlab 實現(xià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給出了卷積

2、結(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)的單位脈沖響應(yīng)。 例1-3 用MATLAB計算例1-2差分方

3、程 所對應(yīng)的系統(tǒng)函數(shù)的DTFT。 解 例1-2差分方程所對應(yīng)的系統(tǒng)函數(shù)為：                                其DTFT為            

4、60;                  用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('幅度

5、')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

6、9;);ylabel('弧度')                   第章 離散傅里葉變換及其快速算法例2-1 解 此時離散序列 ，即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('

7、;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);          &#

8、160;              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é)果，由于截取了兩個整周期，得到單一譜線的頻譜。上述頻

9、譜的誤差主要是由于時域中對信號的非整周期截斷產(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).*y

10、k); 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

11、 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)，我們看到最大值出

12、現(xiàn)在m=4處，正好是y(n)對于x(n)的延遲。2. 3(b)是x(n)的自相關(guān)函數(shù)，他和y(n)的區(qū)別除時間位置外，形狀也略不同，這是由于y(n)受到噪聲的干擾。第3章 無限長單位脈沖響應(yīng)(IIR)濾波器的設(shè)計方法  例3-1 設(shè)采樣周期T=250s（采樣頻率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

13、.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，為雙線性變換法原型低通濾波器的

14、邊界頻率.圖3.1給出了這兩種設(shè)計方法所得到的頻響，虛線為脈沖響應(yīng)不變法的結(jié)果；實線為雙線性變換法的結(jié)果。脈沖響應(yīng)不變法由于混疊效應(yīng)，使得過渡帶和阻帶的衰減特性變差，并且不存在傳輸零點。同時，也看到雙線性變換法，在z=-1即=或f=2000Hz處有一個三階傳輸零點，這個三階零點正是模擬濾波器在=處的三階傳輸零點通過映射形成的。 例3-2 設(shè)計一數(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*t

15、an(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é)果，可以看到模擬濾波器在=處的三階零

16、點通過高通變換后出現(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

17、*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é)果，可

18、以看出數(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

19、9;) ylabel('幅度/dB')            圖3.4為MATLAB的計算結(jié)果                  第4章 有限長單位脈沖響應(yīng)(FIR)濾波器的設(shè)計方法 例2 用凱塞窗設(shè)計一FIR低通濾波器，低通邊界頻率 ，阻帶邊界頻率 ，阻帶衰減 不小于50dB。 解 首先由過渡帶寬和阻帶衰減 來決定凱塞窗

20、的N和         圖4.1給出了以上設(shè)計的頻率特性，(a) 為N=30直接截取的頻率特性(b)為凱塞窗設(shè)計的頻率特性。凱塞窗設(shè)計對應(yīng)的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

21、') ylabel('幅度/dB') 例4-2 利用雷米茲交替算法，設(shè)計一個線性相位低通FIR數(shù)字濾波器，其指標為：通帶邊界頻率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;           

22、60;     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;                        

23、               xlabel('頻率/Hz') ylabel('幅度/dB')                        函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率，數(shù)組mv

24、al是兩個邊界處的幅值，而數(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é)&#

25、39;);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

26、.3630 1.0000 -1.0552 0.7679    系統(tǒng)函數(shù)的二階節(jié)形式為： 極點圖見圖5.1。      例5-2 分析五階橢圓低通濾波器的量化效應(yīng)，其截止頻率為0.4 ，通帶紋波為0.4dB，最小的阻帶衰減為50dB。對濾波器進行截尾處理時，使用函數(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

27、);         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,

28、p2,'o','x','*','+');legend('量化前的零點','量化后的零點','量化前的極點','量化后的極點');       圖5.1（a）表示系數(shù)是無限精度的理想濾波器的頻率響應(yīng)（以實線表示）以及當濾波器系數(shù)截尾到5位時的頻率響應(yīng)（以短線表示）。由圖可知，系數(shù)量化對頻帶的邊緣影響較大，經(jīng)系數(shù)量化后，增加了通帶的波紋幅度，減小了過渡帶寬，并且減小了最小的阻帶衰減。   &#

29、160; 圖5. 1（b）給出了系數(shù)量化以前和系數(shù)量化以后的橢圓低通濾波器的零極點位置。由圖可知，系數(shù)的量化會使零極點的位置與它們的理想的標稱位置相比發(fā)生顯著的改變。在這個例子中，靠近虛軸的零點的位置變動最大，并且移向靠它最近的極點的位置。只要對程序稍作改變就可以分析舍入量化的影響。              為了研究二進制數(shù)量化效應(yīng)對數(shù)字濾波器的影響，首先需要將十進制表示的濾波器系數(shù)轉(zhuǎn)換成二進制數(shù)并進行量化，二進制數(shù)的量化既可以通過截尾法也可以通過舍入法實現(xiàn)。我們提供了如下

30、的兩個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(

31、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降采樣率，并在時域上顯示相應(yīng)的結(jié)果，比較兩者在時域上的特點。 解 用MATLAB計算程序如下：M=3; %down-sampling factor=3;fo=0.042;%sign

32、al 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