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

1、基于MATLA的数字滤波器的设计1引言数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散线性非时变系统，其输入是一组(由模拟信号取样和量化的 )数字量，其输出是经过变换或说处理的另一组数字量。数字滤波器具 有稳定性高、精度高、灵活性大等突出优点。这里所说的数字滤波器是指理想带通，低通等的频率选择数字滤波器。数字滤波器设计的一个重要步骤是确定一个可实现的传输函数H(z),这个确定传输函数H(z)的过程称为数字滤波器设计。数字滤波器的一般设计过程为：(1)按照实际需要，确定滤波器的性能要求(通常在频域内给定数字滤波的性能要求)。(2)寻找一满足预定性能要求的离散时间线性系统。(3)用有

2、限精度的运算实现所设计的系统。(4)通过模拟，验证所设计的系统是否符合给定性能要求。2数字滤波器的设计滤波器分为两种，分别为模拟滤波器和数字滤波器。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化的过程中，使信号按预定的形式变化。数字滤波器有多种分类，从数字滤波器功能上分可分为低通、高通、带阻、带通滤波器，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应滤波器(IIR )和有限长冲激响应滤波器(FIR)。数字滤波器指标：一般来说，滤波器的幅频特性是分段常数的，以低通为例，在通带内逼近于1,阻带内逼近与0,

3、实际设计的滤波器并非是锐截止的通带和阻带两个范围，两者之间总有一个过渡带。在设计滤波器时事先给定幅频特性允许误差，在通带范围内幅度响应以误差仃1逼近于1,在阻带内幅度响应以误差22逼近于0。1 - -1 | H(ejw)|1,w CpQsPpQs.T设计出模拟滤波器0Pii jL *1 1jfj|,-JIO11-2atan(w)20-2-10-55100w图2数字域频率与模拟域频率的对应关系转化公式推导如下j11jK tan - 三 2j-/2-j- ./2K jsin( /2)二 K e_U_K cos( /2)ej /2 e-j /2j-K1 ej因只关心频率转换,故可设s =jG , z

4、 = ejE,则有s = K1z，(称为双线性变换)1 z(9)所以模拟滤波器转换成数字滤波器的公式为(10)H(z) =Ha(s)由双线性变换公式一 1 Ks可得z=,视为两复平面变换，再由1 - Kss 二；： j11 , z 二 rej可得_ (1Kc)2(K)2r 一 (1 K二)2(K)2(11)从上式可得：仃=0时，r=1, s上虚轴二z上单位圆周。仃0时，r 0时，r 1, s上右半平面 u z上单位圆外。故若模拟滤波器稳定，则双线性变换后数字滤波器也稳定。由于双线性变换法是一种单值映射，因此消除了频率混叠的现象。双线性变换法的缺点是模拟频率与数字频率间的非线性，这种非线性关系要

5、求被变换的连续系统的幅度响应是分段常数型的(某一频率范围内幅度响应近似于某一常数)，不然所映射出的数字频率响应相对于原来的模拟频率响应会产生变形。为解决双线性变换中的频率非线性关系，我们采用预畸的方法，即K =Ktan,其中 K=2/To2在Matlab中利用Mt件bilinear可以对模拟传输函数实行双线性变换法。MATLAB中IIR数字滤波器的设计过程包括两步。第一步，根据给定指标，确定滤波器的阶数 闻口频率缩放因子Wn1第二步，利用这些参数和给定的波纹参数，确定传输函数的关系。阶数估计：利用双线性 变换法设计数字滤波器时，首先要对IIR数字滤波器的阶数进行估计，相应的MC件为：butto

6、rd用于巴特沃斯滤波器，cheb1ord用于切比雪夫1型滤波器，cheb20rd用于切比雪夫2型滤波器，ellipord 用于椭 圆滤波器。滤波器的设计：对于基于双线性变换法的IIR滤波器设计，对应于四种逼近技术(即巴特沃斯、切比雪夫1型和2型及椭圆逼近)，MATLAB：具箱中有相应的函数。特别地可以用到下面的M文件：butter用于巴特沃斯滤波器的设计，cheby1用于切比雪夫1型滤波器的设计，cheby2用于切比雪夫2型滤波器的设计，ellip用于椭圆滤波器的设计。这些函数的输出可以是滤波器传输函数分子和分母的系数 向量，也可以是滤波器的零极点向量和标量增益因子。同时，利用 zp2tf可以

