楊柱中　周激流　嚴(yán)斌宇　郭　輝

摘 要：提出一種在不增加分?jǐn)?shù)階微分濾波器復(fù)雜度的前提下，能有效提高分?jǐn)?shù)階微分濾波器性能的方法。該方法利用幾種基于典型微分算子的分?jǐn)?shù)階微分濾波器之間的互補(bǔ)性，通過相互內(nèi)插結(jié)合的方式，用于提高IIR分?jǐn)?shù)階數(shù)字濾波器的性能。改進(jìn)后的分?jǐn)?shù)階微分濾波器頻率響應(yīng)更接近理想分?jǐn)?shù)階微分濾波器，表明所提方法的有效性。

關(guān)鍵詞：分?jǐn)?shù)階微積分；數(shù)字微分器；IIR濾波器；微分算子；連續(xù)分?jǐn)?shù)擴(kuò)充

中圖分類號(hào)：TN713文獻(xiàn)標(biāo)識(shí)碼：B

文章編號(hào)：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 引 言

分?jǐn)?shù)階微積分是一個(gè)既古老又現(xiàn)代的話題。早在整數(shù)階微積分產(chǎn)生的時(shí)候分?jǐn)?shù)階微積分就產(chǎn)生了，該問題曾被許多數(shù)學(xué)家，如Leibniz(1695)，Euler (1738)，Liouville (1850)，Hardy 和Littlewood(1925) 等涉及和探究過［1］。雖然分?jǐn)?shù)階微積分的研究難度很大，但近三百年在眾多科學(xué)家的不懈努力下，分?jǐn)?shù)階微積分作為純數(shù)學(xué)分支已經(jīng)發(fā)展?jié)u成體系，但其物理意義不明確，阻礙了分?jǐn)?shù)維微積分的應(yīng)用，目前在工程技術(shù)界中沒有得到廣泛應(yīng)用。從Mandelbrot提出分形學(xué)說，將Riemann－Liouville分?jǐn)?shù)階微積分用以分析和研究分形媒介中的布朗運(yùn)動(dòng)以來，分?jǐn)?shù)階微積分才在許多學(xué)科，特別是在化學(xué)、電磁學(xué)、控制學(xué)、材料科學(xué)和力學(xué)中引起廣泛關(guān)注并嘗試著應(yīng)用［2，3］。隨信息科學(xué)的變革和迅猛發(fā)展，分?jǐn)?shù)階運(yùn)算在很多問題的處理過程中所擁有整數(shù)階運(yùn)算無(wú)可比擬的優(yōu)點(diǎn)正逐漸顯露出來。

目前分?jǐn)?shù)階濾波器已經(jīng)在分?jǐn)?shù)階控制器、信號(hào)處理、圖像壓縮和處理等領(lǐng)域得到成功應(yīng)用。分?jǐn)?shù)階數(shù)字分?jǐn)?shù)階微分濾波器的設(shè)計(jì)和改進(jìn)，正成為分?jǐn)?shù)階微積分研究領(lǐng)域的一個(gè)熱點(diǎn)［4－8］。數(shù)字微分濾波器的設(shè)計(jì)方法通?？梢詺w為2類:第一種是線性相位FIR 濾波器方法；另一種是IIR濾波器法?？紤]到濾波器設(shè)計(jì)復(fù)雜度因素，F(xiàn)IR微分濾波器階數(shù)會(huì)受到限制，影響了其頻率響應(yīng)對(duì)理想頻率響應(yīng)的逼近效果［9］，因此這里考慮使用IIR分?jǐn)?shù)階微分濾波器來實(shí)現(xiàn)分?jǐn)?shù)階運(yùn)算。

IIR分?jǐn)?shù)階數(shù)字微分濾波器設(shè)計(jì)的重點(diǎn)是實(shí)現(xiàn)分?jǐn)?shù)階算子的離散化［10］，即是找到一個(gè)函數(shù)Gv(z)，使其頻率響應(yīng)無(wú)限逼近理想分?jǐn)?shù)階數(shù)字微分器的頻率響應(yīng)H璿(ω)=(jω)v?；静襟E可以歸納為：首先，找到頻率響應(yīng)接近理想一階微分的算子；然后基于所選用的微分算子，推導(dǎo)出分?jǐn)?shù)階微分濾波器傳輸函數(shù)；最后通過各種展開方法把傳輸函數(shù)的分?jǐn)?shù)階形式轉(zhuǎn)化為整數(shù)階濾波器形式。完成分?jǐn)?shù)階展開的常用方法有冪級(jí)數(shù)展開（PSE）和連續(xù)分?jǐn)?shù)擴(kuò)充（CFE），其中連續(xù)分?jǐn)?shù)擴(kuò)充方法對(duì)函數(shù)的逼近更好，收斂更快［11］。

