数字信号处理实验四 IIR滤波器设计.docx

1、数字信号处理实验四 IIR滤波器设计 实验四 IIR滤波器设计一、教学目的和任务1熟悉用双线性变换法设计IIR数字滤波器的原理和方法；2了解用脉冲响应不变法设计IIR数字滤波器的原理和方法；3掌握双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点；4掌握数字滤波器的计算机仿真方法。二、实验原理介绍IIR数字滤波器的系统函数为的有理分式：设计IIR滤波器的系统函数，就是要确定的阶数N及分子分母多项式的系数和，使其满足指定的频率特性。由于模拟滤波器的设计有许多简单而严谨的设计公式和大量的图表可以利用，因此IIR滤波器设计的方法之一是：先设计一个合适的模拟滤波器

2、，然后将模拟滤波器通过适当的变换转换成满足给定指标的数字滤波器。1、Butterworth模拟低通滤波器幅度平方函数： 其中，N为滤波器的阶数，为通带截止频率。2Chebyshev模拟低通滤波器 3、脉冲响应不变法原理用数字滤波器的单位脉冲响应序列h(n)逼近模拟滤波器的冲激响应，让h(n)正好等于的采样值，即：其中，T为采样间隔。如果以和H(z)分别表示的拉氏变换及h(n)的Z变换，则：4、双线性变换法原理双线性变换法是通过两次映射采用非线性频率压缩的方法，将整个频率轴上的频率范围压缩到 /T之间，再用转换到z平面上，从而使数字滤波器的频率响应与模拟滤波器的频率响应相似。5、设计IIR数字滤

3、波器的步骤1）确定数字滤波器的通带频率、阻带频率，通带最大衰减和阻带最小衰减。2）计算对应的模拟低通滤波器的频率。3）确定模拟低通滤波器的阶数N和3dB截止频率。4）模拟低通滤波器的系统函数H（s）。5）由H（s）经过反归一化、脉冲响应不变法和双线性变换法确定数字低通滤波器的系统函数H(z)。6）设计其它形式的滤波器时，由模拟低通到所需类型滤波器的频率域变换直接得到。6*、MATLAB中用于IIR数字滤波器设计的函数1）滤波器的特性分析1Freqz函数：求解数字滤波器的频率响应h,w=freqz(b,a,n)：返回数字滤波器的n点复频率响应，输入参数b和a分别是滤波器系数的分子和分母向量；输出

4、参数h是复频率响应，w是频率点。输入参数n默认是512。Freqz(b,a,)：没有输出参数，直接在当前窗口中绘制频率响应的幅频响应和相频响应。2Freqs函数：求解模拟滤波器的频率响应h,w=freqz(b,a,n)：返回模拟滤波器的n点复频率响应，输入参数b和a分别是滤波器系数的分子和分母向量；输出参数h是复频率响应，w是频率点。输入参数n默认是512。3Abs、angle函数：分别用于从复频域响应数据中提取幅值信息和相位信息4Zplane函数：绘制系统的零极点图zplane(z,p)：以单位圆为基准；z为系统的零点向量，图中用o表示；p为系统的极点向量，图中用x表示。zplane(b,a

5、)：输入参数为系统传递函数的分子向量和分母向量。2）确定滤波器最小阶数函数函数功能n,wn=Buttord(wp,ws,rp,rs)估计Butterworth滤波器阶数n,wn=Cheb1ord(wp,ws,rp,rs)估计Chebyshev型滤波器阶数n,wn=Cheb2ord(wp,ws,rp,rs)估计Chebyshev型滤波器阶数n,wn=Ellipord(wp,ws,rp,rs)估计椭圆滤波器阶数 wp：归一化的通带截止频率； ws：归一化的阻带截止频率 rp：通带最大衰减量； rs：阻带最小衰减量 n：返回符合要求的滤波器阶数； wn：返回滤波器的截止频率3）模拟低通滤波器的设计函