7、由滤波器的零极点向量和 标量增益因子得到传输函数分子和分母的系数向量。相应地，利用函数zp2sos可以得到传输函数分子和分母系数向量的二次项因子。在计算出传输函数的系数之后，可以利用帔彳freqz来计算频率响应。2.2 FIR数字滤波器设计FIR DF的冲激响应h(n)是有限长的，M阶FIR DF可以表示为：M 1y(n)=h(i)x(n -i)(12)i =0其系统函数为：H(z)八 h(n)z(13)n -0与IIR数字滤波器的设计不同，FIR滤波器的设计与模拟滤波器的设计没有任何联系。因此，FIR滤波器的设计基于对指定幅度响应的直接逼近，并通常要求其具有线性相位响应。为了保证滤波器具有线

8、性相位特性，滤波器系数必须满足条件：h(n)= 士 h(M-1-n)。目前关于FIR滤波器的设计方法主要有三种，即窗函数法，频率取样法和切比雪夫等波纹逼近的最优化设计方法。一般应用较多的是第一种和第三种方法。这是因为窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不严格的情况下市比较灵活的。最优化设计法必须借助计算机计算，但是它能得到最佳的等波纹的线性相位FIR滤波器。目前切比雪夫等波纹的线性相位FIR滤波器的计算机机助设计程序已经比较完善，由于采用了REMEZ代算法，所以设计效率也很高，在应用中越来越占优势。(1)窗函数法一般设计过程总是先给定一理想的滤波器频率响应Hd(ejw),然

9、后设计一个FIR滤波器，用它的频M -1d率响应H(ejw)= h(n)ewn来逼近理想的H d(ejw)。这种逼近中最直接的方法， 是在时域中用FIR滤 波器的单位脉冲响成h(n)去逼近理想的单位脉冲响应 hd(n)。因而，先由Hd(ejw)的IDTFT导出hd(n)hd(n) =Hd(ejw)ejwndw (14)2 二-二由于Hd(ejw)是矩形频率特性，故 hd(n) 一定是无限长的序列，且是非因果的。然而 FIR滤波器是有限 长的，所以用有限长的h(n)来逼近无限长的hd(n),最简单的方法是截取 hd (n)中最重要的一段，将无 限长的hd(n)截取成长度为 M勺有限长序列，等效于

