一种分数阶微分ＩＩＲ滤波器的设计方法和改进

杨柱中　周激流　严斌宇　郭　辉

摘 要：提出一种在不增加分数阶微分滤波器复杂度的前提下，能有效提高分数阶微分滤波器性能的方法。该方法利用几种基于典型微分算子的分数阶微分滤波器之间的互补性，通过相互内插结合的方式，用于提高IIR分数阶数字滤波器的性能。改进后的分数阶微分滤波器频率响应更接近理想分数阶微分滤波器，表明所提方法的有效性。

关键词：分数阶微积分；数字微分器；IIR滤波器；微分算子；连续分数扩充

中图分类号：TN713文献标识码：B

文章编号：1004 373X(2009)02 054 05

Design Method and Improvement of IIR Fractional Order Differential Filter

YANG Zhuzhong1,ZHOU Jiliu2,YAN Binyu1,GUO Hui3

(1.College of Electronics and Information Engineering,Sichuan University,Chengdu,610064,China;

2.College of Computer Science,Sichuan University,Chengdu,610064,China;

3.Xinjiang Polytechnical College,Urumqi,830091,China)

Abstract：This paper proposes a design method which can improve the performance of fractional order differential filter obviously under the premise of not increasing the structure complexity of the filter.This scheme which bases on mutually compensatory characters of typical differentiator and uses interpolated method to improve performance of IIR digital fractional differential filter.The improved frequency response of fractional differential filter is more close to ideal fractional differential filter,it indicates the validity of proposed method.

Keywords：fractional calculus;digital differentiator;IIR filter;differential operator;continuous fractional expansion

0 引 言

分数阶微积分是一个既古老又现代的话题。早在整数阶微积分产生的时候分数阶微积分就产生了，该问题曾被许多数学家，如Leibniz(1695)，Euler (1738)，Liouville (1850)，Hardy 和Littlewood(1925) 等涉及和探究过［1］。虽然分数阶微积分的研究难度很大，但近三百年在众多科学家的不懈努力下，分数阶微积分作为纯数学分支已经发展渐成体系，但其物理意义不明确，阻碍了分数维微积分的应用，目前在工程技术界中没有得到广泛应用。从Mandelbrot提出分形学说，将Riemann－Liouville分数阶微积分用以分析和研究分形媒介中的布朗运动以来，分数阶微积分才在许多学科，特别是在化学、电磁学、控制学、材料科学和力学中引起广泛关注并尝试着应用［2，3］。随信息科学的变革和迅猛发展，分数阶运算在很多问题的处理过程中所拥有整数阶运算无可比拟的优点正逐渐显露出来。

目前分数阶滤波器已经在分数阶控制器、信号处理、图像压缩和处理等领域得到成功应用。分数阶数字分数阶微分滤波器的设计和改进，正成为分数阶微积分研究领域的一个热点［4－8］。数字微分滤波器的设计方法通常可以归为2类:第一种是线性相位FIR 滤波器方法；另一种是IIR滤波器法。考虑到滤波器设计复杂度因素，FIR微分滤波器阶数会受到限制，影响了其频率响应对理想频率响应的逼近效果［9］，因此这里考虑使用IIR分数阶微分滤波器来实现分数阶运算。

IIR分数阶数字微分滤波器设计的重点是实现分数阶算子的离散化［10］，即是找到一个函数Gv(z)，使其频率响应无限逼近理想分数阶数字微分器的频率响应H璿(ω)=(jω)v。基本步骤可以归纳为：首先，找到频率响应接近理想一阶微分的算子；然后基于所选用的微分算子，推导出分数阶微分滤波器传输函数；最后通过各种展开方法把传输函数的分数阶形式转化为整数阶滤波器形式。完成分数阶展开的常用方法有幂级数展开（PSE）和连续分数扩充（CFE），其中连续分数扩充方法对函数的逼近更好，收敛更快［11］。

