數(shù)字信號處理Matlab實(shí)現(xiàn)實(shí)例

1、實(shí)用標(biāo)準(zhǔn)文案數(shù)字信號處理Matlab實(shí)現(xiàn)實(shí)例第1章離散時間信號與系統(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('幅度')；精彩文檔n圖L1用MATLAB計其兩個序列的卷積圖1.1給出了卷積結(jié)果的圖形，求得的結(jié)果存放在數(shù)組c中為：-2 -4 1-3。例1-2用MATLAB十算差分方程了(售)+ 0.7,(差 T) _

2、0,45/(并2) - 0 6y(犀一 35=0&種- 0,44x(h-1) + 036KM -2) + 0 02x(« - 3)當(dāng)輸入序列為 則二帥)時的輸出結(jié)果y(n 0<«<40 。解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)

3、。圖1.2用MATLAB計菖差分方程輸出例1-3用MATLAB十算例1-2差分方程MM + O'MkT)戲-2)_ 0,6/(符=0.取- 044x(-1) + 036x(-2) + 0.02x(«- 3)所對應(yīng)的系統(tǒng)函數(shù)的 DTFT解 例1-2差分方程所對應(yīng)的系統(tǒng)函數(shù)為：_ _j _q_3H(z)=1 0.7z,-0.45z，-0.6z工0.8 -0.44z0.36z0.02z其DTFT為j 0.8 -0.44ej0.36e-20.02e-j3 e 1 0.7e -0.45e-j2 -0.6e-j3 ,用MATLAB十算的程序如下: k=256;num=0.8 -0.44