首先對(duì)Rectangular算子、Tustin算子、Simpson算子這幾種典型微分算子通過連續(xù)分?jǐn)?shù)擴(kuò)充，得到相應(yīng)的0.5階微分濾波器頻率響應(yīng)。通過分析這幾種算子的頻率響應(yīng)表明，基于這幾種典型算子的分?jǐn)?shù)階微分濾波器各有優(yōu)缺點(diǎn)和具有互補(bǔ)性，將這幾種典型算子進(jìn)行結(jié)合可得到更接近理想分?jǐn)?shù)階微分算子頻率響應(yīng)的算子。

1 典型IIR分?jǐn)?shù)階微分濾波器

1.1 基于Simpson算子的IIR分?jǐn)?shù)階數(shù)字微分濾波器

Simpson微分算子表示為:

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

則Simpson分?jǐn)?shù)階微分器傳輸函數(shù)為：

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

在此使用連續(xù)分?jǐn)?shù)擴(kuò)充(CFE)方法完成對(duì)上式的展開，這里簡(jiǎn)要介紹分?jǐn)?shù)階算子實(shí)現(xiàn)過程中使用到的CFE方法。對(duì)于任何一個(gè)函數(shù)D(z)，可以用下面連續(xù)分?jǐn)?shù)的形式來表示：

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

式中，系數(shù)a璱,b璱是關(guān)于變量z的有理函數(shù)或常數(shù)。只需要通過截?cái)嗖僮?，就能得到有限階逼近函數(shù)。下面列出T=0.001 s時(shí)，使用連續(xù)分?jǐn)?shù)擴(kuò)展（CFE）完成上式的展開，得到0.5階微分的Simpson分?jǐn)?shù)階微分濾波器傳遞函數(shù)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表示微分階數(shù)；n表示濾波器階數(shù)。

圖1是基于Simpson算子的0.5階微分濾波器的頻率響應(yīng)曲線圖。

圖1 Simpson算子分?jǐn)?shù)階微分濾波器的頻率響應(yīng)

通過對(duì)比和分析，從誤差和計(jì)算復(fù)雜度兩個(gè)方面均衡考慮分?jǐn)?shù)階微分濾波器階數(shù)的選為5階比較合適。因此這里濾波器的階數(shù)都選為5階。

1.2 基于Rectangular算子的IIR分?jǐn)?shù)階數(shù)字微分濾波器

Rectangular算子表示為：

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

基于Rectangular算子的分?jǐn)?shù)階微分器傳輸函數(shù)可以寫為：

Gv璕(z)=v(6)

這里使用連續(xù)分?jǐn)?shù)擴(kuò)充(CFE)法將展開上式，實(shí)現(xiàn)對(duì)函數(shù)的有限階逼近。下面列出T=0.001 s時(shí)，0.5階微分Rectangular分?jǐn)?shù)階微分濾波器傳遞函數(shù)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表示微分階數(shù)；n表示濾波器階數(shù)。

1.3 基于Tustin算子的IIR分?jǐn)?shù)階數(shù)字微分濾波器

Tustin算子表示為：

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

基于Tustin算子的分?jǐn)?shù)階微分器傳輸函數(shù)可以寫為：

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

使用連續(xù)分?jǐn)?shù)擴(kuò)充(CFE)方法將上式展開，完成對(duì)函數(shù)的有限階逼近。下面列出了T=0.001 s時(shí)，0.5階微分Tustin分?jǐn)?shù)階微分濾波器傳遞函數(shù)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表示微分階數(shù)；n表示濾波器階數(shù)。

