IIR数字滤波器设计.doc

第一章 绪论IIR滤波器的设计可以利用模拟滤波器许多现成的设计公式、数据和表格，计算的工作量较小，所以IIR滤波器的设计设计方法主要有经典设计法、直接设计法和最大平滑滤波器设计方法。本文中主要介绍经典方法设计IIR数字低通滤波器。经典方法首先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换成满足给定技术指标的数字滤波器。在对滤波器进行实际设计时，整个过程的运算量是很大的。设计阶数较高的IIR滤波器时，计算量更大，设计过程中改变参数或滤波器类型时都要重新计算。利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效地设计数字滤波器，大大地简化计算量。MATLAB的信号处理工具箱是专门应用于信号处理领域的专用工具箱，它的两个基本组成就是滤波器的设计与实现部分以及谱分析部分。工具箱提供了丰富而简便的设计，使原来繁琐的程序设计简化成函数的调用。只要以正确的指标参数调用相应的滤波器设计程序或工具箱函数，便可以得到正确的设计结果，使用非常方便。第二章 IIR数字滤波器的设计2.3 数字滤波器的技术要求 滤波器的指标常常在频域给出。数字滤波器的频响特性函数一般为复函数，所以通常表示为 式（2-3）其中，|称为幅频特性函数;称为相频特性函数。幅频特性表示信号通过该滤波器后各频率成分的衰减情况，而相频特性反映各频率通过滤波器后在时间上的延时情况。一般对IIR数字滤波器，通常只用幅频响应函数来描述设计指标，相频特性一般不作要求。而对线性相位特性的滤波器，一般用FIR数字滤波器设计实现。IIR低通滤波器的幅度特性如图2-3所示 图2-3 低通滤波器的幅度特性图中，分别称为通带截止频率和阻带截止频率。通带频率范围为，在通带中，要求（1-）0)；为模拟滤波器，缺省时为数字滤波器；为滤波器类型，为高通滤波器，截止频率；为带阻滤波器，；缺省为低通或带通滤波器。低通、高通滤波器时，截止频率；带通或带阻滤波器时，。a，b分别为滤波器的传递函数分子和分母式项式系数向量；z，p，k分别为滤波器的零极点和增益；A，B，C，D为滤波器状态空间表达式矩阵。滤波器传递函数具有下面形式 式（3-1）状态空间表达式 式（3-2） 式（3-3）若滤波器为带通或带阻型，则滤波器的阶数为2n，否则阶数为n。(2) 切比雪夫型模拟滤波器函数cheby1用于Chebyshev滤波器设计，调用格式为：b,a= cheby1 (n,Rp,)b,a= cheby1 (n,Rp,)z,p,k = cheby1 (n,Rp,)A,B,C,D = cheby1 (n,Rp,)其中，Rp为通带波纹（dB），其余各项意义与函数butter相同。(3) 切比雪夫型模拟滤波器函数cheby2用于Chebyshev滤波器设计，调用格式为：b,a= cheby2 (n,Rs,)b,a= cheby2(n,Rs,)z,p,k = cheby2(n,Rs,)A,B,C,D = cheby2 (n,Rs,)其中，Rs为阻带衰减（dB），其余各项意义与函数butter相同。(4) 椭圆模拟滤波器函数ellip用于椭圆滤波器设计，调用格式为：b,a= ellip (n,Rp, Rs ,)b,a= ellip (n,Rp, Rs ,)z,p,k = ellip (n,Rp, Rs ,)A,B,C,D = ellip (n,Rp, Rs ,)其中，Rp为通带波纹（dB），Rs为阻带衰减（dB），其余各项意义与函数butter相同。(5) 贝塞尔模拟滤波器函数besself用于贝塞尔(Bessel)滤波器设计，调用格式为：b,a= besself (n, )b,a= cheby1 (n, , )z,p,k = cheby1 (n,)A,B,C,D = cheby1 (n, ,)其中，n为滤波器阶数，为滤波器截止频率（-3dB频率），其余各项意义与函数butter相同。3.2.2 最小阶数选择函数Matlab信号处理工具箱提供各种滤波器最小阶数选择工具函数，这些函数和滤波器设计函数连接起来可方便成功地设计各类滤波器。（1）Butterworth模拟滤波器阶数选择函数Buttord，调用格式为 （2）Chebyshev模拟滤波器阶数选择函数Cheb1ord，调用格式为（3）Chebyshev模拟滤波器阶数选择函数Cheb2ord，调用格式为（4）椭圆模拟滤波器阶数选择函数Ellipord，调用格式为 其中，为通带边界频率，单位rad/s；为阻带边界频率，单位rad/s；为通带波动，dB,此值是滤波器在通带0之间允许最大允许幅值损失；为阻带衰减，dB，是滤波器幅值从通带至阻带下降的分贝数；表示模拟滤波器（缺省时该函数用于数字滤波器）；函数返回值n为模拟滤波器最小阶数；为模拟滤波器的截止频率（-3dB），rad/s。这四个函数适用于高通、带通、带阻滤波器。3.2.3 模拟原型滤波器设计函数模拟低通滤波器原型函数有buttap、cheb1ap、cheb2ap、ellipap、besselap。（1） 巴特沃斯模拟低通滤波器原型函数buttap，函数调用格式z,p,k=buttap(n)（2） 切比雪夫模拟低通滤波器原型函数cheb1ap，函数调用格式z,p,k=cheb1ap(n,Rp)（3） 切比雪夫模拟低通滤波器原型函数cheb2ap，函数调用格式z,p,k=cheb2ap(n,Rs)（4） 椭圆模拟低通滤波器原型函数ellipap，函数调用格式 z,p,k=ellipap(n,Rp,Rs)（5） 贝塞尔模拟低通滤波器原型函数besselap，函数调用格式z,p,k=besselap(n)3.2.4 频率转换函数Matlab信号处理工具箱提供的频率转换函数为lp2lp、lp2hp、lp2bp、lp2bs。（1）lp2lp用于实现低通模拟滤波器原型滤波器至低通滤波器的频率转换。调用格式为bt,at=lp2lp(b,a,)At,Bt,Ct,Dt=lp2lp(A,B,C,D,)式中，a,b或A，B，C，D为模拟原型滤波器模型，返回低通滤波器的模型，用传递函数形式（bt,at）或状态空间形式（At,Bt,Ct,Dt）（2）lp2hp用于实现低通模拟滤波器原型滤波器至高通滤波器的频率转换。调用格式为bt,at=lp2hp(b,a,)At,Bt,Ct,Dt=lp2hp(A,B,C,D,)（3）lp2bp用于实现低通模拟滤波器原型滤波器至带通滤波器的频率转换。调用格式为bt,at=lp2bp(b,a,Bw)At,Bt,Ct,Dt=lp2bp(A,B,C,D,Bw)其中，Bw为带通滤波器带宽，rad/s。（4）lp2bs用于实现低通模拟滤波器原型滤波器至带阻滤波器的频率转换。调用格式为bt,at=lp2bs(b,a,Bw)At,Bt,Ct,Dt=lp2bs(A,B,C,D,Bw)其中，Bw为带阻滤波器带宽，rad/s。3.2.5滤波器离散化函数(1) impinvar 基于脉冲响应不变法实现模拟滤波器至数字滤波器的转换。调用格式：bz,az=impinvar(b,a,Fs)其中b，a为模拟滤波器分子、分母多项式系数；Fs为采样频率，缺省值为1Hz；bz，az为数字滤波器的分子、分母多项式系数。(2) bilinear 双线性变换法实现S域至数字Z域的映射，直接用于模拟滤波器变换为数字滤波器。对于不同形式的模拟滤波器模型，函数调用格式不同。分为零极点增益形式、传递函数形式、状态空间形式。格式分别如下： 零极点增益形式：zd,pd,kd=bilinear(z,p,k,Fs) 传递函数形式: numd,dend=bilinear(num,den, Fs)状态空间形式: Ad,Bd,Cd,Dd=(A,B,C,D,Fs)式中，z，p分别为模拟滤波器零、极点列向量；k为模拟滤波器增益；Fs为采样频率，单位为Hz；zd，pd，kd 为数字滤波器零、极点，增益。3.3 MATLAB中 IIR低通数字滤波器的设计方法及设计步骤MATLAB中利用原型模拟滤波器设计IIR低通数字滤波器的步骤为（1） 根据给定的性能指标和方法首先对设计性能指标中的频率指标进行转换，转换后的频率指标作为模拟滤波器原型设计性能指标。（2） 估计模拟滤波器最小阶数和边界频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。（3） 设计模拟低通滤波器原型，利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。（4） 由模拟低通原型经频率变换获得数字滤波器（低通）利用MATLAB函数lp2lp。（5） 将模拟滤波器离散化获得IIR低通数字滤波器，利用MATLAB工具函数blinear、impinvar。确定模拟滤波器的设计性能指标模拟滤波器的最小阶数、边界频率开始设计模拟低通滤波器原型频率变换，获得模拟低通滤波器离散化，获得IIR数字滤波器绘图结束设计流程图为图3-1 Matlab中设计IIR滤波器的流程图3.4 IIR低通数字滤波器的设计实例3.4.1 用脉冲响应不变法设计IIR低通数字滤波器的实例实例1用脉冲响应不变法设计一个BUTTERWORTH低通数字滤波器，使其特性逼近一个BUTTERWORTH低通模拟滤波器的性能指标如下：通带截止频率rad/s，阻带截止频率rad/s，通带波纹Rp小于3dB，阻带衰减Rs大于15dB，采样频率。设计程序wp=2000*2*pi; %技术指标ws=5000*2*pi;Rp=3;Rs=15;Nn=128;Fs=10000;N,Wn=buttord(wp,ws,Rp,Rs,s);%计算模拟滤波器的阶数和-3dB截止频率z,p,k=buttap(N);%设计模拟低通巴特沃斯滤波器原型Bp,Ap=zp2tf(z,p,k);%零极点形式转换成传输函数形式b,a=lp2lp(Bp,Ap,Wn);%把模拟低通滤波器原型转换成低通滤波器bz,az=impinvar(b,a,Fs);%离散化，得到数字滤波器h,w=freqz(bz,az,Nn,Fs);%数字滤小器的频率响应subplot(2,2,2); %绘制零极点图zplane(bz,az);title(系统零极点图);Hmag=abs(h);%系统频率响应的模Hpah=angle(h);%系统频率响应的相位RHpah=unwrap(Hpah);%展开freqz产生的频率wsubplot(2,2,1); %绘制系统幅频图plot(w/(2*pi),20*log10(Hmag);xlabel(频率(Hz);ylabel(|H(jw)|2(dB);title(平方幅频曲线);grid on;subplot(2,2,3);%绘制系统相频图plot(w/(2*pi),RHpah);xlabel(频率(Hz);ylabel(相位(degrees);title(相频曲线);grid on;仿真结果如图3-2所示图3-2 脉冲响应不变法Butterworth低通滤波器的幅频、相频图实例2用脉冲响应不变法设计一个切比雪夫1型低通数字滤波器，使其特性逼近一个切比雪夫1低通模拟滤波器的性能指标如下：通带截止频率rad/s，阻带截止频率rad/s，通带波纹Rp小于3dB，阻带衰减大于15dB，采样频率。程序设计wp=3000*2*pi; %技术指标ws=5000*2*pi;Rp=3;Rs=30;Nn=512;Fs=10000;N,Wn=cheb1ord(wp,ws,Rp,Rs,s) %计算模拟切比雪夫型滤波器的阶数和-3dB截止频率z,p,k=cheb1ap(N,Rp); %设计模拟低通切比雪夫型滤波器原型Bp,Ap=zp2tf(z,p,k); %零极点形式转换成传输函数形式b,a=lp2lp(Bp,Ap,Wn); %把模拟低通滤波器原型转换成低通滤波器bz,az=impinvar(b,a,Fs); %离散化，得到数字滤波器h,w=freqz(bz,az,Nn,Fs); %数字滤小器的频率响应subplot(2,2,2); %绘制零极点图zplane(bz,az);title(系统零极点图);Hmag=abs(h);%系统频率响应的模Hpah=angle(h);%系统频率响应的相位RHpah=unwrap(Hpah); %展开freqz产生的频率wsubplot(2,2,1); %绘制系统幅频图plot(w/(2*pi),20*log10(Hmag);xlabel(频率(Hz);ylabel(|H(jw)|2(dB);title(平方幅频曲线);grid on;subplot(2,2,3);%绘制系统相频图plot(w/(2*pi),RHpah);xlabel(频率(Hz);ylabel(相位(degrees);title(相频曲线);grid on;仿真结果图3-3 脉冲响应不变法切比雪夫1型低通滤波器的幅频、相频图实例3用脉冲响应不变法设计一个10阶数字低通椭圆滤波器，其技术指标为:通带最大衰减3db， 阻带最小衰减60db，截止频率6pi弧度， T为0.1。程序设计n=10;%技术指标Rp=3;Rs=60;w=6*pi;t=0.1;fs=1/t;Nn=512;z,p,k=ellipap(n,Rp,Rs); %设计模拟椭圆滤波器原型A,B,C,D=zp2ss(z,p,k); %表达形式从零极点增益形式转换成状态方程形式At,Bt,Ct,Dt=lp2lp(A,B,C,D,w); %把模拟低通滤波器原型转换成模拟低通滤波器b,a=ss2tf(At,Bt,Ct,Dt); %表达形式从状态方程形式转换成传递函数形式b1,a1=impinvar(b,a,fs); %用脉冲响应不变法将模拟滤波器转换成数字滤波器h,W=freqz(b1,a1,Nn,fs);%数字滤波器频率响应subplot(2,2,2); %绘制零极点图zplane(bz,az);title(系统零极点图);Hmag=abs(h);%系统频率响应的模Hpah=angle(h);%系统频率响应的相位RHpah=unwrap(Hpah); %展开freqz产生的频率wsubplot(2,2,1); %绘制系统幅频图plot(w/(2*pi),20*log10(Hmag);xlabel(频率(Hz);ylabel(|H(jw)|2(dB);title(平方幅频曲线);grid on;subplot(2,2,3);%绘制系统相频图plot(w/(2*pi),RHpah);xlabel(频率(Hz);ylabel(相位(degrees);title(相频曲线);grid on;仿真结果图3-4 脉冲响应不变法椭圆低通滤波器的幅频图3.4.2 用双线性变换法设计IIR低通数字滤波器的实例实例1用双线性变换法设计一个BUTTERWORTH低通数字滤波器，使其特性逼近一个BUTTERWORTH低通模拟滤波器的性能指标如下：通带截止频率rad/s，阻带截止频率rad/s，通带波纹Rp小于3dB，阻带衰减大于15dB，采样频率。程序设计wp=2000*2*pi;%技术指标ws=3000*2*pi;Rp=3;Rs=15;Nn=512;Fs=10000;N,Wn=buttord(wp,ws,Rp,Rs,s)z,p,k=buttap(N);Bp,Ap=zp2tf(z,p,k);b,a=lp2lp(Bp,Ap,Wn);bz,az=bilinear(b,a,Fs);%用双线性变换法将模拟滤波器转换成数字滤波器h,w=freqz(bz,az,Nn,Fs);%数字滤波器频率响应subplot(2,2,2); %绘制零极点图zplane(bz,az);%pzmap(b,a);title(系统零极点图);Hmag=abs(h);%系统频率响应的模Hpah=angle(h);%系统频率响应的相位RHpah=unwrap(Hpah);% 展开freqz产生的频率wsubplot(2,2,1); %系统幅频图plot(w/(2*pi),20*log10(Hmag);xlabel(频率(Hz);ylabel(|H(jw)|2(dB);title(平方幅频曲线);grid on;subplot(2,2,3);%绘制系统的相频图plot(w/(2*pi),RHpah);xlabel(频率(Hz);yla