實(shí)驗(yàn)六-用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器

實(shí)驗(yàn)六用雙線性變換法設(shè)計(jì)IIR數(shù)字濾波器一、實(shí)驗(yàn)?zāi)康膶W(xué)會(huì)運(yùn)用MATLAB設(shè)計(jì)數(shù)字低通、帶通、高通、帶阻濾波器的設(shè)計(jì)方法。二、實(shí)驗(yàn)涉及的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í)驗(yàn)原理下面舉例說(shuō)明用雙線性變換法設(shè)計(jì)各種數(shù)字濾波器的過(guò)程。采用雙線性變換法設(shè)計(jì)一個(gè)巴特沃斯數(shù)字低通濾波器，要求：wp=0.25*pi，rp=1db，ws=0.4*pi，as=15db，濾波器采樣頻率Fs=100hz。MATLAB源程序?yàn)椋?數(shù)字濾波器指標(biāo)wpd=0.25*pi;%濾波器的通帶截止頻率wsd=0.4*pi;%濾波器的阻帶截止頻率Rp=1;As=15;%輸入濾波器的通阻帶衰減指標(biāo)%轉(zhuǎn)換為模擬原型濾波器指標(biāo)Fs=100;T=1/Fs;wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/2);%模擬原型濾波器計(jì)算[n,wc]=buttord(wp,ws,Rp,As,'s')%計(jì)算階數(shù)n和截止頻率[z0,p0,k0]=buttap(n);%歸一化切比雪夫1型原型設(shè)計(jì)ba=k0*poly(z0);%求原型濾波器系數(shù)baa=poly(p0);%求原型濾波器系數(shù)a[ba1,aa1]=lp2lp(ba,aa,wc);%變換為模擬低通濾波器%用雙線性變換法計(jì)算數(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('相位響應(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('零極圖');運(yùn)行結(jié)果為：n=5wc=103.2023bd=0.00720.03620.07250.07250.03620.0072ad=1.0000-1.94341.9680-1.07020.3166-0.0392那么所求濾波器的系統(tǒng)函數(shù)為例2、采用雙線性變換法設(shè)計(jì)一個(gè)橢圓數(shù)字高通濾波器，要求通帶250hz，1db，阻帶150hz，20db，濾波器采樣頻率為Fs=1000hz。MATLAB源程序?yàn)椋?數(shù)字濾波器指標(biāo)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;%輸入濾波器的通阻帶衰減指標(biāo)%轉(zhuǎn)換為模擬濾波器指標(biāo)wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/2);%模擬原型濾波器計(jì)算[n,wc]=ellipord(wp,ws,Rp,As,'s')%計(jì)算階數(shù)n和截止頻率[z0,p0,k0]=ellipap(n,Rp,As);%歸一化橢圓原型設(shè)計(jì)ba=k0*poly(z0);%求原型濾波器系數(shù)baa=poly(p0);%求原型濾波器系數(shù)a[ba1,aa1]=lp2hp(ba,aa,wc);%變換為模擬高通濾波器%用雙線性變換法計(jì)算數(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/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,5]);ylabel('dB');xlabel('頻率(hz)');gridsubplot(2,2,4),zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title('零極圖');運(yùn)行結(jié)果為n=3wc=2.0000e+003bd=0.2545-0.43220.4322-0.2545ad=1.00000.18900.71970.1574例3、采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫1型數(shù)字帶通濾波器，要求：通帶0.3pi~0.7pi，1db，阻帶0.2pi，0.8pi，20db，濾波器采樣周期為Ts=0.1s。MATLAB源程序?yàn)椋?雙線性變換法設(shè)計(jì)數(shù)字帶通%數(shù)字濾波器指標(biāo)wpd1=0.3*pi;wpd2=0.7*pi;%數(shù)字濾波器的通帶截止頻率wsd1=0.2*pi;wsd2=0.8*pi;%數(shù)字濾波器的阻帶截止頻率Rp=1;As=20;%輸入濾波器的通阻帶衰減指標(biāo)%轉(zhuǎn)換為模擬濾波器指標(biāo)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];%模擬濾波器的阻帶截止頻率bw=wp2-wp1;w0=sqrt(wp1*wp2);%模擬通帶帶寬和中心頻率%模擬原型濾波器計(jì)算[n,wn]=cheb1ord(wp,ws,Rp,As,'s')%計(jì)算階數(shù)n和截止頻率[z0,p0,k0]=cheb1ap(n,Rp);%設(shè)計(jì)歸一化的模擬原型濾波器ba1=k0*poly(z0);%求原型濾波器系數(shù)baa1=poly(p0);%求原型濾波器系數(shù)a[ba,aa]=lp2bp(ba1,aa1,w0,bw);%變換為模擬帶通濾波器%用雙線性變換法計(jì)算數(shù)字濾波器系數(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');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('零極圖');運(yùn)行結(jié)果為n=3wn=10.190539.2522bd=0.07360.0000-0.2208-0.00000.2208-0.0000-0.0736ad=1.0000-0.00000.97610.00000.85680.00000.2919例4、采用雙線性變換法設(shè)計(jì)一個(gè)切比雪夫1型數(shù)字帶阻濾波器，要求：阻帶0.3pi~0.7pi，20db，通帶0.2pi，0.8pi，1db，濾波器采樣周期為Ts=0.1s。MATLAB源程序?yàn)椋?數(shù)字濾波器指標(biāo)wsd1=0.3*pi;wsd2=0.7*pi;%數(shù)字濾波器的通帶截止頻率wpd1=0.2*pi;wpd2=0.8*pi;%數(shù)字濾波器的阻帶截止頻率Rp=1;As=20;%輸入濾波器的通阻帶衰減指標(biāo)%轉(zhuǎn)換為模擬濾波器指標(biāo)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];%模擬濾波器的阻帶截止頻率bw=wp2-wp1;w0=sqrt(wp1*wp2);%模擬通帶帶寬和中心頻率%模擬原型濾波器計(jì)算[n,wn]=cheb1ord(wp,ws,Rp,As,'s')%計(jì)算階數(shù)n和截止頻率[z0,p0,k0]=cheb1ap(n,Rp);%設(shè)計(jì)歸一化的模擬原型濾波器ba1=k0*poly(z0);%求原型濾波器系數(shù)baa1=poly(p0);%求原型濾波器系數(shù)a[ba,aa]=lp2bs(ba1,aa1,w0,bw);%變換為模擬帶阻濾波器%用雙線性變換法計(jì)算數(shù)字濾波器系數(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');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('零極圖');運(yùn)行結(jié)果為n=3wn=6.498561.5532bd=0.07360.00000.22080.00000.22080.00000.0736ad=1.00000.0000-0.9761-0.00000.85680.0000-0.2919四、實(shí)驗(yàn)任務(wù)1、采用雙線性變換法設(shè)計(jì)一個(gè)巴特沃斯數(shù)字高通濾波器，要求通帶0.35pi，1db，阻帶0.2pi，15db，濾波器采樣頻率為Fs=10hz。列出系統(tǒng)函數(shù)并做頻率響應(yīng)曲線和零極點(diǎn)分布圖。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);[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)');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,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,'XTickMode','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è)計(jì)一個(gè)橢圓型數(shù)字帶阻濾波器，要求：阻帶0.4pi~0.6pi，20db，通帶0.35pi，0.65pi，1db，濾波器采樣周期為Ts=0.1s。列出系統(tǒng)函數(shù)并做頻率響應(yīng)曲線和零極點(diǎn)分布圖。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(Omgp,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