數(shù)字信號(hào)處理各章節(jié)重點(diǎn)知識(shí)的matlab實(shí)例

1、1數(shù)字信號(hào)處理各章節(jié)重點(diǎn)知識(shí)的matlab實(shí)例第1, 2章離散時(shí)間信號(hào)與系統(tǒng)1的離散電W例1-1用MATLAB計(jì)算序列-2 0 1 1 3和序列1 2 0 - 卷積。解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('幅度')；圖L1用MATLAB計(jì)輯兩個(gè)序列的卷積圖1.1給出了卷積結(jié)果的圖形，求得的結(jié)果存放在數(shù)組c中為-2 -4 1 3 1 5 1 -3。y(輯)+ 0,7y5-1)一 0.45y("

2、 - 2) - 06MH - 3=。以-044m-1) + 0.36x(弁-2) + 0.024-3)當(dāng)輸入序列為市)二訕)時(shí)的輸出結(jié)果陽(yáng)，0<«<40 0解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系統(tǒng)的單位脈沖響應(yīng)。給出了該差分方程的前41個(gè)樣點(diǎn)的輸出，即該圖1.2用MATLAB計(jì)理差分方程輸出例1-3用MATLAB計(jì)

3、算例1-2差分方程了（輯）+ 07丁（"一1）一0.45歹（用一 2） - 0.6y（n - 3）=0以-044秘-1） + 0.36x（” 2） + Q02-3）所對(duì)應(yīng)的系統(tǒng)函數(shù)的DTFT解 例1-2差分方程所對(duì)應(yīng)的系統(tǒng)函數(shù)為：0.8-Q4421+Q,36J+QQ2/億尸不而匚藤石J?"其DTFT為0.8-044,用 +0.36 戶(hù)'+0.02,刖01 + 0 7”、0 4 首初-06”.用MATLA時(shí)算的程序如下：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,d

4、en, w)；subplot(2,2,1);plot(w/pi,real(h );gridtitle(' 實(shí)部')pi');ylabel(' 幅度')subplot(2,2,2);plot(w/pi,imag(h);gridtitle(' 虛部 ')pi');ylabel('Amplitude')subplot(2,2,3);plot(w/pi,abs(h);gridtitle(' 幅度譜 ')pi');ylabel(' 幅值 ')subplot(2,2,4);plot(w/

5、pi,angle( h);grid弧度 ')title(' 相位譜 ')xlabel('omega/pi');ylabel(第3、4章離散傅里葉變換及其快速算法例2-1對(duì)連續(xù)的單一頻率周期信號(hào)按采樣頻率，=”采樣，截取長(zhǎng)度N分別選N =20和N=16 ,觀(guān)察其DFT結(jié)果的幅度譜。解此時(shí)離散序列 咖二版（2時(shí)J力二疝（浙/ 8）即k=8。用MATLA時(shí)算并作圖，函數(shù)fft用于計(jì)算離散傅里葉變換DF【程序如下：k=8;n1=0:1:19;x(n)');k1);k)');x(n)');k2);k)');xa1=sin(2*pi*

