实验六-用双线性变换法设计IIR数字滤波器

实验六用双线性变换法设计IIR数字滤波器一、实验目的学会运用MATLAB设计数字低通、带通、高通、带阻滤波器的设计方法。二、实验涉及的matlab子函数bilinear功能：双线性变换——将s域映射到z域。调用格式：[numd,dend]=bilinear(num,den,Fs)，将模拟域系统函数转换为数字域的系统函数，Fs为采样频率。三、实验原理下面举例说明用双线性变换法设计各种数字滤波器的过程。采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：wp=0.25*pi，rp=1db，ws=0.4*pi，as=15db，滤波器采样频率Fs=100hz。MATLAB源程序为：%数字滤波器指标wpd=0.25*pi;%滤波器的通带截止频率wsd=0.4*pi;%滤波器的阻带截止频率Rp=1;As=15;%输入滤波器的通阻带衰减指标%转换为模拟原型滤波器指标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')%计算阶数n和截止频率[z0,p0,k0]=buttap(n);%归一化切比雪夫1型原型设计ba=k0*poly(z0);%求原型滤波器系数baa=poly(p0);%求原型滤波器系数a[ba1,aa1]=lp2lp(ba,aa,wc);%变换为模拟低通滤波器%用双线性变换法计算数字滤波器系数[bd,ad]=bilinear(ba1,aa1,Fs)%双线性变换%求数字系统的频率特性[H,w]=freqz(bd,ad);dbH=20*log10(abs(H)/max(abs(H)));%化为分贝值subplot(2,2,1),plot(w,abs(H));ylabel('|H|');title('幅度响应');axis([0,pi,0,1.1]);gridsubplot(2,2,2),plot(w,angle(H));ylabel('\phi');title('相位响应');axis([0,pi,-4,4]);gridsubplot(2,2,3),plot(w,dbH);title('幅度响应(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('零极图');运行结果为：n=5wc=103.2023bd=0.00720.03620.07250.07250.03620.0072ad=1.0000-1.94341.9680-1.07020.3166-0.0392那么所求滤波器的系统函数为例2、采用双线性变换法设计一个椭圆数字高通滤波器，要求通带250hz，1db，阻带150hz，20db，滤波器采样频率为Fs=1000hz。MATLAB源程序为：%数字滤波器指标fs=150;fp=250;Fs=1000;T=1/Fs;wpd=fp/Fs*2*pi;%数字滤波器的通带截止频率wsd=fs/Fs*2*pi;%数字滤波器的阻带截止频率Rp=1;As=20;%输入滤波器的通阻带衰减指标%转换为模拟滤波器指标wp=(2/T)*tan(wpd/2);ws=(2/T)*tan(wsd/2);%模拟原型滤波器计算[n,wc]=ellipord(wp,ws,Rp,As,'s')%计算阶数n和截止频率[z0,p0,k0]=ellipap(n,Rp,As);%归一化椭圆原型设计ba=k0*poly(z0);%求原型滤波器系数baa=poly(p0);%求原型滤波器系数a[ba1,aa1]=lp2hp(ba,aa,wc);%变换为模拟高通滤波器%用双线性变换法计算数字滤波器系数[bd,ad]=bilinear(ba1,aa1,Fs)%双线性变换%求数字系统的频率特性[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('幅度响应');axis([0,Fs/2,0,1.1]);gridsubplot(2,2,2),plot(w/2/pi*Fs,angle(H)/pi*180);ylabel('\phi');title('相位响应');axis([0,Fs/2,-180,180]);gridsubplot(2,2,3),plot(w/2/pi*Fs,dbH);title('幅度响应(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('零极图');运行结果为n=3wc=2.0000e+003bd=0.2545-0.43220.4322-0.2545ad=1.00000.18900.71970.1574例3、采用双线性变换法设计一个切比雪夫1型数字带通滤波器，要求：通带0.3pi~0.7pi，1db，阻带0.2pi，0.8pi，20db，滤波器采样周期为Ts=0.1s。MATLAB源程序为：%双线性变换法设计数字带通%数字滤波器指标wpd1=0.3*pi;wpd2=0.7*pi;%数字滤波器的通带截止频率wsd1=0.2*pi;wsd2=0.8*pi;%数字滤波器的阻带截止频率Rp=1;As=20;%输入滤波器的通阻带衰减指标%转换为模拟滤波器指标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);%模拟通带带宽和中心频率%模拟原型滤波器计算[n,wn]=cheb1ord(wp,ws,Rp,As,'s')%计算阶数n和截止频率[z0,p0,k0]=cheb1ap(n,Rp);%设计归一化的模拟原型滤波器ba1=k0*poly(z0);%求原型滤波器系数baa1=poly(p0);%求原型滤波器系数a[ba,aa]=lp2bp(ba1,aa1,w0,bw);%变换为模拟带通滤波器%用双线性变换法计算数字滤波器系数[bd,ad]=bilinear(ba,aa,Fs)%求数字系统的频率特性[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('零极图');运行结果为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、采用双线性变换法设计一个切比雪夫1型数字带阻滤波器，要求：阻带0.3pi~0.7pi，20db，通带0.2pi，0.8pi，1db，滤波器采样周期为Ts=0.1s。MATLAB源程序为：%数字滤波器指标wsd1=0.3*pi;wsd2=0.7*pi;%数字滤波器的通带截止频率wpd1=0.2*pi;wpd2=0.8*pi;%数字滤波器的阻带截止频率Rp=1;As=20;%输入滤波器的通阻带衰减指标%转换为模拟滤波器指标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);%模拟通带带宽和中心频率%模拟原型滤波器计算[n,wn]=cheb1ord(wp,ws,Rp,As,'s')%计算阶数n和截止频率[z0,p0,k0]=cheb1ap(n,Rp);%设计归一化的模拟原型滤波器ba1=k0*poly(z0);%求原型滤波器系数baa1=poly(p0);%求原型滤波器系数a[ba,aa]=lp2bs(ba1,aa1,w0,bw);%变换为模拟带阻滤波器%用双线性变换法计算数字滤波器系数[bd,ad]=bilinear(ba,aa,Fs)%求数字系统的频率特性[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('零极图');运行结果为n=3wn=6.498561.5532bd=0.07360.00000.22080.00000.22080.00000.0736ad=1.00000.0000-0.9761-0.00000.85680.0000-0.2919四、实验任务1、采用双线性变换法设计一个巴特沃斯数字高通滤波器，要求通带0.35pi，1db，阻带0.2pi，15db，滤波器采样频率为Fs=10hz。列出系统函数并做频率响应曲线和零极点分布图。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('幅度响应');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('相位响应');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('幅度响应(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、采用双线性变换法设计一个椭圆型数字带阻滤波器，要求：阻带0.4pi~0.6pi，20db，通带0.35pi，0.65pi，1db，滤波器采样周期为Ts=0.1s。列出系统函数并做频率响应曲线和零极点分布图。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