2、IIR滤波器的MATLAB实现

2.1 IIR滤波器的设计方法及原理

IIR数字滤波器是一种离散时间系统，其系统函数为：

H(Z)bzkk0Nk1Mk1akzky(z) x(z)假设M≤N，当M＞N时,系统函数可以看作一个IIR的子系统和一个(M-N)的FIR子系统的级联。IIR数字滤波器的设计实际上是求解滤波器的系数ak和bk ，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。

2.1.1 用脉冲相应不变法设计IIR数字滤波器

利用模拟滤波器来设计数字滤波器，也就是使数字滤波器能模仿模拟滤波器的特性，这种模仿可以从不同的角度出发。脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)，即将ha(t)进行等间隔采样，使h(n)正好等于ha(t)的采样值，满足 h(n)ha(nT)

式中,T是采样周期。如果令Ha(s)是ha(t)的拉普拉斯变换，H(z)为h(n)的Z变换，利用采样序列的Z变换与模拟信号的拉普拉斯变换的关系得

X(z)zesT112Xa(sjks)XasjkTkTkT(1-1)

则可看出，脉冲响应不变法将模拟滤波器的S平面变换成数字滤波器的Z平面，这个从s到z的变换z=esT是从S平面变换到Z平面的标准变换关系式。

oj3/ TjIm[z]/ T－/ T－3/ T－1o1Re[z]S平面Z平面 图1-1脉冲响应不变法的映射关系

由（1-1）式，数字滤波器的频率响应和模拟滤波器的频率响应间的关系为

12kH(e)HajTkTj(1-2)

这就是说，数字滤波器的频率响应是模拟滤波器频率响应的周期延拓。正如采样定理所讨论的，只有当模拟滤波器的频率响应是限带的，且带限于折叠频率以内时，即

Ha(j)0||Ts2(1-3)

才能使数字滤波器的频率响应在折叠频率以内重现模拟滤波器的频率响应，而不产生混叠失真，即

H(ej)1Ha |jω|<π (1-4) TT但是，任何一个实际的模拟滤波器频率响应都不是严格限带的，变换后就会产生周期延拓分量的频谱交叠，即产生频率响应的混叠失真，如图7-4所示。这时数字滤波器的频响就不同于原模拟滤波器的频响，而带有一定的失真。当模拟滤波器的频率响应在折叠频率以上处衰减越大、越快时，变换后频率响应混叠失真就越小。这时，采用脉冲响应不变法设计的数字滤波器才能得到良好的效果。

图1-2脉冲响应不变法中的频响混叠现象

…－3－2－o23…Ha(jΩ)2πT－πoTπT2πTH(ej)＝T对某一模拟滤波器的单位冲激响应ha(t)进行采样，采样频率为fs，若使fs增加，即令采样时间间隔（T=1/fs）减小，则系统频率响应各周期延拓分量之间相距更远，因而可减小频率响应的混叠效应。 脉冲响应不变法优缺点：

从以上讨论可以看出，脉冲响应不变法使得数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，也就是时域逼近良好，而且模拟频率Ω和数字频率ω之间呈线性关系ω=ΩT。因而，一个线性相位的模拟滤波器（例如贝塞尔滤波器）通过脉冲响应不变法得到的仍然是一个线性相位的数字滤波器。

脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以，脉冲响应不变法只适用于限带的模拟滤波器(例如，衰减特性很好的低通或带通滤波器)，而且高频衰减越快，混叠效应越小。至于高通和带阻滤波器，由于它们在高频部分不衰减，因此将完全混淆在低频响应中。如果要对高通和带阻滤波器采用脉冲响应不变法，就必须先对高通和带阻滤波器加一保护滤波器，滤掉高于折叠频率以上的频率，然后再使用脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计复杂性和滤波器的阶数。

2.1.2 用双线性变换法设计IIR数字滤波器

脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Ｚ平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1-3所示。

S平面S1平面Z平面o/ T－11oRe[z]jj1jIm[z]o－/ T1图1-3双线性变换的映射关系

为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现

式中,T仍是采样间隔。

21TtanT2（1-5）

当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。将式（1-5）写成

2ej1T/2ej1T/2jj1T/2j1T/2Tee2es1T/2es1T/22s1T21es1Tss1T/2s1T/2tanhs1TTeeT2T1e将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得