首先对Rectangular算子、Tustin算子、Simpson算子这几种典型微分算子通过连续分数扩充，得到相应的0.5阶微分滤波器频率响应。通过分析这几种算子的频率响应表明，基于这几种典型算子的分数阶微分滤波器各有优缺点和具有互补性，将这几种典型算子进行结合可得到更接近理想分数阶微分算子频率响应的算子。

1 典型IIR分数阶微分滤波器

1.1 基于Simpson算子的IIR分数阶数字微分滤波器

Simpson微分算子表示为:

H璖(z)=(3/T)(1)

则Simpson分数阶微分器传输函数为：

G璖(z)=(3/T)(2)

在此使用连续分数扩充(CFE)方法完成对上式的展开，这里简要介绍分数阶算子实现过程中使用到的CFE方法。对于任何一个函数D(z)，可以用下面连续分数的形式来表示：

D(z)礱0(z)+b1(z)a1(z)+b2(z)a2(z)+b3(z)a3(z)+…(3)

式中，系数a璱,b璱是关于变量z的有理函数或常数。只需要通过截断操作，就能得到有限阶逼近函数。下面列出T=0.001 s时，使用连续分数扩展（CFE）完成上式的展开，得到0.5阶微分的Simpson分数阶微分滤波器传递函数Gv璖n(z)：

G0.5璖3(z)=54.77(z3+3.303 3z2+0.262 3z-2.725 4)z3+5.303 3z2+5.868 9z-1.504 1

G0.5璖5(z)=54.77(z5+5.867 4z4+7.346 1z3-6.319 7z2-8.379 6z+2.639 8)z5+7.867 4z4+18.808z3+6.505 1z2-12.395 5z-1.420 7

G0.5璖7(z)=54.77(z7+8.447z6+21.364 31z5+6.495 5z4-32.817 8z3-16.668 6z2+13.875 7z+1.262 4)z7+10.447z6+37.258 2z5+44.777 1z4-14.903 4z3-45.694 6z2-0.007 8z+3.092 4(4)

Gv璖n(z)中v表示微分阶数；n表示滤波器阶数。

图1是基于Simpson算子的0.5阶微分滤波器的频率响应曲线图。

图1 Simpson算子分数阶微分滤波器的频率响应

通过对比和分析，从误差和计算复杂度两个方面均衡考虑分数阶微分滤波器阶数的选为5阶比较合适。因此这里滤波器的阶数都选为5阶。

1.2 基于Rectangular算子的IIR分数阶数字微分滤波器

Rectangular算子表示为：

G璕(z)=(1/T)·(5)

基于Rectangular算子的分数阶微分器传输函数可以写为：

Gv璕(z)=v(6)

这里使用连续分数扩充(CFE)法将展开上式，实现对函数的有限阶逼近。下面列出T=0.001 s时，0.5阶微分Rectangular分数阶微分滤波器传递函数Gv璕n(z)：

G0.5璕3(z)=31.62(z3-1.75z2+0.875z2-0.109 4)z3-1.25z2+0.375z-0.015 6

G0.5璕5(z)=31.62(z5-2.75z4+2.75z3-1.203 1z2+0.214 8z-0.010 7)z5-2.25z4+1.75z3-0.546 9z2+0.058 6z-0.000 976 56

G0.5璕7(z)=31.62(z7-3.75z6+5.625z5-4.296 9z4+1.757 8z3-0.369 1z2+0.032 4z-0.000 915 53)z7-3.25z6+4.125z5-2.578 1z4+0.820 3z3-0.123z2+0.006 8z-0.000 061(7)

Gv璕n(z)中v表示微分阶数；n表示滤波器阶数。

1.3 基于Tustin算子的IIR分数阶数字微分滤波器

Tustin算子表示为：

G璗(z)=(2/T)·(8)

基于Tustin算子的分数阶微分器传输函数可以写为：

Gv璗(z)=(2/T)vv(9)

