基于MATLAB的数字滤波器的设计

基于MATLAB的数字滤波器的设计1引言数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散线性非时变系统，其输入是一组（由模拟信号取样和量化的)数字量，其输出是经过变换或说处理的另一组数字量。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。这里所说的数字滤波器是指理想带通，低通等的频率选择数字滤波器。数字滤波器设计的一个重要步骤是确定一个可实现的传输函数H(z)，这个确定传输函数H(z)的过程称为数字滤波器设计。数字滤波器的一般设计过程为：（1）按照实际需要，确定滤波器的性能要求（通常在频域内给定数字滤波的性能要求）。（2）寻找一满足预定性能要求的离散时间线性系统。（3）用有限精度的运算实现所设计的系统。（4）通过模拟，验证所设计的系统是否符合给定性能要求。2数字滤波器的设计滤波器分为两种，分别为模拟滤波器和数字滤波器。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化的过程中，使信号按预定的形式变化。数字滤波器有多种分类，从数字滤波器功能上分可分为低通、高通、带阻、带通滤波器，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应滤波器（IIR）和有限长冲激响应滤波器（FIR）。数字滤波器指标：一般来说，滤波器的幅频特性是分段常数的，以低通为例，在通带内逼近于1，阻带内逼近与0，实际设计的滤波器并非是锐截止的通带和阻带两个范围，两者之间总有一个过渡带。在设计滤波器时事先给定幅频特性允许误差，在通带范围内幅度响应以误差1逼近于1，在阻带内幅度响应以误差2逼近于0。wweHwweHrjwcjw,2|)(|,1|)(|11（1）式中wc和wr分别为通带边界频率和阻带边界频率，wr-wc为过渡带。在具体的技术指标中往往用通带波动来表示1，用最小阻带衰减At来表示2，其具体的对应公式这里就不详述了。2.1IIR数字滤波器设计IIRDF的冲激响应h(n)是无限长的，其输入输出的关系为：)()()(inxihnyi(2)系统函数为nnnznhzH)()(=nkkkmrrrzazb101(3)设计无限长单位脉冲响应（IIR）数字滤波器一般可有三种方法。第一种方法，先设计一个合适的模拟滤波器，然后将其数字话，即将S平面映射到Z平面得到所需的数字滤波器。模拟滤波器的设计技巧非常成熟，不仅得到的是闭合形式的公式，而且设计系数已经表格化了。因此，由模拟滤波器设计数字滤波器的方法准确，简便，得到普遍采用。对于这种方法，工程上有两种常见得变换法——脉冲响应不变法及双线性变换法。第二种方法，在Z平面直接设计IIR数字滤波器，给出闭合形式的公式，或者以所希望的滤波器响应作为依据，直接在Z平面上通过多次选定极点和零点的位置，以逼近该响应。第三种方法，利用最优化技术设计参数，选定极点和零点在Z平面上的合适位置，在某种最优化准则意义上逼近所希望的响应。但一般不能得到滤波器的系数（即零，极点的位置）作为给定响应的闭合形式函数表达式。优化设计需要完成大量的迭代运算，这种设计法实际上也是IIR滤波器的直接设计。本文着重介绍由模拟滤波器设计相应的IIR数字滤波器的方法。（1）脉冲响应不变法脉冲响应不变法是使数字滤波器的单位脉冲响应序列h(n)逼近模拟滤波器的冲激响应()aht，让h(n)正好等于()aht的采样值。设已有满足要求的模拟滤波器,则可()aHs()aht()()ahnhnT()Hz{因为:()()ahthn的图形的图形},公式导出:具体转换如下:设(以一阶极点为例)1()NkakkAHsss(4)作拉氏反变换,得11()[()]()kNstaakkhtFTHsAeut采样得1()()()kNsnTakkhnhnTAeunT作Z变换,得1011()()1kkNNsnTnkksTnkkAHzAeunTzez(5)()Hz与()aHs极点关系为:ksTkze(6)一般对应关系sTjTjTsjjzrezereee,TreT(7)图1S平面到Z平面变换示意图所以,模拟系统稳定因果数字系统稳定因果。按照脉冲响应不变法，从S平面到Z平面的映射不是单值关系，而是先将()aHs在S平面沿虚轴作周期严拓，再按照映射关系将()aHs映射到Z平面，得到()Hz，因此，脉冲响应不变法只适用于带限的滤波器（如低通、带通）。在Matlab中利用M文件impinvar可以对模拟传输函数实行脉冲响应不变法。（2）双线性变换法脉冲响应不变法不适带阻和高通滤波器的设计，因为高频带为通带，前述方法易引起混频。故希望:s平面虚轴z平面单位圆一周,且应有:0:0,:0:0,因为tan(/2),所以选变换tan2K(8)其中K可取任意正常数,后面将导出2/KT.设计思路:tan2,,,,,,pspspspsK.设计出模拟滤波器转化成数字滤波器.图2数字域频率与模拟域频率的对应关系转化公式推导如下:sin(/2)tan2cos(/2)jjjKK/2/2/2/211jjjjjjeeeKKeee因只关心频率转换,故可设sj,jze,则有1111zsKz,(称为双线性变换)(9)所以模拟滤波器转换成数字滤波器的公式为jO11rOs平面00z平面1rTrekskz-10-50510-202w2atan(w)ppO1111()()zasKzHzHs(10)由双线性变换公式,可得11KszKs,视为两复平面变换,再由,jsjzre.可得2222(1)()(1)()KKrKK(11)从上式可得:0时,1r,s上虚轴z上单位圆周。0时,1r,s上左半平面z上单位圆内。0时,1r,s上右半平面z上单位圆外。故若模拟滤波器稳定，则双线性变换后数字滤波器也稳定。由于双线性变换法是一种单值映射，因此消除了频率混叠的现象。双线性变换法的缺点是模拟频率与数字频率间的非线性，这种非线性关系要求被变换的连续系统的幅度响应是分段常数型的（某一频率范围内幅度响应近似于某一常数），不然所映射出的数字频率响应相对于原来的模拟频率响应会产生变形。为解决双线性变换中的频率非线性关系，我们采用预畸的方法，即tan2K，其中K=2/T。在Matlab中利用M文件bilinear可以对模拟传输函数实行双线性变换法。MATLAB中IIR数字滤波器的设计过程包括两步。第一步，根据给定指标，确定滤波器的阶数N和频率缩放因子Wn。第二步，利用这些参数和给定的波纹参数，确定传输函数的关系。阶数估计：利用双线性变换法设计数字滤波器时，首先要对IIR数字滤波器的阶数进行估计，相应的M文件为：buttord用于巴特沃斯滤波器，cheb1ord用于切比雪夫1型滤波器，cheb2ord用于切比雪夫2型滤波器，ellipord用于椭圆滤波器。滤波器的设计：对于基于双线性变换法的IIR滤波器设计，对应于四种逼近技术（即巴特沃斯、切比雪夫1型和2型及椭圆逼近），MATLAB工具箱中有相应的函数。特别地可以用到下面的M文件：butter用于巴特沃斯滤波器的设计，cheby1用于切比雪夫1型滤波器的设计，cheby2用于切比雪夫2型滤波器的设计，ellip用于椭圆滤波器的设计。这些函数的输出可以是滤波器传输函数分子和分母的系数向量，也可以是滤波器的零极点向量和标量增益因子。同时，利用zp2tf可以由滤波器的零极点向量和标量增益因子得到传输函数分子和分母的系数向量。相应地，利用函数zp2sos可以得到传输函数分子和分母系数向量的二次项因子。在计算出传输函数的系数之后，可以利用M文件freqz来计算频率响应。2.2FIR数字滤波器设计FIRDF的冲激响应h(n)是有限长的，M阶FIRDF可以表示为：10)()()(Miinxihny(12)其系统函数为:10)()(MnnznhzH(13)与IIR数字滤波器的设计不同，FIR滤波器的设计与模拟滤波器的设计没有任何联系。因此，FIR滤波器的设计基于对指定幅度响应的直接逼近，并通常要求其具有线性相位响应。为了保证滤波器具有线性相位特性，滤波器系数必须满足条件：h(n)=h(M-1-n)。目前关于FIR滤波器的设计方法主要有三种，即窗函数法，频率取样法和切比雪夫等波纹逼近的最优化设计方法。一般应用较多的是第一种和第三种方法。这是因为窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不严格的情况下市比较灵活的。最优化设计法必须借助计算机计算，但是它能得到最佳的等波纹的线性相位FIR滤波器。目前切比雪夫等波纹的线性相位FIR滤波器的计算机机助设计程序已经比较完善，由于采用了REMEZ迭代算法，所以设计效率也很高，在应用中越来越占优势。（1）窗函数法一般设计过程总是先给定一理想的滤波器频率响应)(jwdeH，然后设计一个FIR滤波器，用它的频率响应jwnMnjwenheH10)()(来逼近理想的)(jwdeH。这种逼近中最直接的方法，是在时域中用FIR滤波器的单位脉冲响应h(n)去逼近理想的单位脉冲响应)(nhd。因而，先由)(jwdeH的IDTFT导出)(nhddweeHnhjwnjwdd)(21)((14)由于)(jwdeH是矩形频率特性，故)(nhd一定是无限长的序列，且是非因果的。然而FIR滤波器是有限长的，所以用有限长的h(n)来逼近无限长的)(nhd，最简单的方法是截取)(nhd中最重要的一段，将无限长的)(nhd截取成长度为M的有限长序列，等效于再)(nhd上施加了一个长度为M的矩形窗口，更为一般的，可以用一个长度为M的窗口函数w(n)来截取)(nhd，即)()()(nhnwnhd(15)这一方法通常称为窗函数法，窗口函数的形状及长度M的选择是窗函数法的关键。下面我们一低通为例，了解一下窗函数法的运用:①提出希望频率响应函数（低通）图3理想低通滤波器的频响线性相位,具有片断特点,即||()0||jjcdceHe②算出-1()()d2jjndhnHeeO0.25|()|jdHe0.2511d2ccjjneesin(())()cnn(无限长)图4理想低通的单位脉冲响应（无限长的一部分）③加窗()wn,长N,得()()()dhnhnwn(*)要线性相位,就要()hn关于(1)/2N偶对称,而()dhn关于偶对称,故要求(1)/2N所以要求()wn关于(1)/2N偶对称.010203000.510102030-0.100.10.20.3图5窗函数图6加窗后的单位脉冲响应再回过来检验()jHe是否满足精度要求.图7图4的脉冲响应的频响图8理想频响与实际频响的对比若基本满足,则依截取的()hn,制硬件,编软件.为便于选择使用,将5种常见的窗函数基本参数如表1所示。表15种常见的窗函数基本参数012300.51O0.25|()|jdHe0.251()jHe012300.51010203000.510102030-0.100.10.20.3类型窗函数的旁瓣峰n过渡带宽度B加窗后滤波器的阻带最小衰减srectwin134/N21bartlet三角258/N25hanning318/N44hamming418/N53blackman5712/N74（2）频率取样法窗口设计法事从时域出发，把理想的)(nhd用一定形状的窗口函数截取成有限长的h(n)，以此h(n)来近似