6、n1/k);subplot(2,2,1)plot(n1,xa1)xlabel('t/T');ylabel('xk1=fft(xa1);xk1=abs(xsubplot(2,2,2)stem(n1,xk1)xlabel('k');ylabel('X(n2=0:1:15;xa2=sin(2*pi*n2/k);subplot(2,2,3)plot(n2,xa2)xlabel('t/T');ylabel('xk2=fft(xa2);xk2=abs(xsubplot(2,2,4)stem(n2,xk2)xlabel('k&#

7、39;);ylabel('X(©解用MATLAB計(jì)算程序如下:圖2.1不同獻(xiàn)取長(zhǎng)度的正弦信號(hào)及其DFT結(jié)果計(jì)算結(jié)果示于圖2.1, (a)和(b)分別是N=20時(shí)的截取信號(hào)和DFT結(jié)果，由 于截取了兩個(gè)半周期，頻譜出現(xiàn)泄漏；(c)和(d)分別是N=16時(shí)的截取信號(hào)和 DFT結(jié)果，由于截取了兩個(gè)整周期，得到單一譜線(xiàn)的頻譜。上述頻譜的誤差主要 是由于時(shí)域中對(duì)信號(hào)的非整周期截?cái)喈a(chǎn)生的頻譜泄漏。二例2-2用FFT計(jì)算兩個(gè)序列則二2 1-1120-1的互相關(guān)函數(shù)刈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

8、);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(' 幅度 ');其計(jì)算結(jié)果如圖 2.2 所示。121086幅4度20- 2- 4- 6-B -6-4-202438圖21兩個(gè)序列的自相關(guān)函數(shù)W例2-3計(jì)算兩個(gè)序列的的互相關(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ī)噪聲，在MATLAB 中可以用隨

9、機(jī)函數(shù)rand解用MATLAB計(jì)算程序如下：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('幅度');計(jì)算結(jié)果如圖

10、2.3(a),我們看到最大值出現(xiàn)在 m=4處，正好是y(n)對(duì)于x(n)y(n)的區(qū)別除時(shí)間位置外，形狀也的延遲。2. 3(b)是x(n)的自相關(guān)函數(shù)，他和 略不同，這是由于y(n)受到噪聲的干擾。幅度8)圖2. 3延遲序列的互相關(guān)函數(shù)色)和自相關(guān)函數(shù)(b)第5章 無(wú)限長(zhǎng)單位脈沖響應(yīng)(IIR)濾波器的設(shè)計(jì)方法2例3-1設(shè)采樣周期T=250巧(采樣頻率fs =4kHz),用脈沖響應(yīng)不變法 和雙線(xiàn)性變換法設(shè)計(jì)一個(gè)三階巴特沃茲濾波器，其 3dB邊界頻率為fc =1kHz。B,A=butter(3,2*pi*1000,'s');num1,den1=impinvar(B,A,4000);

11、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')程序中第一個(gè)butter 的邊界頻率2冗X1000 ,為脈沖響應(yīng)不變法原型低通濾波器的邊界頻率；第二個(gè)butter的邊界頻率2/T=2/0.00025

12、,為雙線(xiàn)性變換法原型低通濾波器的邊界頻率.圖3.1給出了這兩種設(shè)計(jì)方法所得到的頻響，虛線(xiàn)為脈沖響應(yīng)不變法的結(jié)果；實(shí)線(xiàn)為雙線(xiàn)性變換法的結(jié)果。脈沖響應(yīng)不變法由 于混疊效應(yīng)，使得過(guò)渡帶和阻帶的衰減特性變差，并且不存在傳輸零點(diǎn)。同時(shí)， 也看到雙線(xiàn)性變換法，在z=-1即=冗或f=2000Hz處有一個(gè)三階傳輸零點(diǎn)，這 個(gè)三階零點(diǎn)正是模擬濾波器在Q =8處的三階傳輸零點(diǎn)通過(guò)映射形成的。1.1 - 111111»111HliAIIHHN0.94k -1 一4一一»一1一！H19HIIn a11：11H|1111nhU.o11：*1IHIiij|i107.-m-一十三一一一| |) I1iP

13、Ii»N <1111用0.6,一|iH：;sH|1|iH理 0.5中IlH： ： ：H11 . 10 41 - 1s_ »N1 k !|1| ： > ：；0.3一串T-r - »A- = T - =7|i-11X1，0.2三三一ii11-i0.1ii1- - -I-i-11 1 .:xj :0Iin! r'0200400500 800 1000 1200 1400 1SOO 1800 2000頻率ZHz圖3.1三階巴特沃茲濾波器的頻率響應(yīng)逐例3-2設(shè)計(jì)一數(shù)字高通濾波器，它的通帶為 400500Hz,通帶內(nèi)容許 有0.5dB的波動(dòng)，阻帶內(nèi)衰減在小