6、数函数功能z,p,k=Buttap(n)返回Butterworth滤波器的零点、极点、增益z,p,k=Cheb1ap(n,rp)返回Chebyshev型滤波器的零点、极点、增益z,p,k=Cheb2ord(n,rs)返回Chebyshev型滤波器的零点、极点、增益z,p,k=Ellipap(n,rp,rs)返回椭圆滤波器的零点、极点、增益b,a=Butter(n,wn,s)返回Butterworth滤波器的分子分母多项式的系数b,a=Cheby1(n,rp,wn,s)返回Chebyshev型滤波器的分子分母多项式的系数b,a=Cheby2(n,rp,wn,s)返回Chebyshev型滤波器的分

7、子分母多项式的系数b,a=Ellip(n,rp,rs,wn,s)返回椭圆滤波器的分子分母多项式的系数4）模拟滤波器的离散化1Impinvar函数：模拟滤波器变换成数字滤波器的脉冲响应不变法bz,az=impinvar(b,a,fs)：将模拟滤波器(b,a)变换成数字滤波器(bz,az)；输入参数fs是对模拟滤波器频率响应的采样，默认为1。2bilinear函数：模拟滤波器转换为数字滤波器的双线性变换法zd,pd,kd=bilinear(z,p,k,fs)：将采样零极点模型表达的模拟滤波器转换为数字滤波器。列向量zd为零点向量，列向量pd为极点向量，kd为系统增益，fs是指定的采样频率。numd

8、,dend=bilinear(num,den,fs)：将采用传递函数模型表达的模拟滤波器转换为数字滤波器。ad,bd,cd,dd=bilinear(a,b,c,d,fs)：将采用状态空间模型表达的模拟滤波器转换为数字滤波器。5）直接设计IIR数字滤波器函数函数功能b,a=Butter(n,wn)设计巴特沃斯滤波器b,a=Cheby1(n,r,wn)设计切比雪夫型滤波器，r指定通带内波纹大小b,a=Cheby2(n,r,wn)设计切比雪夫型滤波器，r指定通带内波纹大小b,a=Ellip(n,rp,rs,wn)设计椭圆滤波器，rp指定通带内波纹最大衰减，rs指定通带内波纹的最小衰减wn：数字滤波器

9、的截止频率；当wn是一个二元向量w1,w2时，返回一个2n阶的带通滤波器，通带为w1 w2 。返回滤波器的分子分母多项式的系数，b为分子系数向量、a为分母系数向量（降幂排列）。三、实验相关知识准备1模拟滤波器的设计函数为了利用模拟滤波器设计数字滤波器，先必须设计对应的模拟滤波器，常用的模拟滤波器有：bessel滤波器，butterworth滤波器，chebyshev型滤波器，chebyshev型滤波器及椭圆函数滤波器。1模拟滤波器的设计函数1)设计bessel模拟低通滤波器z,p,k=besselap(n)2）设计butterworth模拟低通滤波器z,p,k=buttap(n)3）设计che

10、byshev型模拟低通滤波器z,p,k=cheb1ap(n,Rp)Rp：通带内的波纹系数，单位分贝4）设计chebyshev型模拟低通滤波器z,p,k=cheb2ap(n,Rs)Rs：阻带内的波纹系数低于通带Rs分贝5）设计椭圆模拟滤波器z,p,k=ellipap(n,Rp.Rs)2滤波器阶数的选择下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。1）选择butterworth滤波器阶数数字域：n,Wn=buttord(Wp,Ws,Rp,Rs)模拟域：n,Wn=buttord(Wp,Ws,Rp,Rs,s)2）选择chebyshev型滤波器阶数数字域：n,Wn=cheb1ord

11、(Wp,Ws,Rp,Rs)模拟域：n,Wn=cheb1ord(Wp,Ws,Rp,Rs,s)3）选择chebyshev型滤波器阶数数字域：n,Wn=cheb2ord(Wp,Ws,Rp,Rs)模拟域：n,Wn=cheb2ord(Wp,Ws,Rp,Rs,s)4)选择椭圆滤波器阶数数字域：n,Wn=ellipord(Wp,Ws,Rp,Rs)模拟域：n,Wn=ellipord(Wp,Ws,Rp,Rs,s)注意：n：返回符合要求性能指标的数字滤波器或模拟滤波器的最小阶数Wn：滤波器的截至频率（即3db频率）Wp：通带的截至频率，Ws：阻带的截至频率，单位rad/s。且均为规一化频率，即。1对应弧度。频率规

12、一化：信号处理工具箱中使用的频率为莱奎斯特频率，根据香农定理，它为采样频率的一半，在滤波器设计中的截止频率均使用莱奎斯特频率进行规一化。规一化频率转换为角频率，则将规一化频率乘以。如果将规一化频率转换为Hz，则将规一化频率乘以采样频率的一半。3模拟频率变换1）低通到低通模拟滤波器变换bt,at=lp2lp(b,a,Wo)lp2lp函数将截至频率为1rad/s（规一化截至频率）的模拟低通滤波器原型变换成截至频率为Wo的低通滤波器。 2）低通到带通模拟滤波器变换bt,at=lp2bp(b,a,Wo,Bw)lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型变换成具有指定带宽Bw和中心频率W

13、o的带通滤波器。其中Wo=sqrt(W1*W2)，Bw=W2-W1 ，W2高端截止频率，W1低端截止频率。3）低通到高通模拟滤波器变换bt,at=lp2hp(b,a,Wo)4）低通到带阻模拟滤波器变换bt,at=lp2bs(b,a,Wo,Bw)5）besself模拟滤波器设计 b,a=besself(n,Wn) b,a=besself(n,Wn,ftype) z,p,k=besself(n,Wn) z,p,k=besself(n,Wn,ftype)2由模拟滤波器变换成等效的数字滤波器方法一：双线性变换法实现模拟滤波器到数字滤波器的变换双线性变换将S域映射成Z域，从而将模拟滤波器变换成等效的数字

14、滤波器zd,pd,kd=bilinear(z,p,k,Fs)为零极点增益表示的bilinear函数其中z,p,k为S域传递函数的零点、极点、和增益，Fs为取样频率，zd,pd,kd为双线性变换后Z域传递函数的零点、极点和增益。numd,dend=bilinear(num,den,Fs,Fp)为传输函数表示的bilinear函数方法二：冲激响应不变法实现模拟到数字滤波器的变换bz,az=impinvar(b,a,Fs)bz,az=impinvar(b,a)采用默认值为1hz的Fs。小结：设计IIR数字滤波器的一般步骤：1）把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标2）根据转换后的性能

15、指标，通过滤波器阶数选择函数，来确定滤波器的最小阶数N和固有频率Wn3）由最小阶数N得到低通滤波器原型4）由固有频率Wn把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波器5）运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器MATLAB信号处理工具箱提供了几个用于直接设计IIR数字滤波器的函数，这些函数把上面这些复杂的步骤融合为一个整体，为设计滤波器带来了极大的方便。3直接设计IIR数字滤波器1Butterworth模拟和数字滤波器设计数字域：b,a=butter(n,Wn)可设计出截止频率为Wn的n阶butterworth滤波器 b,a=butter(n,Wn,ftype)当f

16、type=high时，可设计出截止频率为Wn的高通滤波器；当ftypestop时，可设计出带阻滤波器 z,p,k=butter(n,Wn) zp,k=buter(n,Wn,ftype) A,B,C,D=butter(n,Wn) A,B,C,D=butter(n,Wn,ftype)模拟域：b,a=butter(n,Wn,s)可设计出截止频率为Wn的n阶模拟butterworth滤波器， 其余形式类似于数字域的。2chebyshev型滤波器（通带等波纹）设计数字域：b,a=cheby1(n,Rp,Wn)可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定 b,a=c

17、heby1(n,Rp,Wn,ftype)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftypestop时，可设计出带阻滤波器 z,p,k=cheby1(n,Rp,Wn) zp,k= cheby1 (n,Rp,Wn,ftype) A,B,C,D= cheby1 (n,Rp,Wn) A,B,C,D= cheby1 (n,Rp,Wn,ftype)模拟域：b,a= cheby1 (n,Rp,Wn,s)可设计出截止频率为Wn的n阶chebyshevI型模拟滤波器，其余形式类似于数字域的。3chebyshev型滤波器（阻带等波纹）设计数字域：b,a=cheby2(n,Rs,Wn)可设计

18、出n阶chebyshevI滤波器,其截止频率由Wn确定，阻带内的波纹由Rs确定 b,a=cheby2(n,Rs,Wn,ftype)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftypestop时，可设计出带阻滤波器 z,p,k=cheby2(n,Rs,Wn) zp,k= cheby2(n,Rs,Wn,ftype) A,B,C,D= cheby2 (n,Rs,Wn) A,B,C,D= cheby2(n,Rs,Wn,ftype)模拟域：b,a= cheby2(n,Rs,Wn,s)可设计出截止频率为Wn的n阶chebyshev型模拟滤波器，其余形式类似于数字域的。4椭圆滤波器设计

19、数字域：b,a=ellip(n,Rp,Rs,Wn)可设计出n阶椭圆滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定,阻带内的波纹由Rs确定。 b,a=ellip(n,Rp,Rs,Wn,ftype)当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftypestop时，可设计出带阻滤波器 z,p,k=ellip(n,Rp,Rs,Wn) zp,k=ellip(n,Rp,Rs,Wn,ftype) A,B,C,D=ellip(n,Rp,Rs,Wn) A,B,C,D=ellip(n,Rp,Rs,Wn,ftype)模拟域：b,a=ellip(n,Rp,Rs,Wn,s)可设计出截止频率为W

20、n的n阶椭圆型模拟滤波器，其余形式类似于数字域的。5基于幅度特性设计IIR数字滤波器b,a=yulewalk(n,f,m)可得到以b,a表示的n阶IIR滤波器，其幅频特性逼近由频率矢量f和幅度矢量m所给定的特性。矢量f和n的长度必须相同。6基于频率特性设计IIR数字滤波器b,a=invfreqz(h,w,nb,na)可得到以b,a表示的IIR滤波器，其幅频特性逼近由频率矢量w和复频率响应矢量h所给定的特性。矢量w和h的长度必须相同，nb,na分别为指定的分子和分母多项式的阶数。四、实验内容和步骤1、IIR数字滤波器经典设计法 IIR数字滤波器经典设计法的一般步骤为：（1）根据给定的性能指标和方

21、法不同，首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。（2）估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。（3）设计模拟低通滤波器原型。利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。（4）由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。（5）将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bil

22、inear或impinvar。后面的几个步骤在前面的章节中已经论述，这里主要介绍第一步，关于设计性能指标的转换。设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。数字性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。数字频率和的取值范围为0，单位弧度。而MATLAB工具函数常采用归一化频率，和的取值范围为01，对应于0，此时需进行转换。模拟性能指标给出通带截止频率，阻带起始频率，通带波纹Rp，阻带衰减Rs等。模拟频率和单位为弧度/秒（rad/s）。MATLAB信号处理工具箱中，设计性能指标的转换应根据不同设计方法进行不同处理。下面就举例介绍这些方法。冲激（脉冲

