基于Matlab的IIR数字滤波器设计.docx
《基于Matlab的IIR数字滤波器设计.docx》由会员分享,可在线阅读,更多相关《基于Matlab的IIR数字滤波器设计.docx(9页珍藏版)》请在冰豆网上搜索。
基于Matlab的IIR数字滤波器设计
基于Matlab的IIR数字滤波器设计
摘要:
介绍了基于Matlab的IIR数字滤波器设计方法。
先确定数字滤波器的性能指标,再按照一定的映射规则(冲激响应不变法或双线性变换法)变换成模拟滤波器的性能指标,然后采用一定的逼近方法(巴特沃斯型或切比雪夫型)设计模拟滤波器,最后将模拟滤波器按照同样的映射规则转变成数字滤波器。
同时介绍了设计IIR数字滤波器常用的Matlab函数。
通过Matlab实验仿真,利用介绍的数字滤波器的设计方法,成功地设计出了满足预定指标的各型IIR数字滤波器。
关键词:
Matlab;IIR数字滤波器;设计
0引言
常用的数字滤波器主要有两种,无限长单位冲激响应IIR滤波器和有限长单位冲激响应FIR滤波器。
其中IIR数字滤波器主要有两种设计方法:
①利用模拟滤波器的设计资源。
先设计一个合适的模拟滤波器,然后变换成满足预定指标的数字滤波器。
这种方法比较方便,因为模拟滤波器具有很多现成的设计公式,并且设计参数已经表格化,设计起来既方便又准确;②最优化设计方法。
先确定一种最优准则,如实际频率响应幅度与理想频率响应幅度的均方误差最小准则,或是它们的最大误差最小准则等,然后求此准则下滤波器系统函数的系数ai,bi。
这种方法需要进行大量的迭代运算,所以离不开计算机。
本文主要以设计IIR数字低通滤波器为例,介绍基于Matlab的IIR数字滤波器设计方法,其中采用的是利用模拟滤波器设计资源的方法。
1利用模拟滤波器设计IIR数字滤波器的步骤
(1)先将给定的数字滤波器的性能指标,按照某一变换(冲激响应不变法、双线性变换法等)规则转换成相应的模拟滤波器的性能指标。
(2)若设计的不是数字低通滤波器,还需将步骤
(1)变换得到(高通、带通、带阻)模拟滤波器的性能指标转变成模拟低通滤波器的性能指标,因为只有模拟低通滤波器才有图表资源可以利用。
(3)根据得到的模拟低通滤波器的性能指标,利用某种模拟滤波器的逼近方法(巴特沃斯滤波器、切贝雪夫滤波器、椭圆型滤波器、贝塞尔滤波器等),设计并查表求得此模拟低通滤波器的系统函数。
(4)利用与步骤
(1)和步骤
(2)中的同一变换规则,将模拟低通滤波器的系统函数最终转变成所需的数字各型滤波器的系统函数。
2基于Matlab设计IIR数字滤波器的步骤
以设计IIR数字低通滤波器为例,给定滤波器的性能指标:
设计一个数字低通滤波器,通带纹波(最大衰减)δp=1dB,阻带最小衰减δs=25dB,通带截止频率ωp=0.2π,阻带截止频率ωs=0.4π。
根据设计要求,模拟滤波器可以采用巴特沃斯型、切贝雪夫型、椭圆型、贝塞尔型等,而本文介绍巴特沃斯型、切贝雪夫I型。
模拟滤波器到数字滤波器的映射方法可以采用冲激响应不变法或双线性变换法。
2.1确定采样周期T
因为设计的是数字滤波器,必须确定采样周期T。
为了方便计算,可以确定T=1s,实验表明T取1s或0.1s时,设计出的滤波器的数据和图形没什么差别,本文选择的是T=0.05s。
2.2确定数字滤波器的性能指标
ωp=0.2π,ωs=0.4π,δp=1dB,δs=25dB2.3转变成模拟滤波器的性能指标
(1)模拟角频率Ω和数字角频率ω之间的对应关系。
冲激响应不变法:
Ω=ω12T
(1)
两角频率之间呈线性关系。
双线性变换法:
Ω=c·tanω122
(2)
通常使模拟滤波器与数字滤波器在低频处有较好的对应关系,所以c常取2/T。
通过关系式
(1)、
(2)可得模拟滤波器的通带截止频率ΩP和阻带截止频率Ωs。
(2)模拟滤波器的通带波纹与阻带最小衰减仍为δp和δs。
2.4根据模拟滤波器的性能指标,计算出模拟滤波器的阶数N以及特征通带截止频率Ωc2.4.1巴特沃斯滤波器
巴特沃斯低通滤波器幅度平方函数定义为:
|Ha(jΩ)|2=1121+Ω12Ωc2N(3)
其中,N为滤波器的阶数,Ωc是3dB衰减频率。
根据指标可以求得:
N=1g10δs/10-11210δp/10-121gOs12Op(4)
Ωc=Ωp(100.1δp-1)-1122N或Ωc=Ωs(100.1δp-1)-1122N(5)
在Matlab中,可以直接调用函数buttord,具体格式为:
[N,Ωc]=buttord(Ωp,Ωs,δp,δs,′s′)(6)
根据模拟滤波器的性能指标Ωp、Ωs、δp和δs,计算巴特沃斯模拟滤波器的阶数N和3dB衰减频率Ωc。
2.4.2切比雪夫滤波器
以切比雪夫Ⅰ型低通滤波器为例,它的幅度平方函数定义为:
|Ha(jΩ)|2=1121+ε2C2N(Ω12Ωc)(7)
其中,ε为小于1的正数,表示通带波动的程度,ε越大,波纹越大;Ωc为通带截止频率,衰减为210g1+ε2dB;CN(x)是N阶切比雪夫多项式,N也是滤波器的阶数。
CN(x)=cos(Narccosx)|x|≤1
Ch(Narcchx)|x|>1(8)
根据指标可以求得:
N=arch100.1δs-112100.1δp-1arch(Ωs/Ωp)(9)
ε2=10δp/10-1(10)与巴特沃斯滤波器器不同的是,切比雪夫滤波器中的Ωc=Ωp。
在Matlab中,可以直接调用函数cheb1ord。
具体格式为:
[N,Ωc]=cheblord(Ωp,Ωs,δp,δs,′s′)(11)
根据模拟滤波器的性能指标Ωp、δp、和δs,计算切比雪夫I型模拟滤波器的阶数N和特征通带截止频率Ωc。
2.5根据阶数N、特征通带截止频率Ωc或者通带波纹δp,计算出模拟低通滤波器根据N、Ωc或者δp,通过查表法,得到模拟低通滤波器的系统函数Han(s),再对Han(s)进行反归一化ss12Ωc,就可得到模拟低通滤波器。
Ha(s)=Han(s)|ss12Ωc=Han(s12Ωc)(12)
在Matlab中,可以直接调用函数butter或cheby1,具体格式为:
[b,a]=butter(N,Ωc,′s′)(13)
[b,a]=cheby1(N,δp,Ωc,′s′)(14)
根据阶数N、特征通带截止频率Ωc或者通带波纹δp,计算出模拟低通滤波器的系统函数,cs和ds分别是系统函数分子和分母的系数。
Ha(s)=B(s)12A(s)=b
(2)sn+b
(2)sn-1+…+b(n+1)12a
(1)sn+a
(2)sn-1+…+a(n+1)(15)
2.6将模拟滤波器映射成数字滤波器
2.6.1冲激响应不变法
将模拟滤波器的系统函数写成:
Ha(s)=N12k=1Ak12s-sk(16)
则数字滤波器的系统函数写成:
H(z)=N12k=1TA121-eskTz-1(17)
在Matlab中,可以直接调用函数impinvar,具体格式为:
[d,c]=impinvar(b,a,1/T)。
根据模拟滤波器的系统函数和采样频率,利用冲激响应不变法计算数字滤波器的系统函数:
H(z)=B(s)12A(s)=b
(1)+d
(2)z-1+…+d(n+1)z-n12c
(1)+c
(2)z-1+…+c(n+1)z-n(18)
冲激响应不变法使得数字滤波器的单位冲激响应能完全模仿模拟滤波器的单位冲激响应,时域逼近良好,而且模拟角频率Ω和数字角频率ω之间呈线性关系ω=ΩT。
该方法最大的缺点是有频率响应的混叠效应,所以只适用于限带的模拟滤波器(例如,衰减特性很好的低通或带通滤波器),而且阻带衰减越快,混叠效应越小。
2.6.2双线性变换法
双线性变换法中模拟滤波器的系统函数与数字滤波器的系统关系为:
H(z)=Ha(s)|s=c1-z-1121+z-1=Hac1-z-1121+z-1(19)
其中,c取2/T
在Matlab中,可以直接调用函数bilinear,具体格式为:
[d,c]=bilinear(b,a,1/T);
根据模拟滤波器的系统函数和采样频率,利用双线性不变法计算数字滤波器的系统函数。
双线性变换最突出的优点是避免了频率响应的混叠失真,缺点是频率响应的非线性失真,模拟角频率Ω和数字角频率ω之间的关系为Ω=c·tanω122。
在零频率附近,Ω与ω之间的关系近似于线性,随着Ω的增加,Ω与ω之间的关系出现严重非线性,使数字滤波器频率响应不能保真地模仿模拟滤波器频率响应。
双线性变换法的非线性关系要求模拟滤波器的幅频响应必须是分段常数型的,否则变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有较大畸变。
2.7计算数字滤波器的频率响应
数字滤波器的频率响应与系统函数的对应关系为:
H(ejω)=H(z)|z=ejω(20)
数字滤波器的离散频谱与连续频谱和系统函数的对应关系为:
H(k)=H(ejω)|w=2π12Nk=H(z)|z=ej2π12Nk0≤k≤N-1(21)
在Matlab中,可以直接调用函数freqz,具体格式为:
[H,w]=freqz(d,c,m)(22)
计算数字滤波器的频率响应,在0~π的范围内采样m个点,H为频率响应向量,w为对应的频率向量。
可以进一步计算频率响应的幅值、幅值的分贝形式、相位响应、群延时响应:
mag=abs(H)(23)
db=20*log10((mag+eps)/max(mag))(24)
pha=angle(H)(25)
grd=grpdelay(b,a,w)(26)
3Matlab中IIR数字滤波器设计结果
根据上述的设计步骤,可以设计出预定指标为ωp=0.2π,ωs=0.4π,δp=1dB和δs=25dB的IIR数字低通滤波器。
3.1采用巴特沃斯型模拟滤波器以及冲激响应不变法的设计结果数字滤波器的系统函数为:
N=6
b=[-4.9738e-015,0.001099,0.016648,0.025299,0.0061448,0.00014839,0]
a=[1,-3.1025,4.4121,-3.5511,1.6817,-0.44031,0.049487]
图1所示4张图分别是数字滤波器的幅频响应、幅频响应的分贝形式、相频响应以及群延时。
通过前两张图可见,设计结果满足预定指标。
通过后两张图可见,IIR滤波器的相位为非线性。
3.2采用切比雪夫I型模拟滤波器以及双线性变换法的设计结果数字滤波器的系统函数为:
N=3
b=[0.011475,0.034424,0.034424,0.011475]
a=[1,-2.1378,1.7693,-0.53976]
可见,对于同样的性能指标,用切比雪夫型设计的滤波器比巴特沃斯型设计的滤波器阶数低。
图1采用巴特沃斯型模拟滤波器以及冲激响应
不变法设计的IIR数字低通滤波器
图2采用切比雪夫型模拟滤波器以及双线性
变换法设计的IIR数字低通滤波器
从图2中的幅频特性曲线可见,设计结果满足预定指标。
3.3IIR数字高通、带通和带阻滤波器的设计结果
在Matlab中设计IIR数字高通、带通和带阻滤波器也很方便,设计步骤与设计低通滤波器类似,除了以下4个函数(6)、(11)、(13)和(14)参数设置有差别,其它函数使用方式不变。
函数(13)和(14)的完整表达是:
[b,a]=butter(N,Ωc,′ftype′,′s′)(27)
[b,a]=cheby1(N,δp,Ωc,′ftype′,′s′)(28)
若参数'ftype'不设置,默认为设计低通滤波器。
3.3.1IIR数字高通滤波器的设计结果
与IIR数字低通滤波器不同的是,设计IIR数字高通滤波器时,函数(6)和(11)中的Ωp>Ωs,并且函数执行后直接得到的是相应的模拟低通滤波器的阶数和特征通带截止频率(以下带通、带阻滤波器也类似)。
随后调用函数(27)和(28),其中参数′ftype′设为′high′,就能得到相应的模拟高通滤波器。
之后再进行冲激响应不变法或双线性变换法得到数字高通滤波器。
例如要设计一个IIR数字高通滤波器,设计指标为:
Ωp=0.6π,Ωs=0.4π,δp=1dB和δS=25dB。
3.3.2IIR数字带通滤波器的设计结果
设计IIR数字带通滤波器时,函数(6)和(11)中的Ωp=[Ωp1,Ωp2],Ωs=[Ωs1,Ωs2],其中Ωs1<[Ωp1<Ωp2<Ωs2],并且函数(27)和(28)中参数′ftype′为′bandpass′。
例如要设计一个IIR数字带通滤波器,设计指标为:
ωp1=0.3π,ωp2=0.4π,ωs1=0.2π,ωs2=0.5π,δp=1dB和δs=25dB。
图3采用切比雪夫型模拟滤波器以及双线
性变换法设计的IIR数字高通滤波器
图4采用巴特沃斯型模拟滤波器以及冲激响应
不变法设计的IIR数字带通滤波器
3.3.3IIR数字带阻滤波器的设计结果
设计IIR数字带阻滤波器时,函数(6)和(11)中的Ωp[Ωp1,Ωp2],Ωs[Ωs1,Ωs2],其中Ωp1<Ωs2<Ωs2<Ωp2,并且函数(27)和(28)中参数′ftype′为′stop′。
例如要设计一个IIR数字带阻滤波器,设计指标为:
ωp1=0.2π,ωp2=0.8π,ωs1=0.4π,ωs2=0.7π,δp=1dB,和δs=25dB。
图5采用巴特沃斯型模拟滤波器以及双线性变
换法设计的IIR数字带阻滤波器
注意:
双线性变换法可以设计低通、高通、带通和带阻滤波器,而冲激响应不变法只能设计低通和带阻滤波器。
模拟滤波器除了采用巴特沃斯型、切比雪夫型I型以外,还可以采用切比雪夫型II型、椭圆型等,相应求阶数和特征通带截止频率的函数分别为:
cheb2ord,cheby2;ellipord,ellip。
这些函数的使用与巴特沃斯型和切比雪夫型I型类似。
4结语
本文以IIR数字低通滤波器为例,详细介绍了如何在Matlab中进行IIR数字滤波器的设计,该方法同样适用于IIR数字高通、带通和带阻滤波器的设计。
通过Matlab实验仿真,利用本文介绍的设计方法,能有效地设计出各型IIR数字滤波器。
参考文献:
[1]高西全,丁玉美.数字信号处理[M].西安:
西安电子科技大学出版社,2008.
[2]程佩青.数字信号处理教程[M].北京:
清华大学出版社,2007.
[3]吴镇杨.数字信号处理[M].北京:
高等教育出版社,2008.
[4]余成波,陶红艳,杨菁,等.数字信号处理及MATLAB实现[M].第2版.北京:
清华大学出版社,2008.
[5]薛年喜.MATLAB在数字信号处理中的应用[M].第2版.北京:
清华大学出版社,2008.