10、再 hd(n)上施加了一个长度为 M勺矩形窗口，更为一 般的，可以用一个长度为 M勺窗口函数w(n)来截取hd(n),即h(n) =w(n)%(n)(15)这一方法通常称为窗函数法，窗口函数的形状及长度M勺选择是窗函数法的关键。下面我们一低通为例，了解一下窗函数法的运用：提出希望频率响应函数(低通)|Hd(ej|10n-0.25 0.25 n o图3理想低通滤波器的频响线性相位，具有片断特点，即. e| | .Hd(ej)=| |0 c ：二|，偿二算出%(n) = 51：H(j)jnd。2二一csin(，(n - .) c2(无限长)二(n -)0.5图4理想低通的单位脉冲响应(无限长的一部

11、分)加窗w(n),长N ,得h(n) = hd(n)w(n)(*)要线性相位，就要h(n)关于(N -1)/2偶对称，而hd (n)关于T偶对称，故要求所以要求w(n)关于e=(N 1)/2偶对称.= (N-1)/210.5301020图6加窗后的单位脉冲响应图5窗函数再回过来检验H (ej灯是否满足精度要求.|匕(g)|0.5O若基本满足，则依截取的h(n),制硬件,编软件.H(ej )理想频响与实际频响的对比所示。为便于选择使用，将5种常见的窗函数基本参数如表1表1 5种常见的窗函数基本参数窗函数的 旁瓣峰”过渡带宽度好加窗后滤波器的 阻带最小衰减&rectwin_134兀/ N孕bart

12、let 三角58兄NN5hanning8MN-44hamming-418印N-53blackman-5712必N-74(2)频率取样法窗口设计法事从时域出发，把理想的hd(n)用一定形状的窗口函数截取成有限长的h(n),以此h(n)来近似理想的hd(n),从而频率响应 H (ejw)也近似于理想白频率响应Hd(ejw)。我们知道一个有限长序列可以通过其频谱的相同长度的等间隔采样值准确地恢复原有的序列。频率采样法便是从频域出发， 对理想的频率响应 Hd(ejw)加以等间隔采样Hd(ejw)lw 二叫= Hd(k)wHM(16)然后以此Hd(k)作为实际FIR滤波器的频率特性的离散样本H(k),即

13、H(k) = Hd(k) = Hd(ejw)|Ynk，k=0,1,2,.M -1(17)一|M由H(k)通过IDFT可求出有限长序列h(n)为1 Mj 2 nkh(n)H (k)e M,n=0,1，M -1(18)M y利用M频率的离散样本H(k)同样可求出FIR滤波器的系统函数 H(z)及频率响应H(ejw)。M 1H (z) = h(n)z”n =01-z+M,H(k) M ,1-W,z”-j2 二其中W=e M(19)令z =ejw可得到滤波器的频率响应H(ejw)o如果设计的是线性相位的FIR数字滤波器，其采样值 H(k)的相位的幅度一定要满足特定的约束条件，这个设计时一定要注意。(3

14、)最优化设计法最优化设计法事以最佳一致逼近(最大误差最小化)理论为基础，利用雷米兹算法设计的具有等波 纹特性的设计方法。具体设计步骤如下：对设计指标进行归一化处理。确定remezord函数所需要的参数。 包括归一化边界频率、 各频带的幅度要求和波纹要求等。 归一 化边界频率总是从0开始到1结束，故只需递增列出中间的边界频率； 频带幅度要求不含过渡区， 个数是 边界频率个数的一半加1;波纹要求是频带内幅度允许的波动要求，与分贝间的关系是：1-Rp -20log10(-1), Rs = -2010g10(一)：。唠温口 (20)1 - -11 .1利用remezord函数确定remez所需参数。调

15、用remez函数进行设计。利用freqz函数验算技术指标是否满足要求。2.3数字滤波器类型的选择IIR和FIR各有优缺点，在实际运用中如何选择它们，这里做一个简单的比较。表2 IIR与FIR的比较IIRFIR设计方法利用AF的设计图表，可简单， 功效的完成设计一无解析的设计公式，要借助计算机程序完 成设计结果只能得到幅频特性，相频特性 未知(缺点)，如需要线性相位， 需用全通网络校准，但增加滤 波器的阶数和复杂性可得到幅频特性(可以多带)和线性相位(优 点)稳定性有稳定性问题极点全部在原点(永远稳定)，无稳定性问题因果性总是满足，任何一个非因果的有限长序列，总 可以通过一定的延时，转变为因果序

16、列结构递归系统非递归运算误差有反馈，由于运算中的四舍五 入会产生极限环一无反馈，运算误差小快速算法无快速运算方法可用FFTM少运算量从以上简单的比较可以得到，IIR与FIR滤波器各有所长， 所以应根据实际应用要求，从多方面考虑加以选择。3数字滤波器的MATLAB计(1) IIR的直接程序设计法数字通带边缘频率为 0.25冗和0.8此题的MATLAB序为：ws=0.4*pi 0.7*pi;%wp=0.25*pi 0.8*pi;%rp=1 ;%as=40;%例如欲设计一数字(IIR )带阻滤波器，其数字域指标为：数字阻带边缘频率分别为 0.4冗和0.7冗，n ,通带波动为1db最小阻带衰减为40d

17、b。数字阻带边缘频率数字通带边缘频率通带波动(db)阻带衰减n,wn=cheb20rd(wp/pi,ws/pi,rp,as);%艮据给定指标，确定滤波器的阶数用口频率缩放因子 Wnb,a=cheby2(n,as,ws/pi,stop);% 返回的 b,a 分别为 H (z)的分子、分母。 h,w=freqz(b,a,512);% 返回的h,w分别为滤波器的频率响应及其频率 plot(w/pi,abs(h);%画出频率响应(以w/pi为横轴)grid; xlabel(w/pi);ylabel（幅值）；title（频率响应）;程序运行结果为：图9所设计的带阻滤波器的频率响应在设计中如果该滤波器的特

18、性不满足要求,原有的参数必须做相应的调整，在程序中只需对参数做S为含有3个频新的设定就可以得到所需的滤波器。接下来我们来看看此题所设计的滤波器的滤波效果:率成分的信号（归一化频率（w/2 n ）分别为0.1、0.3、0.45 ）,用所设计的滤波器滤除归一化频率为0.3的成分。n=0:100;s1=sin(pi*0.2*n);s2=sin(pi*0.6*n);s3=sin(pi*0.9*n);s4=s1+s3;s=s1+s2+s3;sf=filter(b,a,s);stem(n,s);title(滤波前的信号）；subplot(311)subplot(312);stem(n,sf);title(

19、滤波后的信号）；subplot(313);stem（n,s4）;title（程序运行的结果为:想要保留的信号）；弓潴避前的博号图10采用filter函数进行数字滤波前后信号比较示意图由图可以看出，滤波后的信号与想要保留的信号基本一致(相位有些许偏差，但基本一致)，所以我们可以说该滤波器基本满足了以上所提出的滤波要求。(2) FIR的直接程序设计法例如欲设计一个线性相位数字(FIR)带通滤波器，其数字域指标为：数字通带边界频率为0.35 n和0.65 n ,数字阻带边界频率为 0.2 n和0.8 n ,通带波动为1db,最小阻带衰减为60db。FIR数字滤波器的窗函数法此题的MATLAB序为：w

20、s1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;as=60;tr=min(wp1-ws1),(ws2-wp2);M=ceil(11*pi/tr)+1;犍波器的阶数，程序运行后M=75n=0:1:M-1;r=(M-1)%曲群时延wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=sin(wc2*(n-r)+eps)./(pi*(n-r)+eps)-sin(wc1*(n-r)+eps)./(pi*(n-r)+eps);%h必理想滤波器的脉冲响应w_bla=(blackman(M);张度为M勺blackman窗h=hd.*w_bla;%h滤

21、波器的实际脉冲响应stem(n,h);title(滤波器的实际单位脉冲响应)；freqz(h,1,512);title(幅度响应和相位响应);法希裔的实际用的肝巾坤.下图11所设计的滤波器的实际单位脉冲响应由上图可知滤波器的实际脉冲响应h是偶对称的，即h(n尸h(M-1-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。一气君富fxlSB应D.10.20.3 0J 0.50.E07 O.B 0.9Normaizitd Frfrqubnt (KErad/sainEltf)怕便哨应和相付响晤D M mJ- mswnnEwz.mm 0 D.l 0.203 0J 0.50.5

22、0.7 D.B O.1Nonnaizad Frfcauencv (mx radisafrcMgi图12所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。接下来我们来看看此题所设计的滤波器的滤波效果：助含有3个频率成分的信号(归一化频率(w/2n )分别为0.05、0.2、0.45 ),用所设计的滤波器滤除归一化频率为0.05和0.45的成分。1=0:100;s1=sin(0.1*pi*1);s2=sin(0.4*pi*1);s3=sin(pi*0.9*1);s=s1+s2+s3;sf=fi1ter(h,1,s);subplot(311)stem(l

23、,s);title( 滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);浦,波前的侑号出L1jIIIII01020300506070 BO 90100夏池后的信号I-.- = b -”e rr- rr- W 203040 SO 6070 EO 9D 100想要保留的信号图13采用filter函数进行数字滤波前后信号比较示意图由上图可知滤波后的信号和想要保留的信号的幅度和频率基本一致(滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器

24、留的群延迟r=(M-1)/2=37 ),所以我们可以说该滤波器基本满足了以上所提出的滤波要求。FIR数字滤波器的频率采样法此题的MATLAB序为：M=40;%取滤波器白阶数为40al=(M-1)/2;哪时延n=0:M-1;T2=0.59417456;T1=0.109021;Hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4); %采样值的幅值k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angH=-al*(2*pi)/M*k1,al*(2*pi)/M

25、*(M-k2);麻样值的相位H=Hrs.*exp(j*angH);h=real(ifft(H,M);%长度为M勺单位脉冲响应stem(n,h);title(滤波器的实际单位脉冲响应);freqz(h,1,512);title(幅度响应和相位响应);图14所设计的滤波器的实际单位脉冲响应由图14可知滤波器的实际脉冲响应h是偶对称的，即h(n尸h(M-1-n),故该滤波器满足FIR线性相位的50福度响应和相位响应条件，该滤波器是线性相位的FIR滤波器。Niwiislizfidi Frequency廊.! 口口 DJ 0.405 O& 0 7 口日 口与 1Nio*Tii?lizfldi Frequ

26、ency (kr radsTfph)娟I 100 *150 0.1 口2 D.30.405 OE 07 0日 口日W 0 -500 -1000图15所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线T相位的 FIR滤波器。此滤波器的群延时为al=(M-1)/2=19.5。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号(归一化频率(w/2兀)分别为0.05、0.2、0.45 ),用所设计的滤波器滤除归一化频率为0.05和0.45的成分。1=0:100;s1=sin(0.1*pi*1);s2=sin(0.4*pi*1);s3=sin(pi*0.9*

27、1);s=s1+s2+s3;sf=fi1ter(h,1,s);subplot(311)stem(l,s);title( 滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);淞淞前的信号1DD20 3D 40 SO 6Q 71tn黄星空帕信号flO 5Q田池波后的信号出 W I小 也.图16采用filter函数进行数字滤波前后信号比较示意图同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致(滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延

28、迟引起的，此滤波器留的群延迟(r=(M-1)/2=19.5 ),所以我们可以说该滤波器基本满足了以上所提出的滤波要求。FIR数字滤波器的最优设计法此题的MATLAB序为：姒计指标ws1=0.2*pi;wp1=0.35*pi;wp2=0.65*pi;ws2=0.8*pi;rp=1;as=60;%设置边界频率和幅度要求F=ws1/pi,wp1/pi,wp1/pi,ws2/pi;A=0,1,0；姒置各频带的波纹要求devp=(10A(rp/20)-1)/(10A(rp/20)+1);devs=10A(-as/20);dev=devs,devp,devs;额定remez参数，其中滤波器的阶数为(N+1

29、),程序运行后得到N=26N,Fo,Ao,W=remezord(F,A,dev);%调用remez函数进行设计h=remez(N,Fo,Ao,W);n=0:N;stem(n,h);title(滤波器的单位冲激响应)；freqz(h,1,512);title(幅度响应和相位响应);：引行行二.力|仃 o -)05*1015202530图17所设计的滤波器的实际单位脉冲响应由图17可知滤波器的实际脉冲响应h是偶对称的，即h(n户h(N-n),故该滤波器满足FIR线性相位的条件,该滤波器是线性相位的FIR滤波器。mM 第 h-aE1王幅度响应和相莅响曲 印 0 国 -1D0Nnmielizad Fn

30、eQuency f陋 m出g&mp博5mRO-1nio o - 02 0.3 o.j as o.B o? ae o.sNarmalizeci! FreouBncy 卜事由必.mpisj图18所设计的带通滤波器的幅度和相位响应由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。此滤波器的群延时为S为含有3个频率成分的信号(归al=(N)/2=13 。接下来我们来看看此题所设计的滤波器的滤波效果:化频率(w/2n)分别为0.05、0.2、0.45 ),用所设计的滤波器滤除归一化频率为0.05和0.45的成分。1=0:100;s1=sin(0.1*pi*1);s2=sin(0.4*pi*1

31、);s3=sin(pi*0.9*1);s=s1+s2+s3;sf=filter(h,1,s);subplot(311)stem(l,s);title(滤波前的信号);subplot(312);stem(l,sf);title(滤波后的信号);subplot(313);stem(l,s2);title(想要保留的信号);滤波前的信号1020 3G 0 SO 60700090100速波后的信号IIIII IIII- c 102030405060700090100想要保雷的信号图19采用filter函数进行数字滤波前后信号比较示意图同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致(滤波后的信号相对于 想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器留的群延迟r=(N)/2=