實驗六 用雙線性變換法設(shè)計IIR數(shù)字濾波器

1、實驗六 用雙線性變換法設(shè)計IIR數(shù)字濾波器一、實驗?zāi)康膌 學會運用MATLAB設(shè)計數(shù)字低通、帶通、高通、帶阻濾波器的設(shè)計方法。二、實驗涉及的matlab子函數(shù)bilinear功能：雙線性變換將s域映射到z域。調(diào)用格式：numd,dend= bilinear (num,den,Fs)，將模擬域系統(tǒng)函數(shù)轉(zhuǎn)換為數(shù)字域的系統(tǒng)函數(shù)，F(xiàn)s為采樣頻率。三、實驗原理下面舉例說明用雙線性變換法設(shè)計各種數(shù)字濾波器的過程。例1、 采用雙線性變換法設(shè)計一個巴特沃斯數(shù)字低通濾波器，要求：wp=0.25*pi，rp=1db，ws=0.4*pi，as=15db，濾波器采樣頻率Fs=100hz。MATLAB源程序為：%數(shù)字濾

2、波器指標wpd=0.25*pi; %濾波器的通帶截止頻率wsd=0.4*pi; %濾波器的阻帶截止頻率Rp=1;As=15; %輸入濾波器的通阻帶衰減指標%轉(zhuǎn)換為模擬原型濾波器指標Fs=100;T=1/Fs;wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/2);%模擬原型濾波器計算n,wc=buttord(wp,ws,Rp,As,'s') %計算階數(shù)n和截止頻率z0,p0,k0=buttap(n); %歸一化切比雪夫1型原型設(shè)計ba=k0*poly(z0); %求原型濾波器系數(shù)baa=poly(p0); %求原型濾波器系數(shù)aba1,aa1=lp2lp(

3、ba,aa,wc);%變換為模擬低通濾波器%用雙線性變換法計算數(shù)字濾波器系數(shù)bd,ad=bilinear(ba1,aa1,Fs) %雙線性變換%求數(shù)字系統(tǒng)的頻率特性H,w=freqz(bd,ad); dbH=20*log10(abs(H)/max(abs(H); %化為分貝值subplot(2,2,1),plot(w,abs(H); ylabel('|H|');title('幅度響應(yīng)');axis(0,pi,0,1.1);gridsubplot(2,2,2),plot(w,angle(H);ylabel('phi');title('相位響

4、應(yīng)');axis(0,pi,-4,4);gridsubplot(2,2,3),plot(w,dbH);title('幅度響應(yīng)(dB)'); ylabel('dB');xlabel('頻率');axis(0,pi,-40,5);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零極圖');運行結(jié)果為：n = 5wc = 103.2016bd = 0.0072 0.0362 0.0725 0.0725 0.0362 0.0072ad = 1.0000

5、 -1.9434 1.9680 -1.0702 0.3166 -0.0392則所求濾波器的系統(tǒng)函數(shù)為例2、采用雙線性變換法設(shè)計一個橢圓數(shù)字高通濾波器，要求通帶250hz，1db，阻帶150hz，20db，濾波器采樣頻率為Fs=1000hz。MATLAB源程序為：%數(shù)字濾波器指標fs=150;fp=250;Fs=1000;T=1/Fs;wpd=fp/Fs*2*pi; %數(shù)字濾波器的通帶截止頻率wsd=fs/Fs*2*pi; %數(shù)字濾波器的阻帶截止頻率Rp=1;As=20; %輸入濾波器的通阻帶衰減指標%轉(zhuǎn)換為模擬濾波器指標wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/

6、2);%模擬原型濾波器計算n,wc=ellipord(wp,ws,Rp,As,'s') %計算階數(shù)n和截止頻率z0,p0,k0=ellipap(n,Rp,As); %歸一化橢圓原型設(shè)計ba=k0*poly(z0); %求原型濾波器系數(shù)baa=poly(p0); %求原型濾波器系數(shù)aba1,aa1=lp2hp(ba,aa,wc);%變換為模擬高通濾波器%用雙線性變換法計算數(shù)字濾波器系數(shù)bd,ad=bilinear(ba1,aa1,Fs) %雙線性變換%求數(shù)字系統(tǒng)的頻率特性H,w=freqz(bd,ad); dbH=20*log10(abs(H)/max(abs(H); %化為分貝

7、值%subplot(2,2,1),plot(w/2/pi*Fs,abs(H); ylabel('|H|');title('幅度響應(yīng)');axis(0,Fs/2,0,1.1);gridsubplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180);ylabel('phi');title('相位響應(yīng)');axis(0,Fs/2,-180,180);gridsubplot(2,2,3),plot(w/2/pi*Fs,dbH);title('幅度響應(yīng)(dB)');axis(0,Fs/2,-40

8、,5); ylabel('dB');xlabel('頻率(hz)');gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零極圖');運行結(jié)果為n = 3wc = 2.0000e+003bd = 0.2545 -0.4322 0.4322 -0.2545ad = 1.0000 0.1890 0.7197 0.1574例3、采用雙線性變換法設(shè)計一個切比雪夫1型數(shù)字帶通濾波器，要求：通帶0.3pi0.7pi，1db，阻帶0.2pi，0.8pi，20db，濾波器采樣周期為Ts=

9、0. 1s。MATLAB源程序為：%雙線性變換法設(shè)計數(shù)字帶通%數(shù)字濾波器指標wpd1=0.3*pi;wpd2=0.7*pi; %數(shù)字濾波器的通帶截止頻率wsd1=0.2*pi;wsd2=0.8*pi; %數(shù)字濾波器的阻帶截止頻率Rp=1;As=20; %輸入濾波器的通阻帶衰減指標%轉(zhuǎn)換為模擬濾波器指標Fs=10;T=1/Fs;wp1=(2/T)*tan(wpd1/2);wp2=(2/T)*tan(wpd2/2);wp=wp1,wp2; %模擬濾波器的通帶截止頻率ws1=(2/T)*tan(wsd1/2);ws2=(2/T)*tan(wsd2/2);ws=ws1,ws2; %模擬濾波器的阻帶截止

10、頻率bw=wp2-wp1;w0=sqrt(wp1*wp2);%模擬通帶帶寬和中心頻率%模擬原型濾波器計算n,wn=cheb1ord(wp,ws,Rp,As,'s') %計算階數(shù)n和截止頻率z0,p0,k0=cheb1ap(n,Rp); %設(shè)計歸一化的模擬原型濾波器ba1=k0*poly(z0); %求原型濾波器系數(shù)baa1=poly(p0); %求原型濾波器系數(shù)aba,aa=lp2bp(ba1,aa1,w0,bw); %變換為模擬帶通濾波器%用雙線性變換法計算數(shù)字濾波器系數(shù)bd,ad=bilinear(ba,aa,Fs)%求數(shù)字系統(tǒng)的頻率特性H,w=freqz(bd,ad);

11、dbH=20*log10(abs(H)/max(abs(H); %化為分貝值%subplot(2,2,1),plot(w/pi,abs(H),'k'); ylabel('幅度');xlabel('頻率/pi');axis(0,1,0,1.1);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi,'k');ylabel('相位');xlabel('頻率/pi');axis(0,1,-1,1);gridsubplot(2,2,3),plot(w/pi,dbH,'k&

12、#39;); ylabel('幅度(dB)');xlabel('頻率/pi');axis(0,1,-60,5);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);ylabel('零極圖');運行結(jié)果為n = 3wn = 10.1905 39.2522bd = 0.0736 0.0000 -0.2208 -0.0000 0.2208 -0.0000 -0.0736ad = 1.0000 -0.0000 0.9761 0.0000 0.8568 0.0000 0.2919例4、采用雙線性

13、變換法設(shè)計一個切比雪夫1型數(shù)字帶阻濾波器，要求：阻帶0.3pi0.7pi，20db，通帶0.2pi，0.8pi，1db，濾波器采樣周期為Ts=0. 1s。MATLAB源程序為：%數(shù)字濾波器指標wsd1=0.3*pi;wsd2=0.7*pi; %數(shù)字濾波器的通帶截止頻率wpd1=0.2*pi;wpd2=0.8*pi; %數(shù)字濾波器的阻帶截止頻率Rp=1;As=20; %輸入濾波器的通阻帶衰減指標%轉(zhuǎn)換為模擬濾波器指標Fs=10;T=1/Fs;wp1=(2/T)*tan(wpd1/2);wp2=(2/T)*tan(wpd2/2);wp=wp1,wp2; %模擬濾波器的通帶截止頻率ws1=(2/T)

14、*tan(wsd1/2);ws2=(2/T)*tan(wsd2/2);ws=ws1,ws2; %模擬濾波器的阻帶截止頻率bw=wp2-wp1;w0=sqrt(wp1*wp2);%模擬通帶帶寬和中心頻率%模擬原型濾波器計算n,wn=cheb1ord(wp,ws,Rp,As,'s') %計算階數(shù)n和截止頻率z0,p0,k0=cheb1ap(n,Rp); %設(shè)計歸一化的模擬原型濾波器ba1=k0*poly(z0); %求原型濾波器系數(shù)baa1=poly(p0); %求原型濾波器系數(shù)aba,aa=lp2bs(ba1,aa1,w0,bw); %變換為模擬帶阻濾波器%用雙線性變換法計算數(shù)字

15、濾波器系數(shù)bd,ad=bilinear(ba,aa,Fs)%求數(shù)字系統(tǒng)的頻率特性H,w=freqz(bd,ad); dbH=20*log10(abs(H)/max(abs(H); %化為分貝值%subplot(2,2,1),plot(w/pi,abs(H),'k'); ylabel('幅度');xlabel('頻率/pi');axis(0,1,0,1.1);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi,'k');ylabel('相位');xlabel('頻率/pi'

16、;);axis(0,1,-1,1);gridsubplot(2,2,3),plot(w/pi,dbH,'k'); ylabel('幅度(dB)');xlabel('頻率/pi');axis(0,1,-60,5);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);ylabel('零極圖');運行結(jié)果為n = 3wn = 6.4985 61.5532bd = 0.0736 0.0000 0.2208 0.0000 0.2208 0.0000 0.0736ad = 1.00

17、00 0.0000 -0.9761 -0.0000 0.8568 0.0000 -0.2919四、實驗任務(wù)1、采用雙線性變換法設(shè)計一個巴特沃斯數(shù)字高通濾波器，要求通帶0.35pi，1db，阻帶0.2pi，15db，濾波器采樣頻率為Fs=10hz。列出系統(tǒng)函數(shù)并做頻率響應(yīng)曲線和零極點分布圖。Fs=1000;T=1/Fs;wp=0.35*pi;ws=0.2*pi;fp=wp/(2*pi)*Fs;fs=ws/(2*pi)*Fs;Rp=1;As=15;ripple=10(-Rp/20);Attn=10(-As/20);Omgp=(2/T)*tan(wp/2);Omgs=(2/T)*tan(ws/2);

18、n,Omgc=cheb2ord(Omgp,Omgs,Rp,As,'s') z0,p0,k0=cheb2ap(n,As);ba=k0*real(poly(z0);aa=real(poly(p0);ba1,aa1=lp2hp(ba,aa,Omgc);bd,ad=bilinear(ba1,aa1,Fs)H,w=freqz(bd,ad);dbH=20*log10(abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');ylabel('|H|');title('幅度響應(yīng)

19、9;);axis(0,Fs/2,0,1.1);set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2);set(gca,'YTickMode','manual','YTick',0,Attn,ripple,1);gridsubplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180,'k');ylabel('phi');title('相位響應(yīng)');axis(0,Fs/2,-180,

20、180);set(gca,'XTickMode','manual','XTick',0,fs,fp,Fs/2);set(gca,'YTickMode','manual','YTick',-180,0,180);gridsubplot(2,2,3),plot(w/2/pi*Fs,dbH);title('幅度響應(yīng)( dB)');axis(0,Fs/2,-40,5);ylabel('dB');xlabel('頻率(pi)');set(gca,'XTi

21、ckMode','manual','XTick',0,fs,fp,Fs/2);set(gca,'YTickMode','manual','YTick',-50,-20,-1,0);gridsubplot(2,2,4),zplane(bd,ad);axis(-1.1,1.1,-1.1,1.1);title('零極圖');2、采用雙線性變換法設(shè)計一個橢圓型數(shù)字帶阻濾波器，要求：阻帶0.4pi0.6pi，20db，通帶0.35pi，0.65pi，1db，濾波器采樣周期為Ts=0. 1s。列出系統(tǒng)函

22、數(shù)并做頻率響應(yīng)曲線和零極點分布圖。wp1=0.35*pi;wp2=0.65*pi;ws1=0.4*pi;ws2=0.6*pi;fp=wp/(2*pi)*Fs;fs=ws/(2*pi)*Fs;Rp=1;As=20;T=0.1;Fs=1/T;Omgp1=(2/T)*tan(wp1/2);Omgp2=(2/T)*tan(wp2/2);Omgp=Omgp1,Omgp2;Omgs1=(2/T)*tan(ws1/2);Omgs2=(2/T)*tan(ws2/2);Omgs=Omgs1,Omgs2;bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2);n,Omgn=ellipord(Omg

23、p,Omgs,Rp,As,'s')z0,p0,k0=ellipap(n,Rp,As);ba1=k0*real(poly(z0);aa1=real(poly(p0);ba,aa=lp2bs(ba1,aa1,w0,bw);bd,ad=bilinear(ba,aa,Fs)H,w=freqz(bd,ad);dbH=20*log10(abs(H)+eps)/max(abs(H);subplot(2,2,1),plot(w/2/pi*Fs,abs(H),'k');ylabel('|H|');title('幅度響應(yīng)');axis(0,Fs/2,0,1.1);set(gca,'XTickMo