再将S1平面通过以下标准变换关系映射到Z平面z=es1T 从而得到S平面和Z平面的单值映射关系为：

121z （1-6) s1T1z T21ss2T （1-7） zT21ss 2T式（1-6）与式（1-7）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换.

式（1-5）与式（1-6）的双线性变换符合映射变换应满足的两点要求。 首先,把z=ejω，可得

即S平面的虚轴映射到Z平面的单位圆。 其次，将s=σ+jΩ代入式（1-8），得 因此

21ej2sjtanjT1ejT2（1-8）

2jTz2jT22T|z|222T2

由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。 双线性变换法优缺点

双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单

2tanT2值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（1-8）所示，重写如下：

上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图7-7所示。

由图7-7看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。

－o＝2tan2T

图1-4双线性变换法的频率变换关系

但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（1-8）及图1-4所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正

是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图1-5所示。

Ha(jΩ)oooarg[Ha(jΩ)]ooH(ej)oarg[H(ej)] 图1-5双线性变换法幅度和相位特性的非线性映射

对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。

2.2 Butterworth低通滤波器设计

推导滤波器最小阶数与设计指标的关系：滤波器阶数就是其系统函数的极点个数。为了避免滤波器的复杂程度与我们的要求不匹配，造成不必要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设计要求的最小阶数。以Butterworth滤波器为例，推导其阶数的数学模型。若给出滤波器的设计指标为：通带截止频率p，阻带截止频率s，通带最大纹波Rp(dB)，阻带最小纹波Rs(dB)。因为滤波器幅频特性为：

H(j11()2N （2-1）

c其中c为3dB截止频率，N为滤波器阶数。所以当p以及s时，可得到：

H(jp11(p2N)c、H(js11(s)2N (2-2)

c然后由Rp(dB)与Rs(dB)的关系式可得到：

211R10lgH(j)ppRp2H(jp1010 (2-3) 211R10lgH(j)sRs2sH(js1010由上面的(2-2)式和(2-3)式可以联立求得：

Rpp2NRp1()1010pN10101 (2-4) c ()RsRss102Ns101011()10c通过上面的结果我们就可以求得滤波器阶数N为：

lg NRp1010Rs101011 (2-5)

plg()s所以滤波器的最小阶数就是大于上式所求得的值的最小整数。 IIR滤波器的流程框图:

开始 ↓ 读入数字滤波器技术指标 ↓ 将指标转换成归一化模拟低通滤波器的指标 ↓ 设计归一化的模拟低通滤波器阶数N和3db截止频率 ↓ 模拟域频率变换，将G(P)变换成模拟低通滤波器H(s) ↓ 用双线性变换法将H(s)转换成数字带通滤波器H(z) ↓ 输入信号后显示相关结果 ↓ 结束 2.2.1 冲激响应不变法设计Butterworth低通滤波器 冲激响应不变法是对模拟滤波器的单位冲激响应h(t)等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的系统函数。设定设计指标：通带截止频率p=2000π，阻带截止频率s=3000π，通带最大纹波Rp=3dB，阻带最小纹波Rs=18dB,采样频率Fs=10000Hz。 输入信号：

x(t)sin(2f1t)0.5cos(2f2t) Matlab程序如下：

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('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('输入信号'),grid 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模拟原型滤波器，使用双线性变换法设计数字低

通滤波器的过程。设定与上面不同的设计指标为：通带截止频率p=0.2π,

阻带截止频率s=0.3π，通带最大纹波Rp=1dB，阻带最大纹波Rs=15dB，采样频率Fs=20000Hz。同时，为了检测滤波器的性能，我们设定一样的输入信号

x(t)sin(2f1t)0.5cos(2f2t)，其中

Matlab程序如下：

f1=1000Hz，f2=4000Hz，将该信号

与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。

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,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]=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相频特性');

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('输入信号'),grid 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=4.7433，取大于它的最小整数得到滤波器的最小阶数为5。二者并不一致，但差距只有1，说明我们推导出的计算Butterworth滤波器最小阶数的数学模型存在误差，但误差不大。又因为Matlab计算滤波器的最小阶数也是估算出来的，所以这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。 输入信号和经过滤波后的输出信号图像为：

通过比较两组幅频特性和相频特性曲线我们可以发现，其实双线性变换法和冲激响应不变法在设计较为简单的Butterworth低通滤波器时差别并不大，并且都拥有不错的设计效果。