基于脉冲响应不变法的dsp设计

实验五利用脉冲响应不变法设计IIR数字滤波器一、实验目的1.掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。二、实验设备与环境计算机、MATLAB软件环境。三、实验基础理论1.基本原理从时域响应出发，使数字滤波器的单位脉冲响应()hn模仿模拟滤波器的单位冲击响应(),()ahthn等于()aht的取样值。2.变换方法()()()()zaaahshthnThnH思路：拉式反变换时域采样变换（z)（1）将()aHs进行部分分式展开1()NkakkAHssP（2）对()aHs进行拉式反变换1()()kNptakkhtAeut（3）对()aht时域采样得到()hn11()()()()kkNNpnTpTkkkkhnAeunTAeun（4）对()hn进行Z变换11()1zkNkpTkAhze3.设计步骤（1）确定数字滤波器性能指标1.5stfkHz（2）将数字滤波器频率指标转换成相应的模拟滤波器频率指标ppTststT（3）根据指标,,,pstpsRA设计模拟滤波器()aHs将()aHs展成部分分式形式1()NkakkAHssp（4）把模拟极点kp转换成数字极点kpTe，得到数字滤波器11()1kNkpTkAHzez可见()aHs至()Hz间的变换关系为1111kksTsTkzssezze在MATLAB中有两种方法可以实现上述变换。方法1：利用residue函数和residuez函数实现脉冲响应不变法，这两个函数的使用方法如下：[r,p,k]=residue(b,a)[b,a]=residue(r,p,k)实现多项式形式110110()MMMMNNNNbsbsbHsasasa和部分分式形式1212()()NNrrrHsksspspsp之间的转换。[r,p,k]=residuez(b,a)[b,a]=residuez(r,p,k)实现多项式形式101101()MMNNbbzbzHzaazaz和部分分式形式1121211112()111NNrrrHzkkzpzpzpz之间的转换。方法2：MATLAB提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用方法如下：[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a转换成为数字滤波器系统函数的系统向量bz和az,fs为采样率。[bz,az]=impinvar(b,a)采样频率默认为1的情况向下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。四、实验内容1.设采样频率,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率。实验代码与实验结果：N=3;fs=4000;fc=1000;Wc=0.25*pi;[b,a]=butter(N,Wc,'s');[bz,az]=impinvar(b,a);w=[0:500]*pi/500;[H,w]=freqz(bz,az);subplot221;plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot222;plot(w/pi,20*log10(abs(H)));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot223;plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;subplot224;grd=grpdelay(bz,az,w);plot(w/pi,grd);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;00.5100.51()|H(ej)|(dB)00.51-40-30-20-100()|H(ej)|(dB)00.51-1-0.500.51()|H(ej)|00.51-2024()|H(ej)|2.设采样频率，设计数字低通滤波器，满足如下指标通带截止频率：，通带波动：阻带截止频率：，阻带衰减：要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。实验代码与实验结果：（1）巴特沃斯fs=10000;fp=1000;fst=1500;Wp=0.1*pi;Ws=0.15*pi;Rp=1;As=15;N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));Wc=Wp/((10^(Rp/10)-1)^(1/(2*N)));[b,a]=butter(N,Wc,'s');[bz,az]=impinvar(b,a);w=[0:500]*pi/500;[H,w]=freqz(bz,az);subplot221;plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot222;plot(w/pi,20*log10(abs(H)));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot223;plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;subplot224;grd=grpdelay(bz,az,w);plot(w/pi,grd);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;00.5100.51()|H(ej)|(dB)00.51-150-100-500()|H(ej)|(dB)00.51-1-0.500.51()|H(ej)|00.5105101520()|H(ej)|（2）切比雪夫I型fs=10000;fp=1000;fst=1500;Rp=1;As=15;Ap=Rp;epc=sqrt(10^(Ap/10)-1);Wp=2*pi*fp;Wc=Wp;Wst=2*pi*fst;N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wst/Wc));[ba]=cheby1(N,Rp,Wc,'s');[bzaz]=impinvar(b,a,fs);w=[0:500]*pi/500;[Hw]=freqz(bz,az);subplot221;plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot222;plot(w/pi,20*log10(abs(H)));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot223;plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;subplot224;grd=grpdelay(bz,az,w);plot(w/pi,grd);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;00.5100.511.5()|H(ej)|(dB)00.51-100-50050()|H(ej)|(dB)00.51-1-0.500.51()|H(ej)|00.51051015()|H(ej)|（3）切比雪夫II型fs=10000;fp=1000;fst=1500;Rp=1;As=15;Ap=Rp;epc=sqrt(10^(Ap/10)-1);Wp=2*pi*fp;Wst=2*pi*fst;Wc=Wst;N=ceil(acosh(sqrt(10^(0.1*As)-1)/epc)/acosh(Wst/Wp));[ba]=cheby2(N,As,Wc,'s');[bzaz]=impinvar(b,a,fs)w=[0:500]*pi/500;[Hw]=freqz(bz,az);subplot221;plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot222;plot(w/pi,20*log10(abs(H)));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot223;plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;subplot224;grd=grpdelay(bz,az,w);plot(w/pi,grd);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;00.5100.511.5()|H(ej)|(dB)00.51-10-505()|H(ej)|(dB)00.51-1-0.500.51()|H(ej)|00.51-505()|H(ej)|（4）椭圆模拟原型滤波器fs=10000;fp=1000;fst=1500;Rp=1;As=15;Ap=Rp;epc=sqrt(10^(Ap/10)-1);Wp=2*pi*fp;Wst=2*pi*fst;Wc=Wp;A=10^(As/20);k1=epc/(sqrt(A^2-1));k=Wp/Wst;N=ceil(ellipke(k)*ellipke(sqrt(1-k1^2))/ellipke(k1)/ellipke(sqrt(1-k^2)));[bzaz]=ellip(N,Rp,As,fp/fs*pi)w=[0:500]*pi/500;[Hw]=freqz(bz,az);subplot221;plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot222;plot(w/pi,20*log10(abs(H)));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');gridon;subplot223;plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;subplot224;grd=grpdelay(bz,az,w);plot(w/pi,grd);xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');gridon;00.5100.511.5()|H(ej)|(dB)00.51-100-50050()|H(ej)|(dB)00.51-1-0.500.51()|H(ej)|00.51