使用连续分数扩充(CFE)方法将上式展开，完成对函数的有限阶逼近。下面列出了T=0.001 s时，0.5阶微分Tustin分数阶微分滤波器传递函数Gv璗n(z)：

G0.5璗3(z)=44.72(z3-0.5z2-0.5z+0.125)Z3+0.5z2-0.5z-0.125

G0.5璗5(z)=44.72(z5-0.5z4-z3+0.375z2+0.187 5z-0.031 3)z5+0.5z4-z3-0.375z2+0.187 5z+0.031 3

G0.5璗7(z)=44.72(z7-0.5z6-1.5z5+0.625z4+0.625z3-0.187 5z2-0.062 5z+0.007 813)z7+0.5z6-1.5z5-0.625z4+0.625z3+0.187 5z2-0.062 5z-0.007 813(10)

Gv璗n中v表示微分阶数；n表示滤波器阶数。

图2是基于典型Rectangular算子、Tustin算子和Simpson算子的0.5阶微分滤波器的频率特性曲线，所实现的滤波器阶数都是5阶。从图2中可以看出3种滤波器在低频区域，幅度曲线还能与理想幅度一致，但随着频率增加，特别是在高频区域，误差迅速增大。

图2 三种典型算子分数阶微分滤波器的频率特性

从图2中可以看出，基于Rectangular滤波器的幅度特性最好，但相位特性明显比另两种算子的差。Tustin的优点在于其相位特性非常好，相位曲线绝大部分区域都与理想频率响应相位曲线重合。Tustin和Simpson有很强互补性。因为两者在低频的表现都比较好，虽然在高频都有明显误差，但两者的幅度曲线分别位于理想频率曲线的上下两侧。因此，这里认为通过这3种算子的相互结合，可以得到一种新的、频率特性更好的微分算子。

2 通过内插结合形成新分数阶微分滤波器

2.1 基于Rectangular算子和Tustin算子内插结合的分数阶微分滤波器

通过观察发现矩形（Rectangular）滤波器和梯形(Tustin)滤波器分别具有最好的幅频和相频特性，因此将这两种滤波器通过内插结合，可获得更好的近似理想积分器。

由于微分和积分的互逆性，首先推导新的积分算子H瑼(z)。用下标A表示结合后积分器，用下标R表示矩形积分器，用下标T表示梯形积分器，其积分算子的传输函数由Rectangular算子和Tustin算子按3∶1的比率结合获得。积分器传输函数如下所示：

H瑼(z)=(3/4)H璕(z)+(1/4)H璗(z)(11)

代入相应的传递函数得：

H瑼(z)=(3/4)+

(1/4)(12)

化简得：

H瑼(z)=(T/8)(13)

其零点不在单位圆内将零点z=-7映射到z=-1/7，通过乘以7对幅度进行相应补偿，获得最小相位积分器如下：

H瑼(z)=(7T/8)(14)

通过翻转获得微分器:

G瑼(z)=8(z-1)/7T(z+1/7)(15)

对应的分数阶微分算子Gv瑼(z)为：

Gv瑼(z)={(8/7T)}v(16)

下面是T=0.001 s时，使用该算子实现0.5阶微分的IIR分数阶微分滤波器传递函数Gv瑼n(z)：

G0.5瑼3=33.8(z3-1.571 4z2+0.632 7z-0.037 9)z3-z2+0.142 9z+0.020 4

G0.5瑼5=33.8(z5-2.428 6z4+2z3-0.612 2z2+0.038 7z+0.004)z5-1.857 1z4+1.020 4z3-0.122 4z2-0.021 2z+0.001 4

G0.5瑼7=33.8(z7-3.287 5z6+4.102z5-2.376 1z4+0.597 7z3-0.031 2z2-0.006 8z+0.000 334)z7-2.714 3z6+2.632 7z5-1.035z4+0.097 9z3+0.023 7z2-0.002 6z-0.000 108(17)

2.2 基于Tustin算子和Simpson算子内插结合的分数阶微分滤波器

