1、2、IIR 滤波器的 MATLAB实现2.1 IIR滤波器的设计方法及原理IIR数字滤波器是一种离散时间系统，其系统函数为： )(1)(0zxyabZHNkkM假设MN，当MN时,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数 和 ，它是kab数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。2.1.1 用脉冲相应不变法设计IIR数字滤波器 利用模拟滤波器来设计数字滤波器，也就是使数字滤波器能模仿模拟滤波器的特性
2、，这种模仿可以从不同的角度出发。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满足 )()nTha式中,T是采样周期。如果令Ha(s)是ha(t)的拉普拉斯变换，H(z)为h(n)的Z变换，利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得(1-1)则可看出，脉冲响应不变法将模拟滤波器的 S 平面变换成数字滤波器的 Z 平面，这个从 s 到 z 的变换 z=esT是从 S 平面变换到 Z 平面的标准变换关系式。 kTjsXTjksXTzXkaskaes 21
3、)(1)(图 1-1 脉冲响应不变法的映射关系 由（1-1）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为 (1-2)这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内时，即(1-3)才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即| |0 时，| z|1。也就是说，S 平面的左半平面映射到 Z 平面的单位圆内，S 平面的右半平面映射到 Z 平面的单位圆外，S 平面的虚轴映射到 Z 平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是
4、稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为 S 平面与 Z 平面是单值的一一对应关系。S 平面整个 j轴单值地对应于 Z 平面单位圆一周，即频率轴是单值变换关系。这个关系如式（1-8）所示，重写如下：上式表明，S 平面上 与 Z 平面的 成非线性的正切关系，如图 7-7 所示。由图 7-7 看出，在零频率附近，模拟角频率 与数字频率 之间的变换关系接近于线性关系；但当 进一步增加时， 增长得越来越慢，最后当时， 终止在折叠频率 = 处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现
5、象。图 1-4 双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（1-8）及图 1-4 所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字2tanT- o2tanT滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性） ，不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图1-5 所示。图 1-5 双线性变换法幅度和相位
6、特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。2.2 Butterworth低通滤波器设计推导滤波器最小阶数与设计指标的关系：滤波器阶数就是其系统函数的极点个数。为了避免滤波器的复杂程度与我们的要求不匹配，造成不必要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设计要求的最小阶数。以 Butterworth 滤波器为例，推导其阶数的数学模型。若给出滤波器的设计指标为：通带截止频率
7、，阻带截止频率 ，通带ps最大纹波 ，阻带最小纹波 。因为滤波器幅频特性为：()pRdB()sRdB（2-1）21NcHj其中 为 3dB 截止频率，N 为滤波器阶数。所以当 以及c p时，可得到：so oo)j(aH)(ej ooo )(eargjH)j(arg(2-2)2211()()psNpsNccHjHj、然后由 与 的关系式可得到：()pRdBs(2-3)221010lg()(l psp pRps ssjjHjj由上面的(2-2)式和(2-3)式可以联立求得：(2-4)210101()()ppssRpNRcNss通过上面的结果我们就可以求得滤波器阶数 N 为：(2-5)10lg()p
8、sRps所以滤波器的最小阶数就是大于上式所求得的值的最小整数。IIR 滤波器的流程框图:开始读入数字滤波器技术指标将指标转换成归一化模拟低通滤波器的指标设计归一化的模拟低通滤波器阶数 N 和 3db 截止频率模拟域频率变换，将 G(P)变换成模拟低通滤波器H(s)用双线性变换法将 H(s)转换成数字带通滤波器H(z)输入信号后显示相关结果2.2.1 冲激响应不变法设计 Butterworth低通滤波器冲激响应不变法是对模拟滤波器的单位冲激响应 h(t)等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的系统函数。设定设计指标：通带截止频率 =2000，阻带p截止频率 =3000，通带最
9、大纹波 Rp=3dB，阻带最小纹波sRs=18dB,采样频率 Fs=10000Hz。输入信号： 21)()in2)0.5cos(txtftfMatlab 程序如下：wp=2000*pi;ws=3000*pi;fs=10000;Rp=3;Rs=18;N,wn=buttord(wp,ws,Rp,Rs,s);z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp2lp(Bap,Aap,wn);bz,az=impinvar(b,a,fs);figure(1);h,w=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs(h);title
10、(Butterworth LPF 幅频特性);xlabel(频率(Hz) );ylabel(幅值(dB) );grid on;subplot(2,1,2),plot(w,angle(h);title(Butterworth LPF 相频特性);xlabel(频率(Hz) );ylabel(相位(degree) );grid on;figure(2);f1=1000;f2=4000;N1=100;dt=1/fs;n=0:N1-1;t=n*dt;x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);结束subplot(2,1,1),plot(t,x);title(输入信号),gr
11、id on;y=filter(bz,az,x);subplot(2,1,2),plot(t,y, r-);title(输出信号) ,grid on;得到的 Butterworth 低通滤波器幅频特性和相频特性曲线为：通过语句N,wn=buttord(wp,ws,Rp,Rs,s)可以得到 Matlab 估算出的滤波器最小阶数为 N=6。用式(2-5)出来的 N=5.1925，取大于它的最小整数得到滤波器的最小阶数为 6，两者一致，验证了该数学模型的可靠性。输入信号和经过滤波后的输出信号图像为：2.2.2 双线性变换法设计数字低通滤波器下面是基于 Butterworth 模拟原型滤波器，使用双线性
12、变换法设计数字低通滤波器的过程。设定与上面不同的设计指标为：通带截止频率 =0.2,p阻带截止频率 =0.3，通带最大纹波 Rp=1dB，阻带最大纹波 Rs=15dB，采样s频率 Fs=20000Hz。同时，为了检测滤波器的性能，我们设定一样的输入信号，其中 f1=1000Hz，f2=4000Hz，将该信号21)()in2)0.5co(txtftf与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。Matlab 程序如下：Wp=0.2*pi;Ws=0.3*pi;fs=20000;T=1/fs;Rp=1;Rs=15;wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T;N,W
13、n=buttord(wp,ws,Rp,Rs,s);z,p,k=buttap(N);Bap,Aap=zp2tf(z,p,k);b,a=lp2lp(Bap,Aap,Wn);bz,az=bilinear(b,a,fs);figure(1);h,w=freqz(bz,az,N,fs);subplot(2,1,1),plot(w,abs(h);title(Butterworth LPF 幅频特性);xlabel(频率(Hz) );ylabel(幅值(dB) );grid on;subplot(2,1,2),plot(w,angle(h);title(Butterworth LPF 相频特性);xlabe
14、l(频率(Hz) );ylabel(相位(degree) );grid on;figure(2);f1=1000;f2=4000;N1=100;dt=1/fs;n=0:N1-1;t=n*dt;x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x);title(输入信号),grid on;y=filter(bz,az,x);subplot(2,1,2),plot(t,y, r-);title(输出信号) ,grid on;得到的 Butterworth 低通滤波器幅频特性和相频特性曲线为：由语句N,Wn=buttord(wp,ws
15、,Rp,Rs,s)我们可以用 Matlab得出滤波器的最小阶数 N=6，与上面一个滤波器一样。再利用上面我们推导出的式（2-5） ，将滤波器各项指标代入可以算出N=4.7433，取大于它的最小整数得到滤波器的最小阶数为 5。二者并不一致，但差距只有 1，说明我们推导出的计算Butterworth 滤波器最小阶数的数学模型存在误差，但误差不大。又因为 Matlab 计算滤波器的最小阶数也是估算出来的，所以这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。输入信号和经过滤波后的输出信号图像为：通过比较两组幅频特性和相频特性曲线我们可以发现，其实双线性变换法和冲激响应不变法在设计较为简单的 Butterworth 低通滤波器时差别并不大，并且都拥有不错的设计效果。
