數(shù)字信號(hào)處理實(shí)驗(yàn)4-6

1、實(shí)驗(yàn)4 離散系統(tǒng)的變換域分析一、實(shí)驗(yàn)?zāi)康?、熟悉對(duì)離散系統(tǒng)的頻率響應(yīng)分析方法；2、加深對(duì)零、極點(diǎn)分布的概念理解。 二、實(shí)驗(yàn)原理離散系統(tǒng)的時(shí)域方程為其變換域分析方法如下：頻域:系統(tǒng)的頻率響應(yīng)為:Z域:系統(tǒng)的轉(zhuǎn)移函數(shù)為:分解因式:，其中和稱為零、極點(diǎn)。三、預(yù)習(xí)要求1. 在MATLAB中，熟悉函數(shù)tf2zp、zplane、freqz、residuez、zp2sos的使用，其中：z，p，K=tf2zp（num，den）求得有理分式形式的系統(tǒng)轉(zhuǎn)移函數(shù)的零、極點(diǎn)；zplane（z，p）繪制零、極點(diǎn)分布圖；h=freqz(num,den,w)求系統(tǒng)的單位頻率響應(yīng)；r，p，k=residuez（num，den

2、）完成部分分式展開(kāi)計(jì)算；sos=zp2sos（z，p，K）完成將高階系統(tǒng)分解為2階系統(tǒng)的串聯(lián)。 2. 閱讀擴(kuò)展練習(xí)中的實(shí)例，學(xué)習(xí)頻率分析法在MATLAB中的實(shí)現(xiàn)； 3. 編程實(shí)現(xiàn)系統(tǒng)參數(shù)輸入，繪出幅度頻率響應(yīng)和相位響應(yīng)曲線和零、極點(diǎn)分布圖。 四、實(shí)驗(yàn)內(nèi)容求系統(tǒng)的零、極點(diǎn)和幅度頻率響應(yīng)和相位響應(yīng)。解析：【代碼】num=0.0528 0.0797 0.1295 0.1295 0.797 0.0528;den=1 -1.8107 2.4947 -1.8801 0.9537 -0.2336;z,p,k=tf2zp(num,den);disp('零點(diǎn)');disp(z);disp(

3、9;極點(diǎn)');disp(p);disp('增益系數(shù)');disp(k);figure(1)zplane(num,den)figure(2)freqz(num,den,128) 【圖形】【結(jié)果】零點(diǎn) -1.5870 + 1.4470i -1.5870 - 1.4470i 0.8657 + 1.5779i 0.8657 - 1.5779i -0.0669 + 0.0000i極點(diǎn) 0.2788 + 0.8973i 0.2788 - 0.8973i 0.3811 + 0.6274i 0.3811 - 0.6274i 0.4910 + 0.0000i增益系數(shù) 0.0528五、擴(kuò)展

4、練習(xí)例1: 求下列直接型系統(tǒng)函數(shù)的零、極點(diǎn)，并將它轉(zhuǎn)換成二階節(jié)形式解:用MATLAB計(jì)算程序如下：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”的分別為分子和分

5、母多項(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.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.15520 0.65111.0000 0.28850 0.36300 1.0000 -1.0552 0.7679系統(tǒng)函數(shù)的二階節(jié)形式為：極點(diǎn)圖見(jiàn)圖:圖 系統(tǒng)函數(shù)的零、極點(diǎn)圖例2: 差分方程所對(duì)應(yīng)的系統(tǒng)的頻率響應(yīng)。解:差分方程所對(duì)應(yīng)的系統(tǒng)函

6、數(shù)為：用MATLAB計(jì)算的程序如下：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('實(shí)部')xlabel('omega/pi');ylabel('幅度')subplot(2,2,2);plot(w/pi,imag(h);gridtitle('虛部')xlabel('omega/pi');ylabel('Am

7、plitude')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('弧度') 圖實(shí)驗(yàn)5 有限沖激響應(yīng)數(shù)字濾波器設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康模?1、加深對(duì)數(shù)字濾波器的常用指標(biāo)理解。 2、學(xué)習(xí)數(shù)字濾波器的設(shè)計(jì)方法。 二、實(shí)驗(yàn)原理：低通濾波器的常用指標(biāo)

8、：（1）通帶邊緣頻率；（2）阻帶邊緣頻率；（3）通帶起伏；（4）通帶峰值起伏，（5）阻帶起伏，最小阻帶衰減。三、預(yù)習(xí)要求1. 在MATLAB中，熟悉函數(shù)fir1、kaiserord 、remezord、remez的使用； B = fir1(n,Wn,'high','noscale')設(shè)計(jì)濾波器； n,Wn,beta,ftype = kaiserord(f,a,dev)估計(jì)濾波器階數(shù)； n,fo,ao,w = remezord (f,a,dev,fs)計(jì)算等波紋濾波器階數(shù)n和加權(quán)函數(shù)w()； B=remez(n,f,a)進(jìn)行等波紋濾波器的設(shè)計(jì)。 2. 閱讀擴(kuò)展練習(xí)

9、中的實(shí)例，學(xué)習(xí)FIR濾波器的設(shè)計(jì)方法及其在MATLAB中的實(shí)現(xiàn)； 3. 給出FIR數(shù)字濾波器的沖激響應(yīng)，繪出它們的幅度和相位頻響曲線，討論它們各自的實(shí)現(xiàn)形式和特點(diǎn)。 數(shù)字濾波器有IIR和FIR兩種類型，它們的特點(diǎn)和設(shè)計(jì)方法不同。四、實(shí)驗(yàn)內(nèi)容：利用MATLAB編程,分別用窗函數(shù)法和等波紋濾波器法設(shè)計(jì)兩種FIR數(shù)字濾波器，指標(biāo)要求如下：通帶邊緣頻率：，通帶峰值起伏：。阻帶邊緣頻率：，最小阻帶衰減：。 l 用窗函數(shù)法設(shè)計(jì)： 【解析】調(diào)用函數(shù)n,wn,bta,ftype=kaiserord（f，a，dev，fs）參數(shù)：f=0.3 0.45 0.65 0.8為對(duì)應(yīng)

10、數(shù)字頻率a=0 1 0為由f指定的各個(gè)頻帶上的幅值向量，一般只有0和1表示；和f長(zhǎng)度關(guān)系為（2*a的長(zhǎng)度）2=（f的長(zhǎng)度） devs=0.01 0.1087 0.01用于指定各個(gè)頻帶輸出濾波器的頻率響應(yīng)與其期望幅值之間的最大輸出誤差或偏差，長(zhǎng)度與a相等，計(jì)算公式：阻帶衰減誤差=，通帶波動(dòng)衰減誤差=fs缺省值為2HZ【代碼】n,wn,bta,ftype=kaiserord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用kaiserord函數(shù)估計(jì)出濾波器階數(shù)n和beta參數(shù)h1=fir1(n,wn,ftyp

11、e,kaiser(n+1,bta),'noscale');hh1,w1=freqz(h1,1,256);figure(1)subplot(2,1,1)plot(w1/pi,20*log10(abs(hh1)gridxlabel('歸一化頻率w');ylabel('幅度/db')subplot(2,1,2)plot(w1/pi,angle(hh1)gridxlabel('歸一化頻率w');ylabel('相位/rad');【圖形】l 用等波紋法設(shè)計(jì)： 【解析】調(diào)用函數(shù)n,fpts,mag,wt=remezo

12、rd(f，a，dev) f=0.3 0.45 0.65 0.8 a=0 1 0 dev=0.01 0.1087 0.01 其含義同函數(shù)n,wn,bta,ftype=kaiserord（f，a，dev，fs）中的參數(shù)相同?！敬a】n,fpts,mag,wt=remezord(0.3 0.45 0.65 0.8,0 1 0,0.01 0.1087 0.01);%用remezord函數(shù)估算出remez函數(shù)要用到的階n、歸一化頻帶邊緣矢量fpts、頻帶內(nèi)幅值響應(yīng)矢量mag及加權(quán)矢量w，使rem

13、ez函數(shù)設(shè)計(jì)出的濾波器滿足f、a及dev指定的性能要求h2=remez(n,fpts,mag,wt);%設(shè)計(jì)出等波紋濾波器hh2,w2=freqz(h2,1,256);figure(2)subplot(2,1,1)plot(w2/pi,20*log10(abs(hh2)gridxlabel('歸一化頻率w');ylabel('幅度/db')subplot(2,1,2)plot(w2/pi,angle(hh2)gridxlabel('歸一化頻率w');ylabel('相位/rad');【圖形】五、擴(kuò)展練習(xí)例1: 用凱塞窗設(shè)計(jì)一FIR

14、低通濾波器，通帶邊界頻率，阻帶邊界頻率，阻帶衰減不小于50dB。解: 首先由過(guò)渡帶寬和阻帶衰減來(lái)決定凱塞窗的N和上圖給出了以上設(shè)計(jì)的頻率特性，(a) 為N=30直接截取的頻率特性(b)為凱塞窗設(shè)計(jì)的頻率特性。凱塞窗設(shè)計(jì)對(duì)應(yīng)的MATLAB程序?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);或者：b = fir1(29,0.4,kaiser(30,4.55);h1,w1=freqz(b,1);plot(w1/pi,20*

15、log10(abs(h1);axis(0,1,-80,10); grid;xlabel('歸一化頻率/p');ylabel('幅度/dB');還可以使用n,Wn,beta,ftype=kaiserord(f,a,dev)函數(shù)來(lái)估計(jì)濾波器階數(shù)等，得到凱塞窗濾波器：fcuts = 0.3 0.5; %歸一化頻率omega/pimags = 1 0;devs = 0.05 10(-2.5);n,Wn,beta,ftype = kaiserord(fcuts,mags,devs); %計(jì)算出凱塞窗N，beta的值hh = fir1(n,Wn,ftype,kaiser(n

16、+1,beta),'noscale');freqz(hh);實(shí)際中，一般調(diào)用MATLAB信號(hào)處理工具箱函數(shù)remezord來(lái)計(jì)算等波紋濾波器階數(shù)N和加權(quán)函數(shù)W()，調(diào)用函數(shù)remez可進(jìn)行等波紋濾波器的設(shè)計(jì)，直接求出濾波器系數(shù)。函數(shù)remezord中的數(shù)組fedge為通帶和阻帶邊界例2：利用雷米茲交替算法設(shè)計(jì)等波紋濾波器，設(shè)計(jì)一個(gè)線性相位低通FIR數(shù)字濾波器，其指標(biāo)為：通帶邊界頻率fc=800Hz，阻帶邊界fr=1000Hz，通帶波動(dòng)阻帶最小衰減At=40dB，采樣頻率fs=4000Hz。解:在MATLAB中可以用remezord 和remez兩個(gè)函數(shù)設(shè)計(jì)，其結(jié)果如圖2，MAT

17、LAB程序如下：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í)驗(yàn)6 無(wú)限沖激響應(yīng)數(shù)字濾波器設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康?、掌握雙線性變換法及脈沖相應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的具體設(shè)計(jì)方

18、法;2、熟悉用雙線性變換法及脈沖響應(yīng)不變法設(shè)計(jì)低通、高通和帶通IIR數(shù)字濾波器的計(jì)算機(jī)編程。二、實(shí)驗(yàn)原理在MATLAB中，可以用下列函數(shù)輔助設(shè)計(jì)IIR數(shù)字濾波器：1)利用buttord和cheb1ord可以確定低通原型巴特沃斯和切比雪夫?yàn)V波器的階數(shù)和截止頻率；2)num,den=butter（N,Wn）(巴特沃斯)和num,den=cheby1（N,Wn）,num,den=cheby2（N,Wn）(切比雪夫1型和2型)可以進(jìn)行濾波器的設(shè)計(jì)；3）lp2hp，lp2bp，lp2bs可以完成低通濾波器到高通、帶通、帶阻濾波器的轉(zhuǎn)換；4)使用bilinear可以對(duì)模擬濾波器進(jìn)行雙線性變換，求得數(shù)字濾波

19、器的傳輸函數(shù)系數(shù)；5)利用impinvar可以完成脈沖響應(yīng)不變法的模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換。三、預(yù)習(xí)要求1. 在MATLAB中，熟悉函數(shù)butter、cheby1、cheby2的使用，其中：num,den=butter（N,Wn）巴特沃斯濾波器設(shè)計(jì)；num,den=cheby1（N,Wn）切比雪夫1型濾波器設(shè)計(jì)；num,den=cheby2（N,Wn）切比雪夫2型濾波器設(shè)計(jì)。2. 閱讀擴(kuò)展練習(xí)中的實(shí)例，學(xué)習(xí)在MATLAB中進(jìn)行數(shù)字濾波器的設(shè)計(jì)； 3. 給出IIR數(shù)字濾波器參數(shù)和濾波器的沖激響應(yīng)，繪出它們的幅度和相位頻響曲線，討論它們各自的實(shí)現(xiàn)形式和特點(diǎn)。 四、實(shí)驗(yàn)內(nèi)容利用MATLAB編程，

20、用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)一個(gè)數(shù)字帶通濾波器，指標(biāo)要求如下：通帶邊緣頻率：，通帶峰值起伏：；阻帶邊緣頻率：，最小阻帶衰減： ?！窘馕觥縧 巴特沃茲原型1 雙線性變換法【代碼】ws1=2*8000*tan(0.3*pi/2);ws2=2*8000*tan(0.8*pi/2);wp1=2*8000*tan(0.45*pi/2);wp2=2*8000*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,

21、A=bilinear(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1);plot(f,20*log10(abs(h);axis(0,4000,-60,10);gridxlabel('頻率/Hz');ylabel('幅度/dB');subplot(2,1,2)plot(f,angle(h);gridxlabel('頻率/Hz');ylabel('相位');【圖形】2 脈沖響應(yīng)不變法【代碼】fs=8000;ws1=0.3*pi*fs;ws2=0.8*pi*fs;wp1=0.

22、45*pi*fs;wp2=0.65*pi*fs;ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=impinvar(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1)plot(f,20*log10(abs(h);axis(0,4000,-80,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB');subplot(

23、2,1,2)plot(f,angle(h)grid;xlabel('頻率/Hz');ylabel('相位');【圖形】fs=2;ws1=2*fs*tan(0.3*pi/2);ws2=2*fs*tan(0.8*pi/2);wp1=2*fs*tan(0.45*pi/2);wp2=2*fs*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=buttord(wp,ws,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=bilinear(num,den,f

24、s);h1,w=freqz(B,A);ws11=0.3*pi*fs;ws22=0.8*pi*fs;wp11=0.45*pi*fs;wp22=0.65*pi*fs;ws0=ws11 ws22;wp0=wp11 wp22; N,Wn=buttord(wp0,ws0,Rp,Rs,'s');num,den=butter(N,Wn,'s');B,A=impinvar(num,den,fs);h2,w=freqz(B,A);figure(1)plot(w/pi,20*log10(abs(h1),'-.',w/pi,20*log10(abs(h2),'

25、-');axis(0,1,-80,0);grid;xlabel('頻率/Hz');ylabel('幅度/dB');l 切必雪夫原型1 雙線性變換法【代碼】ws1=2*8000*tan(0.3*pi/2);ws2=2*8000*tan(0.8*pi/2);wp1=2*8000*tan(0.45*pi/2);wp2=2*8000*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2;Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s&

26、#39;);B,A=bilinear(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1)plot(f,20*log10(abs(h);axis(0,4000,-60,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB');subplot(2,1,2)plot(f,angle(h);grid;xlabel('頻率/Hz');ylabel('相位');【圖形】2 脈沖響應(yīng)不變法【代碼】fs=8000;ws1=0.3*pi*fs;ws2=0.8*pi*

27、fs;wp1=0.45*pi*fs;wp2=0.65*pi*fs;ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s');B,A=impinvar(num,den,8000);h,w=freqz(B,A);f=w/pi*4000;subplot(2,1,1);plot(f,20*log10(abs(h);axis(0,4000,-90,10);grid;xlabel('頻率/Hz');ylabel('幅度/dB&

28、#39;);subplot(2,1,2);plot(f,angle(h);grid;xlabel('頻率/Hz');ylabel('相位');【圖形】fs=8000;ws1=2*fs*tan(0.3*pi/2);ws2=2*fs*tan(0.8*pi/2);wp1=2*fs*tan(0.45*pi/2);wp2=2*fs*tan(0.65*pi/2);ws=ws1 ws2;wp=wp1 wp2; Rp=1;Rs=40;N,Wn=cheb1ord(wp,ws,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s')

29、;B,A=bilinear(num,den,fs);h1,w=freqz(B,A);ws11=0.3*pi*fs;ws22=0.8*pi*fs;wp11=0.45*pi*fs;wp22=0.65*pi*fs;ws0=ws11 ws22;wp0=wp11 wp22;N,Wn=cheb1ord(wp0,ws0,Rp,Rs,'s');num,den=cheby1(N,1,Wn,'s');B,A=impinvar(num,den,fs);h2,w=freqz(B,A);figure(1)plot(w/pi,20*log10(abs(h1),'-.',w/

30、pi,20*log10(abs(h2),'-');axis(0,1,-80,10);grid;xlabel('頻率/Hz')ylabel('幅度/dB');【結(jié)論】 雙線性變換法通過(guò)將數(shù)字頻率的取值范圍從0到p對(duì)應(yīng)到模擬頻率的范圍0到的范圍，也就對(duì)應(yīng)于模擬域中所有可能的頻率值。雙線性變換法不會(huì)出現(xiàn)頻率混疊，但非線性關(guān)系卻導(dǎo)致數(shù)字濾波器的頻率響應(yīng)不能逼真地模仿模擬濾波器的頻率響應(yīng)。脈沖響應(yīng)不變法通過(guò)選擇滿足設(shè)計(jì)要求的模擬濾波器沖激響應(yīng)h（t）的采樣值的數(shù)字脈沖響應(yīng)hn得到的被采樣的沖激響應(yīng)將給出與原模擬濾波器非常相近的濾波器形狀。由于該方法不可避免

31、的要發(fā)生頻率混疊現(xiàn)象，所以只適合設(shè)計(jì)低通和帶通濾波器。從實(shí)驗(yàn)結(jié)果可以看出：雙線性變換法所設(shè)計(jì)的巴特沃茲濾波器最符合設(shè)計(jì)指標(biāo)，而用脈沖響應(yīng)不變法設(shè)計(jì)的濾波器（無(wú)論是巴特沃茲還是切比雪夫）有一定的誤差，主要是由于混疊所引起。五、擴(kuò)展練習(xí)例1：設(shè)采樣周期T=250s（采樣頻率fs =4kHz），用脈沖響應(yīng)不變法和雙線性變換法設(shè)計(jì)一個(gè)三階巴特沃茲濾波器，其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.0002

32、5,'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×1000，為脈沖響應(yīng)不變法原型低通濾波器的邊界頻率；第二個(gè)butter的邊界頻率2/T=2/0.00025，為雙線性變換法原型低通濾波器的邊界頻率.圖1給出了這兩種設(shè)計(jì)方法所得到的頻響，虛線為脈沖響應(yīng)不變法的結(jié)果；實(shí)線為雙線性變換法的結(jié)果。脈沖響應(yīng)不變法由于混疊效應(yīng)，使得過(guò)渡帶和阻帶的衰減特性變差，并且不存在傳輸零點(diǎn)。同時(shí)，也看到雙線性變換法，在z=-1即=或f=2000Hz處有一個(gè)三階傳輸零點(diǎn)，這個(gè)三階零點(diǎn)正是模擬濾波器在=處的三階傳輸零點(diǎn)通過(guò)映射形成的。下圖給出了MATLAB計(jì)算的結(jié)果。 例2: