基于Matlab的IIR数字低通滤波器的计

陕西理工学院课程设计引言当今，数字信号处理技术正飞速发展，它不但自成一门学科，更是以不同的形式渗透到其他学科。它与国民经济息息相关，与国防建设紧密相连；它影响和改变着我们的生产、生活方式，因此受到人们普遍的关注。数字化、智能化和网络化是当代信息技术发展的大趋势，而数字化是智能化和网络化的基础。数字滤波技术是数字信号分析、处理技术的重要分支。无论是信号的获取、传输，还是信号的处理和交换都离不开滤波技术，它对于信号安全可靠和有效灵活地传输是至关重要的。在所有的电子系统中，使用最多技术最复杂的要算数字滤波器了。数字滤波器的优劣直接决定产品的优劣。因此，在教学中，我们也要学会低通滤波器的设计1 IIR数字滤波器1.1数字滤波器的概念滤波器是指用来对输入信号进行滤波的硬件和软件。数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波器可以理解为是一个计算程序或算法将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列并在转化过程中使信号按预定的形式变化。数字滤波实质上是一种运算过程实现对信号运算处理。数字滤波器和模拟滤波器相比因为信号的形式和实现滤波的方法不同数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。输入数字信号数字序列通过特定的运算转变为输出的数字序列因此数字滤波器本质上是一个完成特定运算的数字计算过程也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信滤波器提供运算规则使其按照这个规则完成对输入数据的处理。1.2 IIR数字滤波器结构 IIR数字滤波器的系统函数可以写成封闭函数的形式。IIR数字滤波器采用递归型结构，即结构上带有反馈环路。运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。IIR数字滤波器在计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。1.3 数字低通滤波器的设计原理数字滤波器（Digital Filter，简称DF）是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。DF根据其实现的网络结构或者从单位脉冲响应函数的时域特性分类，可分为两种，即无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。数字滤波器分为有限脉冲响应数字滤波器(FIR)和无限脉冲响应数字滤波器(IIR)。我们可以借助于模拟滤波器的理论和设计方法来设计IIR滤波器的方法，即根据给定的滤波器技术指标，确定滤波器的系统函数，使其尽可能地逼近滤波器的指标，并用MATLAB进行仿真。数字滤波器的最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了相当成熟的技术和方法，有完整的设计公式，还有比较完整的图表可以查询，因此设计数字滤波器可以充分利用这些丰富的资源来进行。 数字滤波器是指输入、输出都是离散时间信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器在数字信号处理中起着非常重要的作用，在信号的过滤、检测与参数的估计等方面，是使用最为广泛的一种线性系统。 为了保证转换后H(z)稳定，必须满足以下条件：第一：因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的，s平面的左半平面映射z平面的单位圆内部；第二：数字滤波器的频率响应模仿模拟滤波器的频响，s平面的虚轴映射z平面的单位圆，相应的频率之间成线性关系。运用数字滤波器的方法有两种：一是采用计算机软件进行，就是把所要完成的工作通过程序让计算机来实现;二是设计专用的数字处理硬件。这个地方主要用到的就是第一种方法。即是用Matlab提供的信号处理工具箱来实现数字滤波器。Matlab信号处理程序提供了丰富的设计方法，可以使得繁琐的程序设计简化成函数的调用，只要以正确的指标参数调用函数，就可以正确快捷地得到设计结果.它们的系统函数分别为： （1.1） （1.2）1.1式中H(z)称为N阶IIR数字滤波器系统函数；1.2式中H(z)称为N-1阶FIR数字滤波器系统函数。典型的模拟低通滤波器的指标如下：分别为通带频率和阻带频率，分别为通带和阻带容限（峰波纹值）。在通带内要求,有时指标由通带最大衰减和阻带最小衰减给出，定义如下： 和第二种常用指标是用参数和A表示通带和阻带要求，如图1.1所示图1.1 二者之间的关系为：和，根据这几个参数可导出另外两个参数d，k，分别称为判别因子和选择性因子。 BUTTERWORTH低通滤波器：幅度平方函数定义为，N为滤波器阶数，为截止频率。当时，有，为3DB带宽。BUTTERWORTH低通滤波器系统函数有以下形式：由模拟滤波器设计IIR数字滤波器，必须建立好s平面和z平面的映射关系。使模拟系统函数变换成数字滤波器的系统函数，通常采用冲激相应不变法和双线性变换法。冲激相应不变法存在频谱混叠现象，双线性变换法消除了这一线象，在IIR数字滤波器的设计中得到了更广泛的应用。 s平面和Z平面的映射关系为,将和待入数字频率和等效的模拟频率之间的映射关系：,由于二者不是线性关系,所以称为预畸变。1.4 方案选择一般滤波器的设计有两种方法：直接法或者间接法，间接法是通过模拟滤波器的设计方法进行的，设计步骤：先设计过渡模拟滤波器得到系统函然后将按照某种方法装换成数字滤波器的系统函数H(z)。由于这种方法已经比较成熟所以我采用种种方法设计：利用巴特沃斯设计数字低通滤波器。成熟的模拟滤波器设计方法主要有脉冲响应不变法和双线性变换法。2.设计结果及分析2.1用脉冲相应不变法设计IIR数字滤波器脉冲响应不变法使得数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应也就是时域逼近良好而且模拟频率和数字频率之间呈线性关系=T。因而，一个线性相位的模拟滤波器。脉冲响应不变法的最大缺点是有频率响应的混叠效应。所以脉冲响应不变法只适用于限带的模拟滤波器(例如衰减特性很好的低通或带通滤波器),而且高频衰减越快,混叠效应越小。至于高通和带阻滤波器,由于它们在高频部分不衰减因此将完全混淆在低频响应中。如果要对高通和带阻滤波器采用脉冲响应不变法就必须先对高通和带阻滤波器加保护滤波器,滤掉高于折叠频率以上的频率然后再使用脉冲响应不变法转换为数字滤波器。当然这样会进一步增加设计复杂性和滤波器的阶数。程序设计：fp=2100;fs=8000; Fs=20000;Rp=0.5;Rs=30;T=1/Fs; %设计指标W1p=fp/Fs*2;W1s=fs/Fs*2;%求归一化频率N,Wn=buttord(W1p,W1s,Rp,Rs,s);%确定butterworth的最小介数N和频率参数Wnz,p,k=buttap(N); %设计模拟低通原型的零极点增益参数bp,ap=zp2tf(z,p,k); %将零极点增益转换成分子分母参数bs,as=lp2lp(bp,ap,Wn*pi*Fs);%将低通原型转换为模拟低通bz,az=impinvar(bs,as,Fs); %用脉冲响应不变法进行模数变换sys=tf(bz,az,T); %给出传输函数H(Z)H,W=freqz(bz,az,512,Fs); %生成频率响应参数subplot(2,1,1);plot(W,20*log10(abs(H); %绘制幅频响应grid on; %加坐标网格xlabel(频率/Hz);ylabel(振幅/dB);subplot(2,1,2);plot(W,abs(H); grid on;xlabel(频率/Hz);ylabel(振幅/H);运行结果：N =4bz = 0.0000 0.0999 0.1914 0.0252az= 1.0000 -1.4336 1.0984 -0.4115 0.0627程序图如图2.1：图2.1 脉冲响应不变法2.2.用双线性变换法设计IIR数字滤波器双线性变换的这个特点是靠频率的严重非线性关系而得到的由于这种频率之间的非线性变换关系就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤器，不再保持原有的线性相位了。其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数，一般典型的低通、高通、带通、带阻型滤波器的响应特性不然变换产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变。2.2.1用双线性变换法设计椭圆数字低通滤波器程序设计：fs=20000;wp=2*pi*2100/fs;ws=2*pi*8000/fs;Rp=0.5;Rs=30;Ts=1/fs;Wp=2/Ts*tan(wp/2);Ws=2/Ts*tan(ws/2); %按频率转换公式进行转换N,Wn=ellipord(Wp,Ws,Rp,Rs,s); %计算模拟滤波器的最小阶数z,p,k=ellipap(N,Rp,Rs);%设计模拟原型滤波器Bap,Aap=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式b,a=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器的频率转化bz,az=bilinear(b,a,fs); %运用双线性变换法得到数字滤波器传递函数H,f=freqz(bz,az,512,fs);subplot(2,1,1);plot(f,20*log10(abs(H);title(N=2 频率响应);grid on;xlabel(频率/Hz);ylabel(振幅/dB);subplot(2,1,2);plot(f,abs(H); grid on;xlabel(频率/Hz);Ylabel（振幅/H）;运行结果N=2bz= 0.1213 0.1662 0.1213az= 1.0000 -0.9889 0.4218图2.2.1双线性变换法的椭圆数字低通滤波器2.2.2.切比雪夫设计数字低通滤波器设计Chebyshev I型和Chebyshev II型数字低通滤波器，要求通带边界频率fp=2.1kHZ通带最大衰减Rp=0.5dB，阻带边界频率fs=8kH，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。Chebyshev I型的M程序如下:Fs=20000; %抽样频率20KHzFlp=2100;Fls=8000;Wp=2*Flp/Fs; %归一化的通带截止频率Ws=2*Fls/Fs; %归一化的阻带截止频率Rp=0.5; %通带最大衰减单位dBRs=30; %阻带最小衰减单位dBN,Wn=cheb1ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率b,a=cheby1(N,Rp,Wn); %返回H(z)的分子分母系数hw,w=freqz(b,a);subplot(2,1,1);plot(w/pi,20*log10(abs(hw);grid on;xlabel(/);ylabel(幅度dB)title(切比雪夫I型幅频响应);subplot(2,1,2);plot(w/pi,abs(hw);grid on;xlabel(/);ylabel(幅度H);程序结果：N=2b= 0.1007 0.2014 0.1007a= 1.0000 -0.9872 0.4140程序图形：图2.2.2.1 双线性变换发的切比雪夫I数字低通滤波Chebyshev II的程序如下：Fs=20000; %抽样频率20KHz Flp=2100; Fls=8000; Wp=2*Flp/Fs; %归一化的通带截止频率 Ws=2*Fls/Fs; %归一化的阻带截止频率 Rp=0.5; %通带最大衰减单位dB Rs=30; %阻带最小衰减单位dB N,Wn=cheb2ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率 b,a=cheby2(N,Rs,Wn); %返回H(z)的分子分母系数 hw,w=freqz(b,a); subplot(2,1,1); plot(w/pi,20*log10(abs(hw);grid on; xlabel(/);ylabel(幅度dB); title(切比雪夫II型幅频响应); subplot(2,1,2);plot(w/pi,abs(hw); grid on; xlabel(/);ylabel(幅度H); 运行结果 N=2 b= 0.2357 0.4241 0.2357 a= 1.0000 -0.2996 0.1950程序响应图2.2.2.2：图2.2.2.2 双线性变换发的切比雪II夫数字低通滤波从频率响应图中可以看出:巴特沃斯滤波器具有单调下降的幅频特性通带内平滑；chebyshev I型滤波器的幅频特性在通带内有波动阻带内单调；chebyshev II型滤波器的幅频特性在阻带内有波动通带内单调；椭圆滤波器的选择性相对前三种是最好的。下降斜度比较大，通带和阻带内均为等波纹同样的性能指标，椭圆滤波器可以最低的阶数来实现。这样根据不同的要求可以选用不同类滤波器。3.总结（1）在实验过程MATLAB函数命令有数字滤波器函数N,Wn=buttord(Wp,Ws,Rp,Rs)、模拟滤波器函数N,Wn=buttord(Wp,Ws,Rp,Rs,s)。式中Wp,Ws分别上通带和阻带的截止频率，实际上它们是归一化频率，其值在01之间，Rp,Rs分别是通带和阻带的衰减，单位为dB。N是求出的相应低通滤波器的阶次，Wn是求出的3dB频率，第二个式子的单位为rad/s，因此，它们实际上是频率。用来设计模拟低通原型滤波器，其调用格式是z,p,k=buttap(N) N是欲设计的低通原型滤波器的阶次，z、p和k分别是设计出的的极点、零点及增益。以下4个文件用来将模拟低通原型滤波器分别转换为低通、高通、带通、及带阻滤波器。其调用格式分别是：（1）B,A=lp2lp(b,a,Wo) 或 B,A=lp2hp(b,a,Wo)（2）B,A=lp2bp(b,a,Wo,Bw) 或 B,A=lp2bs(b,a,Wo,Bw)式中b，a分别是模拟低通原型滤波器有分子、分母多项式的系数向量，B，A分别是转换后的有分子、分母多项