14、于 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(&

15、#39;')ylabel(' 幅度/dB')圖3.2給出了 MATLAB計(jì)算的結(jié)果，可以看到模擬濾波器在Q =8處的三階零點(diǎn) 通過(guò)高通變換后出現(xiàn)在 必=0 (z=1)處，這正是高通濾波器所希望得到的。圖3.2切比雪夫高通濾波器底例3-3設(shè)計(jì)一巴特沃茲帶通濾波器,其3 dB邊界頻率分別為f2=110kHz和f1=90kHz,在阻帶f3 = 120kHz處的最小衰減大于1 0 dB,采樣頻率fs=400kHz。10圖3.3巴特沃茲帶通痣波器w1=2*400*tan(2*pi*90/(2*400);w2=2*400*tan(2*pi*110/(2*400);wr=2*400*

16、tan(2*pi*120/(2*400);'s');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(40,160,-30,10);grid;xlabel(' 頻率/kHz')ylabel(' 幅度/dB')圖3.3給出了 MATLAB計(jì)算的結(jié)果，可以看出數(shù)字濾波器將無(wú)窮遠(yuǎn)點(diǎn)的二階零 點(diǎn)映射為z=+的二階零點(diǎn)，數(shù)

17、字帶通濾波器的極點(diǎn)數(shù)是模擬低通濾波器的極點(diǎn) 數(shù)的兩倍。W例3-4 數(shù)字濾波器采樣頻率fs= 1kHz,要求濾除100Hz的干擾，其3 dB的邊界頻率為95Hz和105Hz,原型歸一化低通濾波器為1 + 5w1=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的計(jì)算結(jié)果10加

18、LIJJLLKIJL3090703090100110120113014)130頻率ZH,圖3.4巴特沃茲帶阻濾波器第5章 有限長(zhǎng)單位脈沖響應(yīng)（FIR）濾波器的設(shè)計(jì)方法孤例4.1用凱塞窗設(shè)計(jì)一 FIR低通濾波器，低通邊界頻率二0一%，阻 帶邊界頻率 嗎;。，加，阻帶衰減從 不小于50dB。解 首先由過(guò)渡帶寬和阻帶衰減 丸 來(lái)決定凱塞窗的N和£= -iu = 0 2開(kāi)r1N= 50-8 附 302 285x0 2江 = 0.1102(50-8.7) = 4.55UN fl A 口鼻 H9 OB Ojr J J 歸一北期品，真不號(hào)勝磔歸一優(yōu)西塞/宣圖4 1凱塞窗設(shè)計(jì)舉例圖4.1給出了以上設(shè)

19、計(jì)的頻率特性，(a)為N=30直接截取的頻率特性(b)為 凱塞窗設(shè)計(jì)的頻率特性。凱塞窗設(shè)計(jì)對(duì)應(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,10);grid;xlabel('歸一化頻率/兀')ylabel(' 幅度/dB')2例4-2利用雷米茲交替算法，設(shè)計(jì)一個(gè)線(xiàn)性相位低通FIR數(shù)字濾波

20、器， 其指標(biāo)為：通帶邊界頻率fc=800Hz;阻帶邊界f=1000Hz;通帶波動(dòng)3二05dB阻 帶最小衰減 At=40dB,采樣頻率fs=4000Hz- 一一解.,在MATLA呻可以用remezord和remez兩個(gè)函數(shù)設(shè)計(jì)，其結(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);g

21、rid;xlabel('頻率/Hz')ylabel('幅度/dB')ID圖42Remez交替法設(shè)計(jì)舉例函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界頻率，數(shù)組 mval是兩個(gè)邊界處的幅值，而數(shù)組 dev是通帶和阻帶的波動(dòng)，fs是采樣頻率單 位為Hz。6章 數(shù)字信號(hào)處理系統(tǒng)的實(shí)現(xiàn)咆例5-1求下列直接型系統(tǒng)函數(shù)的零、極點(diǎn)，并將它轉(zhuǎn)換成二階節(jié)形式、l-O.lz-OJz-O.-O 4H =;571 + O.lz4 +0.2Z-3 +0.2z +0.5z*解用MATLA臥算程序如下：num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0