圖2是基于典型Rectangular算子、Tustin算子和Simpson算子的0.5階微分濾波器的頻率特性曲線，所實(shí)現(xiàn)的濾波器階數(shù)都是5階。從圖2中可以看出3種濾波器在低頻區(qū)域，幅度曲線還能與理想幅度一致，但隨著頻率增加，特別是在高頻區(qū)域，誤差迅速增大。

圖2 三種典型算子分?jǐn)?shù)階微分濾波器的頻率特性

從圖2中可以看出，基于Rectangular濾波器的幅度特性最好，但相位特性明顯比另兩種算子的差。Tustin的優(yōu)點(diǎn)在于其相位特性非常好，相位曲線絕大部分區(qū)域都與理想頻率響應(yīng)相位曲線重合。Tustin和Simpson有很強(qiáng)互補(bǔ)性。因?yàn)閮烧咴诘皖l的表現(xiàn)都比較好，雖然在高頻都有明顯誤差，但兩者的幅度曲線分別位于理想頻率曲線的上下兩側(cè)。因此，這里認(rèn)為通過這3種算子的相互結(jié)合，可以得到一種新的、頻率特性更好的微分算子。

2 通過內(nèi)插結(jié)合形成新分?jǐn)?shù)階微分濾波器

2.1 基于Rectangular算子和Tustin算子內(nèi)插結(jié)合的分?jǐn)?shù)階微分濾波器

通過觀察發(fā)現(xiàn)矩形（Rectangular）濾波器和梯形(Tustin)濾波器分別具有最好的幅頻和相頻特性，因此將這兩種濾波器通過內(nèi)插結(jié)合，可獲得更好的近似理想積分器。

由于微分和積分的互逆性，首先推導(dǎo)新的積分算子H瑼(z)。用下標(biāo)A表示結(jié)合后積分器，用下標(biāo)R表示矩形積分器，用下標(biāo)T表示梯形積分器，其積分算子的傳輸函數(shù)由Rectangular算子和Tustin算子按3∶1的比率結(jié)合獲得。積分器傳輸函數(shù)如下所示：

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

代入相應(yīng)的傳遞函數(shù)得：

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

(1/4)(12)

化簡(jiǎn)得：

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

其零點(diǎn)不在單位圓內(nèi)將零點(diǎn)z=-7映射到z=-1/7，通過乘以7對(duì)幅度進(jìn)行相應(yīng)補(bǔ)償，獲得最小相位積分器如下：

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

通過翻轉(zhuǎn)獲得微分器:

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

對(duì)應(yīng)的分?jǐn)?shù)階微分算子Gv瑼(z)為：

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

下面是T=0.001 s時(shí)，使用該算子實(shí)現(xiàn)0.5階微分的IIR分?jǐn)?shù)階微分濾波器傳遞函數(shù)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算子內(nèi)插結(jié)合的分?jǐn)?shù)階微分濾波器

同樣通過觀察發(fā)現(xiàn)Tustin算子和Simpson算子雖然在高頻都有明顯誤差，但兩者的幅度曲線分別位于理想頻率曲線的上下兩側(cè)，以期通過內(nèi)插結(jié)合相互抵消，而獲得性能更好的濾波器。新的積分算子H瑽(z)傳輸函數(shù)通過梯形（Tustin）算子和辛普森（Simpson）算子按2∶3比例結(jié)合獲得。

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

代入相應(yīng)的傳輸函數(shù)化簡(jiǎn)得:

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

通過翻轉(zhuǎn)可以得到相應(yīng)的微分算子

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

式（20）中積分算子的零點(diǎn)為r1=(-3+5)/2和r2=(-3-5)/2，零點(diǎn)r1和r2互為倒數(shù)且r2零點(diǎn)不在單位圓內(nèi)。為了構(gòu)造最小相位系統(tǒng)，將零點(diǎn)r2映射到其倒數(shù)r1上。同時(shí)為了使幅度保持不變，引入補(bǔ)償因子-r2。獲得的積分算子如下:

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

同樣，改進(jìn)得微分算子為:

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

對(duì)應(yīng)的分?jǐn)?shù)階微分算子Gv瑽(z)為：

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

積分算子的極點(diǎn)是1和-1，在單位圓上，不滿足系統(tǒng)穩(wěn)定性，但經(jīng)過后面連續(xù)分?jǐn)?shù)擴(kuò)充方法截?cái)嗪?，可以使極點(diǎn)都在單位圓內(nèi)。

