脉冲响应不变法在IIR滤波器设计中的应用(实例程序).doc

脉冲响应不变法在IIR滤波器设计中的应用（实例程序）摘自：张登奇,彭鑫. 脉冲响应不变法在IIR滤波器设计中的应用J. 湖南理工学院学报(自然科学版). 2015(04)摘 要：数字滤波器的设计是“数字信号处理”的重要教学内容，脉冲响应不变法是间接设计IIR数字滤波器的基本方法。文章以巴特沃斯数字低通滤波器的设计为例，根据脉冲响应不变法的基本原理，归纳了该法的设计步骤，并对设计过程进行了深入分析，最后实例列出了MATLAB实现的程序。3 MATLAB实现程序根据实例1确定的滤波器设计指标，用脉冲响应不变法设计巴特沃斯数字低通滤波器的MATLAB程序如下：%技术指标：fp=3000Hz,Ap=2dB，fs=6000Hz，As=7db,Fs=20000Hz。clc; clear; close all; format compact;wp=0.3*pi, Ap=2, ws=0.6*pi, As=7, Ts=1/20000;%输入数字指标T=1,%改变采样周期不影响H(z)的设计结果Wp=wp/T; Ws=ws/T;%将数字指标转化为过渡模拟指标M=log10(10(0.1*Ap)-1)/(10(0.1*As)-1) / (2*log10(Wp/Ws),%计算滤波器阶数N=ceil(M),%滤波器阶数向上取整Wcp=Wp*(10(0.1*Ap)-1)(-1/(2*N);%通带边界指标精确满足的3dB截止频率Wcs=Ws*(10(.1*As)-1)(-1/(2*N);%阻带边界指标精确满足的3dB截止频率Wc=Wcp,%模拟滤波器通带边界指标精确满足%Wc=Wcs,%模拟滤波器阻带边界指标精确满足%Wc=(Wcp+Wcs)/2,%设计的模拟滤波器会使通带阻带边界指标都有余量bp,ap=butter(N,1,s);%求归一化低通原型系统函数Ga(p)tf(bp,ap,variable,p),%显示Ga(p)bs,as=lp2lp(bp,ap,Wc);%去归一化得模拟滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)Ak,sk=residue(bs,as);%将Ha(s)按部分分式形式展开ak=T*Ak;zk=exp(sk*T);%将Ha(s)的部分分式参数转换为H(z)的部分分式参数bz,az=residuez(ak,zk,0);%将H(z)的部分分式形式化为分子分母形式tf(bz,az,Ts,variable,z-1),%显示系统函数freqz(bz,az,100),%绘频响曲线检验设计指标附图及程序：%用脉冲响应不变法设计butterworth数字低通滤波器 clc; clear; close all; format compact;%程序初始化%输入数字指标wp=0.3*pi; Ap=2; ws=0.6*pi; As=7;%数字频率及分贝数%利用归一化频率及倍数绘数字指标图wp1=wp/pi,AP=10(-Ap/20),ws1=ws/pi,AS=10(-As/20),%指标归一化频率及倍数T=1;%用于设计过渡模拟滤波器，改变采样周期不影响H(z)的设计结果subplot(2,6,1),stem(0,1.5,k.),hold on,stem(0,-0.5,k.),plot(-wp1,wp1,1,1,:),plot(-wp1,wp1,AP,AP,:),plot(-1,-ws1,AS,AS,:),plot(ws1,1,AS,AS,:),stem(-ws1,-wp1,wp1,ws1,AS,AP,AP,AS,r:.),axis(-1,1,-0.2,1.2)%绘Ha(s)的幅频曲线图Wp=wp/T; Ap; Ws=ws/T; As;%将数字指标转化为过渡模拟指标M=log10(10(0.1*Ap)-1)/(10(0.1*As)-1) / (2*log10(Wp/Ws),%计算滤波器阶数N=ceil(M),%滤波器阶数向上取整%N=2,Wcp=Wp*(10(0.1*Ap)-1)(-1/(2*N);%通带边界精确满足的截止频率Wcs=Ws*(10(.1*As)-1)(-1/(2*N);%阻带边界精确满足的截止频率Wc=Wcp,%截止频率用通带边界精确满足的截止频率%Wc=(Wcp+Wcs)/2,%通带阻带边界都有余量的截止频率%Wc=Wcs,%截止频率用阻带边界精确满足的截止频率bp,ap=butter(N,1,s);%求归一化原型滤波器系统函数Ga(p)P157tf(bp,ap,variable,p),%显示Ga(p)bs,as=lp2lp(bp,ap,Wc),%去归一化得过渡滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)，分子不足前面补0%利用符号运算和数值运算计算模拟滤波器幅频曲线syms s W t %定义符号对象 A_s=poly2sym(as,s);B_s=poly2sym(bs,s);Has=B_s/A_s;s=j*W;%计算系统函数，令s=j*wHjw=compose(Has,s); %利用复合函数计算并显示频率响应解析解W_num=linspace(-3,3,601);%在-3到3 范围内取300个频点Hjw_num=subs(Hjw,W,W_num*pi/T);%利用置换函数求频率响应各频点的数值解Ha_mag=abs(Hjw_num);subplot(2,6,2:4),stem(0,1.5,k.),hold on,stem(0,-0.5,k.),plot(-wp1,wp1,1,1,:),plot(-wp1,wp1,AP,AP,:),plot(-1,-ws1,AS,AS,:),plot(ws1,1,AS,AS,:),stem(-ws1,-wp1,wp1,ws1,AS,AP,AP,AS,r:.),plot(W_num,Ha_mag),axis(-3.6,3.6,-0.2,1.2),%绘Ha(s)的冲激响应ha(t)图hat=ilaplace(Has);ht1=vpa(hat,4);%求冲激响应,保留4位有效数字,输出便读式 t1=linspace(0,10,101);%在0到 秒取100个时点ht=subs(ht1,t,t1);%利用置换函数求冲激响应各时点的数值解subplot(2,6,5:6),stem(0,1.5,k.),hold on,stem(0,-0.5,k.),plot(t1,ht),axis(0,10,-0.1,0.6),%绘H(z)的冲激响应h(n)图n=0:10;hn=T*subs(ht1,t,n*T);%利用置换函数求冲激响应各时点的数值解subplot(2,6,11:12),plot(t1,ht,:),hold on,stem(n*T,hn,.),axis(0,10,-0.1,0.6),%绘Ha(s)的幅频曲线图subplot(2,6,8:10),stem(0,1.5,k.),hold on,stem(0,-0.5,k.),plot(-wp1,wp1,1,1,:),plot(-wp1,wp1,AP,AP,:),plot(-1,-ws1,AS,AS,:),plot(ws1,1,AS,AS,:),stem(-ws1,-wp1,wp1,ws1,AS,AP,AP,AS,r:.),plot(W_num,Ha_mag),plot(W_num-2,Ha_mag,:),plot(W_num+2,Ha_mag,:),axis(-3.6,3.6,-0.2,1.2),%绘数字指标及H(z)的幅频曲线图subplot(2,6,7),stem(0,1.5,k.),hold on,stem(0,-0.5,k.),plot(-wp1,wp1,1,1,:),plot(-wp1,wp1,AP,AP,:),plot(-1,-ws1,AS,AS,:),plot(ws1,1,AS,AS,:),stem(-ws1,-wp1,wp1,ws1,AS,AP,AP,AS,r:.),Ak,sk=residue(bs,as),%将Ha(s)按部分分式形式展开增益极点形式ak=T*Ak,zk=exp(sk*T),%将Ha(s)的增益极点转换为H(z)的增益极点bz,az=residuez(ak,zk,0),%将H(z)的部分分式形式化为分子分母等阶形式tf(bz,az,variable,z-1),%显示系统函数bz1,az1 = impinvar(bs,as,1/T)%可调用impinvar函数验证tf(bz1,az1,T,variable,z-1),%显示验证系统函数,z-1式的分子不足是后面补0w=linspace(-1,1,201); %确定频点z=exp(j*w*pi); %求频点对应的z点Hz=polyval(bz,z)./polyval(az,z); %求各频点的频响plot(w,abs(Hz) %绘制幅频曲线axis(-1,1,-0.2,1.2)T=1.3;%用于设计过渡模拟滤波器，改变采样周期不影响H(z)的设计结果Wp=wp/T; Ap; Ws=ws/T; As;%将数字指标转化为过渡模拟指标M=log10(10(0.1*Ap)-1)/(10(0.1*As)-1) / (2*log10(Wp/Ws),%计算滤波器阶数N=ceil(M),%滤波器阶数向上取整%N=2,Wcp=Wp*(10(0.1*Ap)-1)(-1/(2*N);%通带边界精确满足的截止频率Wcs=Ws*(10(.1*As)-1)(-1/(2*N);%阻带边界精确满足的截止频率Wc=Wcp,%截止频率用通带边界精确满足的截止频率%Wc=(Wcp+Wcs)/2,%通带阻带边界都有余量的截止频率%Wc=Wcs,%截止频率用阻带边界精确满足的截止频率bp,ap=butter(N,1,s);%求归一化原型滤波器系统函数Ga(p)P157tf(bp,ap,variable,p),%显示Ga(p)bs,as=lp2lp(bp,ap,Wc),%去归一化得过渡滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)，分子不足前面补0%利用符号运算和数值运算计算模拟滤波器幅频曲线syms s W t %定义符号对象 A_s=poly2sym(as,s);B_s=poly2sym(bs,s);Has=B_s/A_s;%绘Ha(s)的冲激响应ha(t)图hat=ilaplace(Has);ht1=vpa(hat,4);%求冲激响应,保留4位有效数字,输出便读式 t1=linspace(0,10,101);%在0到 秒取100个时点ht=subs(ht1,t,t1);%利用置换函数求冲激响应各时点的数值解subplot(2,6,5:6),hold on,plot(t1,ht,r:),%绘H(z)的冲激响应h(n)图n=0:10;hn=T*subs(ht1,t,n*T);%利用置换函数求冲激响应各时点的数值解subplot(2,6,11:12),hold on,plot(t1,T*ht,r:),stem(n*T,hn,r:.),%用脉冲响应不变法设计butterworth数字低通滤波器 clc; clear; close all; format compact;%程序初始化%输入数字指标wp=0.3*pi; Ap=2; ws=0.6*pi; As=7;%数字频率及分贝数%利用归一化频率及倍数绘数字指标图wp1=wp/pi,AP=10(-Ap/20),ws1=ws/pi,AS=10(-As/20),%指标归一化频率及倍数T=1;%用于设计过渡模拟滤波器，改变采样周期不影响H(z)的设计结果stem(wp1,AP,:.),hold on,stem(ws1,AS,r:.),%用倍率绘数字指标Wp=wp/T; Ap; Ws=ws/T; As;%将数字指标转化为过渡模拟指标M=log10(10(0.1*Ap)-1)/(10(0.1*As)-1) / (2*log10(Wp/Ws),%计算滤波器阶数N=ceil(M),%滤波器阶数向上取整Wcp=Wp*(10(0.1*Ap)-1)(-1/(2*N);%通带边界精确满足的截止频率Wcs=Ws*(10(.1*As)-1)(-1/(2*N);%阻带边界精确满足的截止频率WCP=Wcp/pi,WCS=Wcs/pi,%归一化3dB模拟截止频率Wc=Wcp,%截止频率用通带边界精确满足的截止频率%Wc=(Wcp+Wcs)/2,%通带阻带边界都有余量的截止频率%Wc=Wcs,%截止频率用阻带边界精确满足的截止频率bp,ap=butter(N,1,s);%求归一化原型滤波器系统函数Ga(p)P157tf(bp,ap,variable,p),%显示Ga(p)bs,as=lp2lp(bp,ap,Wc),%去归一化得过渡滤波器系统函数Ha(s)tf(bs,as),%显示Ha(s)，分子不足前面补0%利用符号运算和数值运算计算模拟滤波器幅频曲线syms s W t %定义符号对象 A_s=poly2sym(as,s);B_s=poly2sym(bs,s);Has=B_s/A_s;s=j*W;%计算系统函数，令s=j*wHjw=compose(Has,s); %利用复合函数计算并显示频率响应解析解W_num=linspace(0,3,301);%在-3到3 范围内取300个频点Hjw_num=subs(Hjw,W,W_num*pi/T);%利用置换函数求频率响应各频点的数值解Ha_mag=abs(Hjw_num);plo