23、）响应不变法采用来将和变换为和【例1】用脉冲响应不变法设计一个Butterworth低通数字滤波器，使其特征逼近一个低通Butterworth模拟滤波器的下列性能指标：通带截止频率，通带波纹Rp小于3dB，阻带边界频率为，阻带衰减大于15dB，采样频率Fs=10000Hz。假设一个信号，其中f1=1000Hz,f2=4000Hz。试将原信号与通过该滤波器的输出信号进行比较。%Samp6_5Wp=2000*2*pi;Ws=3000*2*pi; %滤波器截止频率Rp=3;Rs=15; %通带波纹和阻带衰减Fs=10000; %采样频率Nn=128; %调用freqz所用的频率点数N,Wn=butt

24、ord(Wp,Ws,Rp,Rs,s); %模拟滤波器的最小阶数z,p,k=buttap(N); %设计模拟低通原型Butterworth滤波器Bap,Aap=zp2tf(z,p,k); %将零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap,Wn); %进行频率转换bz,az=impinvar(b,a,Fs); %运用脉冲响应不变法得到数字滤波器的传递函数figure(1)H,f=freqz(bz,az,Nn,Fs); %求解数字滤波器的幅频特性和相频特性subplot(2,1,1),plot(f,20*log10(abs(H)xlabel(频率/Hz);ylabel(振幅/

25、dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(2)f1=1000;f2=4000; %输入信号的频率N=100; %数据长度dt=1/Fs;n=0:N-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(输入信号) %绘制输入信号y=filtfilt(bz,az,x); %用函数filtfilt对输入信号进行

26、滤波y1=filter(bz,az,x); %用filter函数对输入信号滤波subplot(2,1,2),plot(t,y,t,y1,:),title(输出信号),xlabel(时间/s)legend( filtfilt , filter) %加图例由图可知，在小于2000Hz处的衰减小于3dB，而在大于3000Hz处衰减大于15dB，满足滤波器的设计指标。由图6-5可见滤波器对含有1000Hz和4000Hz频率成分的信号进行了滤波，滤除了4000Hz的信号。由程序的输出还可以看出，采用filtfilt函数，输出的1000Hz信号（实线）与输入1000Hz的信号相位一致，即经过滤波后并没有改

27、变信号波形形状。而运用filter函数滤波后（虚线）有一些延迟，改变了信号的形状。图 6-4 例6-5所设计Butterworth滤波器的频率响应。上图：幅频特性；下图：相频特性图 设计滤波器的输入和输出信号输出信号中实线表示运用filtfilt函数的输出，虚线表示运用filter函数的输出【例2】用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为，通带波纹小于1dB，阻带在内，幅度衰减大于15dB，采样周期T=0.01s。假设一个信号，其中f1=5Hz,f2=30Hz。试将原信号与经过该滤波器的输出信号进行比较。注意该题给出的通带边界频率和阻带边界频率均为数字频率，因此

28、设计时首先要将其转换为模拟频率。%Samp6_6wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15; %数字滤波器截止频率、通带波纹和阻带衰减T=0.01;Nn=128; %采样间隔Wp=wp/T;Ws=ws/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=impinvar(b

29、,a,1/T); %脉冲响应不变法设计数字滤波器传递函数figure(1)H,f=freqz(bz,az,Nn,1/T); %输出幅频响应和相频响应subplot(2,1,1),plot(f,20*log10(abs(H);xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;figure(2)f1=5;f2=30; %输入信号含有的频率N=100; %数据点数n=0:N-1;t=n*T; %时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号subplot(2,1,1),plot(t,x),title(输入信号)y=filtfilt(bz,az,x); %对信号进行滤波subplot(2,1,2),plot(t,y),title(输出信号),xlabel(时间/s) 该例所要求的通带频率为，而该滤波器的采样间隔为0.01s，采样频率为100Hz，即2对应于100Hz，则0.2对应于10Hz，即该滤波器的通带范围为010Hz。观看图，其通带范围最大衰减小于1dB，与该分析一致。题意要求阻带在内，幅度衰减大于15dB。其中0.3对应于15Hz，幅度衰减大于15