22、.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ù)。計(jì)算求得零、 極點(diǎn)增益系數(shù)和二階節(jié)的系數(shù)：零點(diǎn)0.9615-0.5730- 0.1443 + 0.5850i- 0.1443 - 0.5850i極點(diǎn)0

23、.5276 + 0.6997i0.5276 - 0.6997i- 0.5776 + 0.5635i- 0.5776 - 0.5635i增益系數(shù)二階節(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é)形式為:月1一 03 鴕 5-0.5509/1 + 1.1552/+ 0.651121 + 0 2885/+0 3630/1-1 0552/+076792“極點(diǎn)圖見(jiàn)圖5.1。圖5.1系統(tǒng)函數(shù)的零、極點(diǎn)圖2例5-2分析五階橢圓低通濾波器的量化效應(yīng)，其截止頻率為0.4

24、乃,通帶紋波為0.4dB,最小的阻帶衰減為50dB。對(duì)濾波器進(jìn)行截尾處理時(shí)，使用函 數(shù) a2dT.m.。解用以下MATLABi序分析量化效應(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');yla

25、bel('Gain, dB');legend('量化前,'量化后')；figurez1,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ù)是無(wú)限精度的理想濾波器的頻率響應(yīng)(以實(shí)線(xiàn)表示)以及當(dāng)濾波器系數(shù)截尾到5位時(shí)的頻率響應(yīng)(以短線(xiàn)表示)。由圖可知，系數(shù)量化 對(duì)頻

26、帶的邊緣影響較大，經(jīng)系數(shù)量化后，增加了通帶的波紋幅度，減小了過(guò)渡帶 寬，并且減小了最小的阻帶衰減。圖5. 1 (b)給出了系數(shù)量化以前和系數(shù)量化以后的橢圓低通濾波器的零極點(diǎn)位置。由圖可知，系數(shù)的量化會(huì)使零極點(diǎn)的位置與它們的理想的標(biāo)稱(chēng)位置相比 發(fā)生顯著的改變。在這個(gè)例子中，靠近虛軸的零點(diǎn)的位置變動(dòng)最大， 并且移向靠 它最近的極點(diǎn)的位置。只要對(duì)程序稍作改變就可以分析舍入量化的影響。系數(shù)量化前后的頻率響應(yīng)14$00.51實(shí)部(b)系數(shù)量化前后的零極點(diǎn)分布 P量化前的零點(diǎn)量化后的零點(diǎn)， R量化前的極點(diǎn)量化后的極點(diǎn)圖5五階橢圓低通濾波器的量化效應(yīng)為了研究二進(jìn)制數(shù)量化效應(yīng)對(duì)數(shù)字濾波器的影響，首先需要將十進(jìn)

27、制表示的濾波器系數(shù)轉(zhuǎn)換成二進(jìn)制數(shù)并進(jìn)行量化，二進(jìn)制數(shù)的量化既可以通過(guò)截尾法也可 以通過(guò)舍入法實(shí)現(xiàn)。我們提供了如下的兩個(gè) MATLAB程序：a2dT.m和a2dR.m, 這兩段程序分別將向量d中的每一個(gè)數(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)>0d1=abs(d)/(2Am);m=m+1;endbeq=fix(d1*2Ab);b

28、eq=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章多采樣率信號(hào)處理9例7-1在時(shí)域上顯示一個(gè)N二5。，信號(hào)頻率為0.042星的正弦信 號(hào)，然后以抽取因子3降采樣率，并在時(shí)域上顯示相應(yīng)的結(jié)果，比較兩者在時(shí)域 上的特點(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);subplot(2,1,1)stem(n,x(1:N);title('輸入序列)；xlabel(' 時(shí)間/n');ylabel('幅度'