同样通过观察发现Tustin算子和Simpson算子虽然在高频都有明显误差，但两者的幅度曲线分别位于理想频率曲线的上下两侧，以期通过内插结合相互抵消，而获得性能更好的滤波器。新的积分算子H瑽(z)传输函数通过梯形（Tustin）算子和辛普森（Simpson）算子按2∶3比例结合获得。

H瑽(z)=(2/5)H璗(z)+(3/5)H璖(z)(18)

代入相应的传输函数化简得:

H瑽=(2T/5)(19)

通过翻转可以得到相应的微分算子

G瑽(z)=(5/2T)(20)

式（20）中积分算子的零点为r1=(-3+5)/2和r2=(-3-5)/2，零点r1和r2互为倒数且r2零点不在单位圆内。为了构造最小相位系统，将零点r2映射到其倒数r1上。同时为了使幅度保持不变，引入补偿因子-r2。获得的积分算子如下:

H瑽(z)=(-2Tr2/5)(21)

同样，改进得微分算子为:

G瑽(z)=(-5r1/2T)(22)

对应的分数阶微分算子Gv瑽(z)为：

Gv瑽(z)={(-5r1/2T)}v(23)

积分算子的极点是1和-1，在单位圆上，不满足系统稳定性，但经过后面连续分数扩充方法截断后，可以使极点都在单位圆内。

下面是T=0.001 s时，使用新算子B实现0.5阶微分的IIR分数阶微分滤波器函数Gv瑽n(z)：

G0.5瑽3=30.9(z3-0.383 6z2-0.853 5z+0.327 2)z3-0.001 6z2-0.5z+0.000 414

G0.5瑽5(z)=30.9(z5-0.001 6z4-1.499 4z3+0.097 3z2+0.525 4z-0.081 8)z5+0.380 4z4-z3-0.285 3z2+0.187 5z+0.023 8

G0.5瑽7(z)=30.9(z7+0.908 6z6-2.347 1z5-1.361 9z4+1.77z3+0.453z2-0.423 1z+0.020 5)z7+1.290 6z6-1.5z5-1.613 2z4+0.625 1z3+0.484z2-0.062 5z-0.020 2(24)

2.3 基于Rectangular算子和Simpson算子内插结合的分数阶微分滤波器

同样将Rectangular算子和Simpson算子结合也可以形成新算子。新的积分算子H瑿(z)传输函数通过矩形（Rectangular）算子和辛普森（Simpson）算子按5∶3比例结合获得：

H瑿(z)=(5/8)H璕(z)+(3/8)H璖(z)(25)

代入相应的传输函数化简得:

H瑿(z)=(6T/8)(26)

新积分算子的零点为r1=(-9+57)/12和r2=(-9-57)/12(r2零点不在单位圆内)。为了构造最小相位系统，将零点r2映射到其倒数1/r2上。同时为了使幅度保持不变，引入补偿因子-r2。获得的积分算子如下:

H瑿(z)=(-3Tr2/4)(27)

同样，改进得微分算子为:

G瑿(z)=(-4/3Tr2)(28)

积分算子的极点是1和-1，在单位圆上，不满足系统稳定性，但经过后面连续分数扩充方法截短后，可以使极点都在单位圆内。

下面是T=0.001 s时，使用新算子C实现0.5阶微分的IIR分数阶微分滤波器函数Gv瑿n(z)：

G0.5瑿3=31.09(z3-0.350 6z2-0.888 9z+0.335 3)z3+0.072 3z2-0.582 9z+0.030 8

G0.5瑿5=31.09(z5+0.148 6z4-1.696 1z3+0.013 8z2+0.706 1z-0.132 9)z5+0.571 6z4-1.178 8z3-0.405 2z2+0.318 6z+0.004 7