4、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('實(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

5、(h);gridtitle(' 幅度譜)xlabel('omega/pi');ylabel('幅值')subplot(2,2,4);plot(w/pi,angle(h);gridtitle(' 相位譜)xlabel('omegaApi');ylabel('弧度')第2章離散傅里葉變換及其快速算法例2-1對連續(xù)的單一頻率周期信號按采樣頻率 1=8/ 采樣，截取長度N分別選N=20和N=16 ,觀察其DFT結(jié)果的幅度譜。解此時離散序列= sin( 2喇/工)=渺(2就/8) ,即k=8。用MATLAB十算并作圖，函數(shù)f

6、ft用于計算離散傅里葉變換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');yl

7、abel('x(n)');xk2=fft(xa2);xk2=abs(xk2);subplot(2,2,4)stem(n2,xk2)xlabel('k');ylabel('X(k)');106642005101510|8o6 W-*受42051015圖2.1不同威取長度的正弦信號及其DFT結(jié)果計算結(jié)果示于圖2.1 , (a)和(b)分別是N=20時的截取信號和 DFT結(jié)果,由于截取了兩個半周期,頻譜出現(xiàn)泄漏；(c)和(d)分別是N=16時的截取信號和 DFT結(jié)果,由于截取了兩個整周期，得到單一譜線的頻譜。非整周期截斷產(chǎn)生的頻譜泄漏。電工例2-2用F

8、FT計算兩個序列上述頻譜的誤差主要是由于時域中對信號的x(n) = 1 3-112帥):(2 1-1的互相關(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所示。圖22兩個序列的自相關(guān)函數(shù)例

9、2-3計算兩個序列的的互相關(guān)函數(shù)，其中x(n尸2 3 5 2 1- 1 0 0 12 3 5 3 0- 1 - 2 0 1 2y(n)=x(n-4)+e(n), e(n)為一隨機(jī)噪聲，在MATLAB3可以用隨機(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(

10、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)受到噪聲的干擾。圖2. 3延遲序列的互相關(guān)函數(shù)值)和自相關(guān)函數(shù)(b)第3章 無限長單位脈沖響應(yīng)(IIR)濾波器的設(shè)計方法例3-1設(shè)采樣周期T=250s (采樣頻率fs=4kHz),用脈沖響應(yīng)不變法和雙線性變換 法設(shè)計一個三

11、階巴特沃茲濾波器，其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 ')ylabe

12、l('幅值 /dB')程序中第一個butter的邊界頻率2兀x 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即3=?；騠=2000Hz處有一個 三階傳輸零點(diǎn)，這個三階零點(diǎn)正是模擬濾波器在Q=8處的三階傳輸零點(diǎn)通過映射形成的。1, 1 - 1-1Illi|i

13、i1ddIIJHH09;:i -1 一一一 一J=-T一 一ii1ghiin a11i1nIi1Iiiiiiii11111KHlidi07.-三胃一 T三三十三三 | |11HliNM1' 18N限06 HH1 1!;II|1HH111%iNIi理 05HH1b11 1,1110.41 . -1 _ iih11 1_i_!|l|i| :、；0.3-|irrA- - - T - - -X-1Ii11 1、j HN1i 11J n0.2Ih1 11 X 11'i aiiii0.1l!l!-一-1-IlII11111ii-.iiri- -ii0 -11IIl>0|l|l1i|(

14、1|l1rih020040060080口1000 1200 1400 1600 1300 2000頻率圖3 1三階巴特沃茲濾波器的頻率響應(yīng) 動，阻帶內(nèi)衰減在小于 317Hz的頻帶內(nèi)至少為19dB,采樣頻率為1,000Hz。例3-2 設(shè)計一數(shù)字高通濾波器，它的通帶為400500Hz,通帶內(nèi)容許有 0.5dB的波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' , 

15、9;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é)果，可以看到模擬濾波器在=8處的三階零點(diǎn)通過高通變換后出現(xiàn)在3 =0 (z=1)處，這正是高通濾波器所希望得到的。頻率/取圖12切比雪夫高通濾波器例3-3設(shè)計一巴特沃茲帶通濾波器，其3dB邊界頻率分別為f2=110kHz和f i=90kHz,在阻帶f3 =

16、120kHz處的最小衰減大于1 0dB,采樣頻率fs=400kHz。10頻率/kHz§ -s媽 里io圖3.3巴特沃茲帶通法波器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,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

17、(40,160,-30,10);'s')；grid；xlabel('ylabel('圖3.3給出了1的二階零點(diǎn),頻率/kHz')幅度/dB')MATLAB十算的結(jié)果，可以看出數(shù)字濾波器將無窮遠(yuǎn)點(diǎn)的二階零點(diǎn)映射為 數(shù)字帶通濾波器的極點(diǎn)數(shù)是模擬低通濾波器的極點(diǎn)數(shù)的兩倍。z= ±總工例3-4 一數(shù)字濾波器采樣頻率fs = 1kHz,要求濾除100Hz的干擾，其3 dB的邊界頻率為95Hz和105Hz,原型歸一化低通濾波器為母二41 + Sw1=95巧00;w2=105/500;B,A=butter(1,w1, w2, 'stop

18、9;);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é)果1050i 1i i ii i i ii ii 11 II 1 IIi i i p i h11 iiI ii i 1 1 1 1i i 1 i I iiiii一 i.piiii圖ii i1- i i -r 入4一|i i1 iH liHII1H旭 轡-10i li iI-1 |1fl ii rII-15i 丁 1

19、.1 1引-r111iI iH 111201-* -1寸11-1- ii-H-1|i一 X .Ii11 I三4=11i11 1i:11I ：1 1 ki-I - i iiI 一一i1 iIi一 X =, IJ hSO 磕?0 #90 iw 110 1X M） 130頻率/取圖3.4巴特沃茲帶阻濾波器第4章 有限長單位脈沖響應(yīng)（FIR）濾波器的設(shè)計方法例2用凱塞窗設(shè)計一FIR低通濾波器，低通邊界頻率值 =Q一%,阻帶邊界頻率叫二。.玩,阻帶衰減人不小于50dB。解 首先由過渡帶寬和阻帶衰減人 來決定凱塞窗的N和6Q) = Qy-Q)= 0 2打r150-82 2g5x0 2戀解3010 = 0,

20、1102(50-8.7) = 4.55T I I:-:n不號勝磔歸-4t戰(zhàn)賽,怎UN fl 5 口鼻 H9 OB Ojr J J歸一優(yōu)期嵬，真圖4 1凱塞窗設(shè)計舉例圖4.1給出了以上設(shè)計的頻率特性，(a)為N=30直接截取的頻率特性(b)為凱塞窗設(shè)計的頻率特性。凱塞窗設(shè)計對應(yīng)的MATLA翼序?yàn)椋簑n=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

