标准文案实验5脉冲响应不变法设计IIR数字滤波器一、实验目的1.掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。二、实验设备及环境计算机、matlab软件环境。三、实验基础理论1.基本原理从时域响应出发，是数字录波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲击响应ha(t)，h(n)等于ha(t)的取样值。2变换方法Ha(s)拉氏反变换→ha(t)时域采样→?a(nT)=h(n)z变换→H(z)将Ha(s)进行部分分式展开Ha(s)=∑Aks-pkNk=1(2)对Ha(s)进行拉式变换ha(t)=∑Akepku(t)Nk=1标准文案(3)对ha(t)时域采样得到h(n)h(n)=∑AkepknTu(nT)Nk=1=∑Ak(epkT)nu(n)Nk=1(4)对h(n)进行z变换H(z)=∑Ak1-epkz-1Nk=13.设计步骤(1)确定数字滤波器性能指标ωp，ωst，Rp和As。(2)将数字滤波器频率指标转换成相应的模拟滤波器频率指标Ωp=ωpTΩst=ωstT(3)根据指标ωp，ωst，Rp和As设计模拟滤波器Ha(s)。(4)将Ha(s)展成部分分式形式Ha(s)=∑Aks-pkNk=1(5)将模拟极点Pk转换成数字极点ePkT，得到数字滤波器Ha(s)=∑Ak1-epkTz-1Nk=1可见Ha(s)至H(z)间的变换关系为1s-sk<=>11-eskTz-1=zz-eskT在MATLAB中有两种方法可以实现上述变换。方法1：利用residue函数和residuez函数实现脉冲响应不变法，这两个函数的使用方法如下标准文案[r,p,k]=residue(b,a)[b,a]=reisdue(r,p,k)实现多项式形式H(s)=bMsM+bM-1sM-1+?+b0aNsN+aN-1sN-1+?+a0和部分分式形式H(s)=r1s-p1+r2s-p2+?+rNs-pN+k(s)之间的转换。[r,p,k]=residuez(b,a)[b,a]=reisduez(r,p,k)实现多项式形式H(z)=b0+b1z-1+?+bMz-Ma0+a1z-1+?+aNz-N和部分分式形式H(z)=r11-p1z-1+r21-p2z-1+?+rk1-pkz-1+k1+k2z-1+?之间的转换。方法2：MATLAB提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其使用方法如下[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系统向量bz和az，fs为采样频率。标准文案[bz,az]=impinvar(b,a)采样频率默认为1的情况下，采用脉冲响应不变法将模拟滤波器变换为数字滤波器。四、实验内容1.设采用率为fs=4kHz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=1kHz。实验代码及结果：>>[b,a]=butter(3,0.5*pi,'s')b=0003.8758a=1.00003.14164.93483.8758>>[bz,az]=impinvar(b,a)bz=00.58130.21140az=1.0000-0.39840.2475-0.0432>>w=[0:500]*pi/500;>>[H,w]=freqz(bz,az);>>subplot(211);>>plot(w/pi,abs(H));标准文案>>xlabel('\omega(\pi)');>>ylabel('|H(e^j^\omega)|');>>subplot(212);>>plot(w/pi,angle(H)/pi);>>xlabel('\omega(\pi)');>>ylabel('PhaseofH(e^j^\omega)(\pi)');2.设采用频率为fs=10kHz，设计数字低通滤波器，满足如下指标：通带截止频率：fp=1kHz，通带波动：Rp=1dB阻带截止频率：fst=1.5kHz,阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫1型、切比雪夫2型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。实验代码及结果：00.10.20.30.40.50.60.70.80.910.20.40.60.81()|H(ej)|00.10.20.30.40.50.60.70.80.91-1-0.500.51()PhaseofH(ej)()标准文案Butterworth滤波器：>>[N,Wn]=buttord(0.2*pi,0.3*pi,1,15,'s')N=6Wn=0.7087>>[b,a]=butter(N,Wn,'s')b=0000000.1266a=1.00002.73803.74843.25331.88240.69050.1266>>[bz,az]=impinvar(b,a)bz=-0.00000.00070.01050.01670.00420.00010az=1.0000-3.34435.0183-4.21902.0725-0.56000.0647标准文案>>w=[0:500]*pi/500;[H,w]=freqz(bz,az);subplot(211);plot(w/pi,abs(H));xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(212);plot(w/pi,angle(H)/pi);xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');(2)ChebyshevI型>>[N,Wn]=cheb1ord(0.2*pi,0.3*pi,1,15,'s')N=400.10.20.30.40....