G0.5瑿7=31.09(z7-2.228 4z6-0.803 3z5+3.662 0z4-0.809 1z3-1.415 9z2+0.630 1z+0.043 4)z7-1.805 5z6-1.291 5z5+2.540 5z4+0.203 4z3-0.870 9z2+0.135 8z+0.011 2(29)

图3显示的是通过相互结合的3种新算子的分数阶微分滤波器频率响应。可以看出，新算子中A相比B和C具有更好的频率特性。其幅度特性曲线从低频到高频都基本接近理想频率响应曲线。新算子中A的相位特性随频率的增大，相位延迟近似线性增加，可以引入分数阶延迟滤波器来进一步改进相位特性。

3 结 语

主要从频域角度出发，对分数阶微分IIR滤波器的设计以及实现进行了深入分析。分数阶微分IIR滤波器的实现有两个重要的步骤。首先，找到合适的微分算子，所选算子的频率响应逼近理想分数阶微分频率响应的程度直接影响到所实现滤波器的表现；其次，要使用合适的展开方法把传输函数从分数阶形式转化成整数阶滤波器的形式，连续分数扩充(CFE)方法是一种广泛使用并有良好效果的方法。这里通过将几种典型算子进行内插结合获得了一种整体更接近理想频率响应的算子，使用连续分数扩充(CFE)方法，完成了分数阶微分IIR滤波器的数字实现，通过新算子频率响应的对比分析，分数阶微分滤波器的性能获得了明显的提高。

图3 三种新算子分数阶微分滤波器的频率响应

参考文献

［1］Adam Loverro,Fractional Calculus: History,Definitions and Applications for the Engineer ［EB/OL］.http://www.nd.edu/~msen/Teaching/UnderRes/FracCalc.pdf,2007.

［2］JennSen Leu,Papamarcou A.On Estimating the Spectral Exponent of Fractional Brownian Motion ［J］.IEEE Trans.Information Theory,1995,41(1):233－244.

［3］SzuChu Liu,Shyang Chang.Dimension Estimation of Discrete－time Fractional Brownian Motion with Applications to Image Texture Classification［J］.IEEE Trans.on Image Processing,1997,6(8):1 176－1 184.

［4］Al－Alaoui.Novel Digital Integrator and Differentiator［J］.Electronics Letters,1993,29(4):376－378.

［5］ChienCheng Tseng.Design of FIR and IIR Fractional Order Simpson Digital Integrators［J］.Signal Processing,2006,87(5):1 045－1 057.

［6］ChienCheng Tseng.Digital Integrator Design Using Simpson Rule and Fractional Delay Filter［J］.Image and Signal Processing,2006,153(1):79－86.

［7］Zhao Hui,Qiu Gang,Yao Lirnin.Design of Fractional Order Digital FIR Differentiators Using Frequency Response Approximation［A］.International Conference on Communications,Circuits and Systems.2005(2):27－30.

［8］ChienCheng Tseng.Improved Design of Digital Fractional－order Differentiators Using Fractional Sample Delay［J］.IEEE Trans.On Circuits and Systems I:Regular Papers,2006,5(1):193－203.

［9］ChienCheng Tseng.Design of Fractional Order Digital FIR Differentiator［J］.Signal Processing Letters,2001,8(3):77－79.

［10］Chen Y Q,Moore K L.Discretization Schemes for Fractional－order Differentiators and Integrators［J］.IEEE Trans.on Circuits and Systems,2002,49(3):363－367.

［11］Chen Y Q,Vinagre B M.Continued Fraction Expansion Approaches to Discretizing Fractional Order Derivatives－an Expository Review［J］.Nonlinear Dynamics,2005,38:155－170.

作者简介

杨柱中 四川大学博士研究生。研究方向为数字信号处理及计算智能等。

周激流 四川大学教授，博士生导师。研究方向为计算智能、信号处理以及生物特征识别等。

严斌宇 博士研究生。研究方向为数字信号处理及容延迟移动传感器网络等。

郭 辉 新疆工业高等专科学校副教授。研究方向为自动控制及数字信号处理等。