下面是T=0.001 s時(shí)，使用新算子B實(shí)現(xiàn)0.5階微分的IIR分?jǐn)?shù)階微分濾波器函數(shù)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算子內(nèi)插結(jié)合的分?jǐn)?shù)階微分濾波器

同樣將Rectangular算子和Simpson算子結(jié)合也可以形成新算子。新的積分算子H瑿(z)傳輸函數(shù)通過矩形（Rectangular）算子和辛普森（Simpson）算子按5∶3比例結(jié)合獲得：

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

代入相應(yīng)的傳輸函數(shù)化簡(jiǎn)得:

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

新積分算子的零點(diǎn)為r1=(-9+57)/12和r2=(-9-57)/12(r2零點(diǎn)不在單位圓內(nèi))。為了構(gòu)造最小相位系統(tǒng)，將零點(diǎn)r2映射到其倒數(shù)1/r2上。同時(shí)為了使幅度保持不變，引入補(bǔ)償因子-r2。獲得的積分算子如下:

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

同樣，改進(jìn)得微分算子為:

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

積分算子的極點(diǎn)是1和-1，在單位圓上，不滿足系統(tǒng)穩(wěn)定性，但經(jīng)過后面連續(xù)分?jǐn)?shù)擴(kuò)充方法截短后，可以使極點(diǎn)都在單位圓內(nèi)。

下面是T=0.001 s時(shí)，使用新算子C實(shí)現(xiàn)0.5階微分的IIR分?jǐn)?shù)階微分濾波器函數(shù)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顯示的是通過相互結(jié)合的3種新算子的分?jǐn)?shù)階微分濾波器頻率響應(yīng)。可以看出，新算子中A相比B和C具有更好的頻率特性。其幅度特性曲線從低頻到高頻都基本接近理想頻率響應(yīng)曲線。新算子中A的相位特性隨頻率的增大，相位延遲近似線性增加，可以引入分?jǐn)?shù)階延遲濾波器來進(jìn)一步改進(jìn)相位特性。

3 結(jié) 語(yǔ)

主要從頻域角度出發(fā)，對(duì)分?jǐn)?shù)階微分IIR濾波器的設(shè)計(jì)以及實(shí)現(xiàn)進(jìn)行了深入分析。分?jǐn)?shù)階微分IIR濾波器的實(shí)現(xiàn)有兩個(gè)重要的步驟。首先，找到合適的微分算子，所選算子的頻率響應(yīng)逼近理想分?jǐn)?shù)階微分頻率響應(yīng)的程度直接影響到所實(shí)現(xiàn)濾波器的表現(xiàn)；其次，要使用合適的展開方法把傳輸函數(shù)從分?jǐn)?shù)階形式轉(zhuǎn)化成整數(shù)階濾波器的形式，連續(xù)分?jǐn)?shù)擴(kuò)充(CFE)方法是一種廣泛使用并有良好效果的方法。這里通過將幾種典型算子進(jìn)行內(nèi)插結(jié)合獲得了一種整體更接近理想頻率響應(yīng)的算子，使用連續(xù)分?jǐn)?shù)擴(kuò)充(CFE)方法，完成了分?jǐn)?shù)階微分IIR濾波器的數(shù)字實(shí)現(xiàn)，通過新算子頻率響應(yīng)的對(duì)比分析，分?jǐn)?shù)階微分濾波器的性能獲得了明顯的提高。

圖3 三種新算子分?jǐn)?shù)階微分濾波器的頻率響應(yīng)

參考文獻(xiàn)

［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.

作者簡(jiǎn)介

楊柱中 四川大學(xué)博士研究生。研究方向?yàn)閿?shù)字信號(hào)處理及計(jì)算智能等。

周激流 四川大學(xué)教授，博士生導(dǎo)師。研究方向?yàn)橛?jì)算智能、信號(hào)處理以及生物特征識(shí)別等。

嚴(yán)斌宇 博士研究生。研究方向?yàn)閿?shù)字信號(hào)處理及容延遲移動(dòng)傳感器網(wǎng)絡(luò)等。

郭 輝 新疆工業(yè)高等專科學(xué)校副教授。研究方向?yàn)樽詣?dòng)控制及數(shù)字信號(hào)處理等。