21、,10);grid；xlabel(' 歸一化頻率/元)ylabel(' 幅度 /dB')例4-2利用雷米茲交替算法，設(shè)計一個線性相位低通FIR數(shù)字濾波器，其指標(biāo)為：通帶邊界頻率fc=800Hz,阻帶邊界fr=1000Hz,通帶波動3二05dB 阻帶最小衰減 At=40dB, 采樣頻率fs=4000Hz。解 = 1-10 = 0.0559在MATLA沖可以用remezord和remez兩個函數(shù)設(shè)計，其結(jié)果如圖4.2 , MATLA翼序如下:fedge=800 1000;mval=1 0;dev=0.0559 0.01;fs=4000;N,fpts,mag,wt=remez

22、ord(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')1002030-40如,即70ISOO 20000200400600800100012001400 KOO頻率生圖42 Remez交替法設(shè)計舉例函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率，數(shù)組 mval是兩個 邊界處的幅值，而數(shù)組dev是通帶和阻帶的波動，fs是采樣頻率單位為 Hz。第5

23、章數(shù)字信號處理系統(tǒng)的實(shí)現(xiàn)例5-1求下列直接型系統(tǒng)函數(shù)的零、極點(diǎn)，并將它轉(zhuǎn)換成二階節(jié)形式H =；-7TT1 + 01，+0 2/+0 叱 +0 5/解用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);zpl

24、ane(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.6997i0.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.65111.0000 0.2885 0.3630 1.0000 -1.0552 0.7679系統(tǒng)函數(shù)的二階節(jié)形式為:1-0.3885r】-0,5

25、509只1 + 0一2就一】+01%30 產(chǎn)1 + 1,1552z-1 +0.65Hz-2 1-1 0552rl+0.7679z-3極點(diǎn)圖見圖5.1 。0.6圖5.1系統(tǒng)函數(shù)的零、極點(diǎn)圖0.4x，通帶紋波為例5-2分析五階橢圓低通濾波器的量化效應(yīng)，其截止頻率為0.4dB,最小的阻帶衰減為50dB。對濾波器進(jìn)行截尾處理時，使用函數(shù)a2dT.m.。解用以下MATLA翼序分析量化效應(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,51

26、2);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','*','+');l

27、egend('量化前的零點(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)的位置變動最大, 程

28、序稍作改變就可以分析舍入量化的影響。并且移向靠它最近的極點(diǎn)的位置。只要對為了研究二進(jìn)制數(shù)量化效應(yīng)對數(shù)字濾波器的影響，首先需要將十進(jìn)制表示的濾波器系數(shù)轉(zhuǎn)換成二進(jìn)制數(shù)并進(jìn)行量化，二進(jìn)制數(shù)的量化既可以通過截尾法也可以通過舍入法實(shí)現(xiàn)。我們提供了如下的兩個MATLA叫序：a2dT.m和a2dR.m,這兩段程序分別將向量d中的每一個數(shù)按二進(jìn)制數(shù)進(jìn)行截尾或舍入量化, 返回的向量為beq。量化的精度是小數(shù)點(diǎn)以后保留 b位，量化后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);w

29、hile fix(d1)>0d1=abs(d)/(2Am);m=m+1;endbeq=fix(d1*2Ab)；beq=sign(d).*beq.*2A(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)>0d1=abs(d)/(2Am)；m=m+1;endbeq=fix(d1*2Ab+.5)；beq=sign(d).*beq.*2A(m-b-1)；第7章多采樣率信號處理電工例7-1在時域上顯示一個 N二50，信號頻率為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 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);subpl