- Email: [email protected]
Two-dimensional fractional-order digital differentiator design by using differential evolution algorithm
Digital Signal Processing 19 (2009) 660–667
Contents lists available at ScienceDirect
Digital Signal Processing www.elsevier.com/locate/dsp
Two-dimensional fractional-order digital differentiator design by using differential evolution algorithm Wei-Der Chang Department of Computer and Communication, Shu-Te University, Kaohsiung 824, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Available online 20 January 2009 Keywords: Two-dimensional digital filter Fractional-order differentiator Differential evolution
Designing a fractional-order digital differentiator often requires considerably complex mathematical operations and numerical approximations. Thus this paper will propose a simple method to achieve the fractional-order digital differentiator design, particularly for two-dimensional fractional-order differentiators. A two-dimensional finite impulse response (FIR) digital filter structure is utilized and designed so that its corresponding magnitude response can satisfy that of a desired fractional-order differentiator of two variables. The algorithm used to design such two-dimensional digital differentiator is the differential evolution (DE), which is one of evolutionary computations and has excellent searching capacity. The efficiency of the proposed scheme can be confirmed by some illustrative examples. © 2009 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the research topic regarding fractional calculus has attracted much attention and successfully introduced to a variety of engineering applications, such as digital signal processing [1–5], chaos systems [6], control system design [7,8], and electrical circuit devices [9]. Generally speaking, the notion of fractional calculus is to enlarge integer order to fractional order in numerical representations. In this paper, we will focus our effort on the issue of fractional-order digital differentiator design, especially for multiple dimensions. About this issue, researchers have presented various design schemes in designing the fractional-order digital differentiator, including time-domain Taylor series expansions [3], Newton series expansions [4], least-squares method in the time domain [5], etc. Basically, the fractional-order digital differentiator is an extended version of general digital differentiator, which generalizes integer order to fractional case, to provide a more flexibility in the real applications. Its role is to give a fractional-order differentiation on a given digital input signal. In comparison with designing the general digital differentiator, designing such a fractional-order case is considerably complicated and difficult. On the other hand, another topic discussed in this paper is about the two-dimensional digital filter design [10–12]. Twodimensional digital filter is to tackle the two-dimensional digital signal, for example, digital image signal. It is obvious that efforts in designing a two-dimensional digital filter are much beyond that in designing a one-dimensional case. The reasons are as follows [10]: 1. Two-dimensional digital filter basically involves more data than one-dimensional case. 2. Mathematics for tackling the two-dimensional digital filter is more incomplete than that for tackling the onedimensional case. 3. In the design, two-dimensional digital filter has more degrees of freedom to be evaluated than one-dimensional case.
E-mail address: [email protected]. 1051-2004/$ – see front matter doi:10.1016/j.dsp.2008.12.004
© 2009 Elsevier Inc.
All rights reserved.
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
661
Some of approaches used in one-dimensional digital filter design may straightforward be extended to two-dimensional case; whereas most familiar procedures of designing one-dimensional digital filter do not readily applied to the two-dimensional case. Because of these factors, the aim of this paper is to propose a simple method considered in the frequency domain to solve the design problem of a two-dimensional fractional-order digital differentiator, which utilizes an evolutionary computation called the differential evolution (DE). The DE is a population-based search algorithm and directly uses real-valued manipulations to obtain numerical optimal solutions. We first give a desired frequency response prototype for two-dimensional fractional-order differentiator and a finite impulse response (FIR) digital filter with two dimension structure is employed and designed using the DE algorithm such that its magnitude response may satisfy the desired fractional-order one. This completes the two-dimensional fractional-order digital differentiator design. The rest of this paper can be organized as follows. Section 2 gives a simple problem formulation. Section 3 describes the detailed DE algorithm and presents a DE-based design step for two-dimensional fractional-order digital differentiator. Section 4 illustrates the efficiency of the proposed method with two simulation examples. Finally, a brief conclusion is given in Section 5. 2. Problem formulation In the time-domain analysis, there are two commonly used approximation methods for deriving the fractional-order derivative. One is the Grünwald–Letnikov method and the other is the Riemann–Liouville method. These two approximations are equivalent for a wide class of functions [5]. To perform the fractional-order derivative on a given function f (t ), the definition of the Grünwald–Letnikov method is given by [8] ∞ 1
D α f (t ) = lim
τ →0 t =N τ
α k
=
τα
(−1)k
k=0
α k
f (t − kτ ) and
Γ (α + 1) , Γ (k + 1)Γ (α − k + 1)
(1a)
(1b)
where α ∈ is the arbitrary order, τ is the time increment, Γ (x) is the Gamma function and Γ (k + 1) = k! for a positive integer k. Furthermore, by this definition and assuming f (t ) = 0 for t < 0, the fractional derivative D α f (t ) evaluated at t = nτ can be calculated by the following expression: D α f (nτ ) =
n 1
τα
(−1)k
k=0
α k
f (nτ − kτ ).
(2)
From a signal processing point of view, this definition seems to be intuitive and is considerably suitable for a discrete-time implementation [5]. Another approximation for deriving a fractional derivative is the Riemann–Liouville method defined by D α f (t ) =
1
dn
Γ (n − α ) dt n
t a
f (h) dh, (t − h)α +1−n
α > 0,
(3)
where a is the lower limit of integral and n − 1 < α < n. Recently, in the literature [3] the author developed an alternative expression which is based on using Taylor series expansion for handling the fractional-order derivative. It is defined by the following expression: D α f (t ) =
∞ Γ (n + 1) n−α 1 (n) f (0) t . n! Γ (n − α + 1)
(4)
n=0
The definitions above are the fundamental of designing digital differentiators and they all are considered in the time domain. On the other hand, a two-dimensional FIR digital filter structure can be generally described by y (n1 , n2 ) =
K2 K1
h(k1 , k2 )x(n1 − k1 , n2 − k2 ),
(5)
k1 =0 k2 =0
where x denotes a two-dimensional digital input signal and h represents a unit impulse response of two-dimensional filter, y is the output of the filter, K 1 and K 2 are the numbers of past inputs required, respectively. Eq. (5) is also referred to as the two-dimensional convolution sum. If the input signal x is assumed to be a complex sinusoid x(n1 , n2 ) = exp( jn1 Ω1 + jn2 Ω2 ),
(6)
where Ω1 and Ω2 are called the horizontal and vertical digital frequencies, respectively, then the corresponding output signal y can be expressed by
662
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
y (n1 , n2 ) =
K2 K1
h(k1 , k2 ) exp j (n1 − k1 )Ω1 + j (n2 − k2 )Ω2
k1 =0 k2 =0
= exp( jn1 Ω1 + jn2 Ω2 )
K2 K1
h(k1 , k2 ) exp(− jk1 Ω1 − jk2 Ω2 )
k1 =0 k2 =0
= x(n1 , n2 ) H (Ω1 , Ω2 ), where
K2 K1
H (Ω1 , Ω2 ) ≡
(7)
h(k1 , k2 ) exp(− jk1 Ω1 − jk2 Ω2 )
(8)
k1 =0 k2 =0
is the frequency response of two-dimensional FIR digital filter. In order for application to the DE algorithm, let
Θ = [θ1 , θ2 , . . . , θ( K 1 +1)×( K 2 +1) ] = h(0, 0), . . . , h(0, K 2 ), h(1, 0), . . . , h(1, K 2 ), . . . , h( K 1 , 0), . . . , h( K 1 , K 2 )
(9)
be a manipulated parameter vector with total number of ( K 1 + 1) × ( K 2 + 1). In other words, the impulse response h(k1 , k2 ) is rearranged from two dimension to one dimension. This defined vector Θ is utilized in the DE algorithm. Traditionally, the notation of partial derivative for a function with two variables, f (t 1 , t 2 ), is expressed by
∂ n1 ∂ n2 f (t 1 , t 2 ), ∂ t n11 ∂ t n22
(10)
where n1 and n2 are integers. However, in this study we will explore the fractional-order partial derivative of two variables. Expanding the integer order in Eq. (10), therefore, to the fractional order yields
∂ α 1 ∂ ε2 f (t 1 , t 2 ), ∂ t 1α1 ∂ t 2α2
(11)
where α1 and α2 are arbitrary orders. In fact, numerical analyses and approximations considered in the time domain for solving this kind of derivative are much complicated and also seem to be rarely studied and developed. Thus this paper is to propose an alternative frequency method to achieve the design of fractional-order digital differentiator of two variables using two-dimensional FIR digital filter structure. The frequency prototype of a desired fractional-order differentiator of two variables can be characterized by D (Ω1 , Ω2 ) = ( j Ω1 )α1 ( j Ω2 )α2 .
(12)
The frequency response H (Ω1 , Ω2 ) of Eq. (8) must be designed to meet the differentiation characteristic. The fractional orders of α1 and α2 are always given by the designer for certain actual specifications. In order to accomplish this goal, the DE algorithm is used and it is clearly introduced in the next section. 3. Differential evolution and design steps In 1997, Storn and Price initially presented the differential evolution (DE) algorithm [13] and it has been proven to be an efficient means in solving engineering optimization problems [14–17]. This algorithm is somewhat similar to the genetic algorithm (GA) because it also contains three evolutionary operations: mutation, crossover, and selection. To tackle the two-dimensional digital differentiator design problem, the cost function for evaluating the performance of a manipulated parameter vector as in Eq. (9) is first defined by
E=
D (Ω1 , Ω2 ) − H (Ω1 , Ω2 ) 2 dΩ2 dΩ1 ,
(13)
Ω1 Ω2
where D (Ω1 , Ω2 ) is the required fractional-order differentiator prototype of two variables as in Eq. (12) and H (Ω1 , Ω2 ) is the frequency response of a two-dimensional FIR digital filter described by Eq. (8). This function significantly guides the algorithm converging to the optimal solution. Minimizing Eq. (13) means that the designed filter magnitude response may more meet the required derivative specification. 3.1. Mutation operation Mutation operation begins with picking three parameter vectors Θα , Θβ , and Θγ at random from the population. These three vectors are used to generate a new parameter vector called the donor vector V = [ v 1 , v 2 , . . . , v (k1 +1)×( K 2 +1) ] by the following formula: V = Θα + F · (Θβ − Θγ ),
(14)
where F is called the mutation constant factor. This resulting donor vector V will further make a crossover with the target vector Θ = [θ1 , θ2 , . . . , θ( K 1 +1)×( K 2 +1) ].
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
663
3.2. Crossover operation Crossover operation simply means exchanging some components between vectors V and Θ . To achieve that, it is required to generate a set of ( K 1 + 1) × ( K 2 + 1) random numbers as {r1 , r2 , . . . , r( K 1 +1)×( K 2 +1) }, where r i is selected from the interval (0, 1), and another set of binary sequences {b1 , b2 , . . . , b( K 1 +1)×( K 2 +1) } is derived by
bi =
1 if r i < C 0
otherwise
for i = 1, 2, . . . , ( K 1 + 1) × ( K 2 + 1),
(15)
where C ∈ (0, 1) represents the crossover rate and is always set to 0.5. Then a trial vector W = [ w 1 , w 2 , . . . , w ( K 1 +1)×( K 2 +1) ] is obtained by
wi =
θi
if b i = 1
vi
if b i = 0
for i = 1, 2, . . . , ( K 1 + 1) × ( K 2 + 1).
(16)
This completes the crossover operation. 3.3. Selection operation Having determined the trial vector W , a selection operation is further performed. The cost function of the resulting trial vector is evaluated and then compared with that of the original target vector. Once the cost function of trial vector is less than that of target vector (i.e., the trial vector is better than target vector), then the trial vector take the place of the target vector; otherwise, this trial vector must be rejected and the target vector still survives in the next generation. 3.4. Stopping criterion In the DE algorithm, a complete execution of the mutation, crossover and selection operations on each parameter vector is referred to as “one generation”. The DE algorithm is terminated when the pre-specified number of generations G has been attained. The steps for two-dimensional fractional-order FIR digital differentiator design using the DE algorithm are listed below. Data: Filter parameters K 1 and K 2 in Eq. (5), frequency response prototype of a desired fractional-order differentiator of two variables in Eq. (12), population size P , mutation constant factor F in Eq. (14), crossover rate C in Eq. (15), and number of generations G. Goal: Minimize the defined cost function in Eq. (13) using the DE algorithm to obtain the two-dimensional FIR digital filter. Step 1. Generate N parameter vectors to form an initial population from the interval [−1, 1] at random. Step 2. If the pre-specified number of generations G is attained, then stop the algorithm. for i = 1 to P Evaluate the cost function of target vector Θi using Eq. (13). Obtain the donor vector V using Eq. (14). Derive a set of binary sequences {b1 , b2 , . . . , b( K 1 +1)×( K 2 +1) } using Eq. (15). Apply the crossover formula of Eq. (16) to obtain the trial vector W . Evaluate the cost function of trial vector W also using Eq. (13). Perform the selection operation on both W and Θi to generation a new offspring Θinew . If E ( W ) < E (Θi ), then Θinew = W , otherwise Θinew = Θi . End for i = 1 to P Θi = Θinew . End Step 3. Go back to Step 2. 4. Illustrative examples We here illustrate the efficiency of the proposed method with two different fractional-order examples. The related values assigned to variables of the DE algorithm are listed in Table 1. Besides, the values of K 1 and K 2 are about the orders of two-dimensional FIR digital filter and are always given by the designer. In general, if we require the better design performance, the larger values for K 1 and K 2 are also needed. Here we choose K 1 = 4 and K 2 = 4 for simulations. The DE-based design algorithm is implemented by using the software of Borland C++ 5.02 under PC environments. In this study, the first example is to set fractional orders being α1 = 0.5 and α2 = 0.5 for the desired frequency response specification, i.e., D (Ω1 , Ω2 ) = ( j Ω1 )0.5 ( j Ω2 )0.5 with 0 Ω1 , Ω2 3. The second example is to choose the fractional orders as α1 = 1.2 and α2 = 1.2 for simulations. Simulation results are given in Figs. 1–6 and in Tables 2 and 3, respectively. For the first example, the desired fractional-order magnitude response of two variables is displayed in Fig. 1 and the magnitude response of
664
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
Table 1 The related values assigned to variables of the DE algorithm. Population size Mutation constant factor Crossover rate Number of generations
P = 60 F = 0.05 C = 0.5 G = 1000
Fig. 1. The desired fractional-order magnitude response of two variables, | D (Ω1 , Ω2 )| = |( j Ω1 )0.5 ( j Ω2 )0.5 |.
Fig. 2. The designed magnitude response of two-dimensional FIR digital differentiator, | H (Ω1 , Ω2 )|, for the first example.
two-dimensional FIR digital differentiator designed by the proposed method is then shown in Fig. 2. It is clearly seen that a better approximation can be achieved. Fig. 3 further demonstrates the convergence of the optimal cost function with respect to number of generations to verify the utility of the DE algorithm. Also, a satisfactory convergence is obtained. Table 2 lists the resulting impulse responses of the two-dimensional FIR digital filter and these values are always located in the interval [−1, 1]. Similarly, Figs. 4–6 and Table 3 then show the relative results to the second example. 5. Conclusions In this paper, we have successfully presented a novel design method for designing the fractional-order two-dimensional FIR digital differentiator. The impulse responses of two-dimensional FIR digital filter are solved by using the DE algorithm such that its corresponding magnitude response may satisfy that of the desired fractional-order differentiator of two variables. This scheme is considered in the frequency domain. Finally, two illustrative examples with different fractional orders
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
665
Fig. 3. Convergence trajectory of cost function for the first example.
Table 2 The designed impulse responses of two-dimensional FIR digital filter for the first example. h(k1 , k2 ) k1 k1 k1 k1 k1
=0 =1 =2 =3 =4
k2 = 0
k2 = 1
k2 = 2
k2 = 3
k2 = 4
0.057463 0.277533 −0.266738 −0.121528 0.074812
−0.080047 0.617100 0.048996 −0.428500 −0.186711
−0.273819 −0.165218 −0.216396 0.527541 0.082124
0.214724 −0.742744 0.776760 −0.145418 −0.082790
0.068500 0.045818 −0.375454 0.2.4639 0.048852
Table 3 The designed impulse responses of two-dimensional FIR digital filter for the second example. h(k1 , k2 ) k1 k1 k1 k1 k1
=0 =1 =2 =3 =4
k2 = 0
k2 = 1
k2 = 2
k2 = 3
k2 = 4
−0.549183 0.397936 −0.051481 −0.072136 0.239871
0.884649 −1.758188 1.220823 −0.963963 0.453469
−0.626167 2.165458 −1.542556 0.490688 −0.490421
0.411015 −1.416074 0.115456 0.899382 0.012883
−0.188830 0.614723 0.209984 −0.423311 −0.191813
Fig. 4. The desired fractional-order magnitude response of two variables, | D (Ω1 , Ω2 )| = |( j Ω1 )1.2 ( j Ω2 )1.2 |.
666
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
Fig. 5. The designed magnitude response of two-dimensional FIR digital differentiator, | H (Ω1 , Ω2 )|, for the second example.
Fig. 6. Convergence trajectory of cost function for the second example.
are provided to examine the validity of the proposed method. It can be concluded from simulation results that the DE-based method is considerably valid and efficient in designing the two-dimensional fractional-order FIR digital differentiator. Acknowledgment This work was partially supported by the National Science Council of Taiwan under Grant NSC 96-2221-E-366-001. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
J.J. Shyu, S.C. Pei, A generalized approach to the design of variable fractional-order FIR digital filters, Signal Process. 88 (2008) 1428–1435. C.C. Tseng, Design of variable and adaptive fractional order FIR differentiators, Signal Process. 86 (2006) 2554–2566. C.C. Tseng, Design of fractional order digital FIR differentiators, IEEE Signal Process. Lett. 8 (2001) 77–79. S. Samadi, M. Omair Ahmad, M.N.S. Swamy, Exact fractional-order differentiators for polynomial signals, IEEE Signal Process. Lett. 11 (2004) 529–532. R.S. Barbosa, J.A. Tenreiro Machado, M.F. Silva, Time domain design of fractional differintegrators using least-squares, Signal Process. 86 (2006) 2567– 2581. W.M. Ahmad, Hyperchaos in fractional order nonlinear systems, Chaos Solitons Fractals 26 (2005) 1459–1465. I. Podlubny, Fractional-order systems and P I λ D μ -controllers, IEEE Trans. Automat. Control 44 (1999) 208–214. C. Hwang, J.F. Leu, S.Y. Tsay, A note on time-domain simulation of feedback fractional-order systems, IEEE Trans. Automat. Control 47 (2002) 625–631. H. Samavati, A. Hajimiri, A. Shahani, G. Nasserbakht, T. Lee, Fractal capacitors, IEEE J. Solid-State Circ. 33 (1998) 2035–2041. D.E. Dudgeon, R.M. Mersereau, Multidimensional Digital Signal Processing, Prentice Hall, New Jersey, 1984. V. Singh, New LMI condition for the nonexistence of overflow oscillations in 2-D state-space digital filters using saturation arithmetic, Digit. Signal Process. 17 (2007) 345–352. J.A. Cadzow, Multi-dimensional recursive linear system synthesis, Digit. Signal Process. 8 (1998) 139–157.
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
667
[13] R. Storn, K. Price, Differential evolution – A simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim. 11 (1997) 341–359. [14] M. Varadarajan, K.S. Swarup, Network loss minimization with voltage security using differential evolution, Electr. Power Syst. Res. 78 (2008) 815–823. [15] B. Liu, L. Wang, Y.H. Jin, D.X. Huang, F. Tang, Control and synchronization of chaotic systems by differential evolution algorithm, Chaos Solitons Fractals 34 (2007) 412–419. [16] R. Angira, A. Santosh, Optimization of dynamic systems: A trigonometric differential evolution approach, Comput. Chem. Eng. 31 (2007) 1055–1063. [17] W.D. Chang, Parameter identification of Chen and Lü systems: A differential evolution approach, Chaos Solitons Fractals 32 (2007) 1464–1471.
Wei-Der Chang was born in Kaohsiung, Taiwan, in 1970. He received the Ph.D. degree in electrical engineering from National Sun Yat-Sen University, Kaohsiung, Taiwan, in 2002. He is now an associate professor at the Department of Computer and Communication, Shu-Te University, Kaohsiung, Taiwan. His research interests are in digital signal processing, digital filter design, evolutionary computations, chaotic secure communication, and control engineering.
Contents lists available at ScienceDirect
Digital Signal Processing www.elsevier.com/locate/dsp
Two-dimensional fractional-order digital differentiator design by using differential evolution algorithm Wei-Der Chang Department of Computer and Communication, Shu-Te University, Kaohsiung 824, Taiwan
a r t i c l e
i n f o
a b s t r a c t
Article history: Available online 20 January 2009 Keywords: Two-dimensional digital filter Fractional-order differentiator Differential evolution
Designing a fractional-order digital differentiator often requires considerably complex mathematical operations and numerical approximations. Thus this paper will propose a simple method to achieve the fractional-order digital differentiator design, particularly for two-dimensional fractional-order differentiators. A two-dimensional finite impulse response (FIR) digital filter structure is utilized and designed so that its corresponding magnitude response can satisfy that of a desired fractional-order differentiator of two variables. The algorithm used to design such two-dimensional digital differentiator is the differential evolution (DE), which is one of evolutionary computations and has excellent searching capacity. The efficiency of the proposed scheme can be confirmed by some illustrative examples. © 2009 Elsevier Inc. All rights reserved.
1. Introduction In recent years, the research topic regarding fractional calculus has attracted much attention and successfully introduced to a variety of engineering applications, such as digital signal processing [1–5], chaos systems [6], control system design [7,8], and electrical circuit devices [9]. Generally speaking, the notion of fractional calculus is to enlarge integer order to fractional order in numerical representations. In this paper, we will focus our effort on the issue of fractional-order digital differentiator design, especially for multiple dimensions. About this issue, researchers have presented various design schemes in designing the fractional-order digital differentiator, including time-domain Taylor series expansions [3], Newton series expansions [4], least-squares method in the time domain [5], etc. Basically, the fractional-order digital differentiator is an extended version of general digital differentiator, which generalizes integer order to fractional case, to provide a more flexibility in the real applications. Its role is to give a fractional-order differentiation on a given digital input signal. In comparison with designing the general digital differentiator, designing such a fractional-order case is considerably complicated and difficult. On the other hand, another topic discussed in this paper is about the two-dimensional digital filter design [10–12]. Twodimensional digital filter is to tackle the two-dimensional digital signal, for example, digital image signal. It is obvious that efforts in designing a two-dimensional digital filter are much beyond that in designing a one-dimensional case. The reasons are as follows [10]: 1. Two-dimensional digital filter basically involves more data than one-dimensional case. 2. Mathematics for tackling the two-dimensional digital filter is more incomplete than that for tackling the onedimensional case. 3. In the design, two-dimensional digital filter has more degrees of freedom to be evaluated than one-dimensional case.
E-mail address: [email protected]. 1051-2004/$ – see front matter doi:10.1016/j.dsp.2008.12.004
© 2009 Elsevier Inc.
All rights reserved.
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
661
Some of approaches used in one-dimensional digital filter design may straightforward be extended to two-dimensional case; whereas most familiar procedures of designing one-dimensional digital filter do not readily applied to the two-dimensional case. Because of these factors, the aim of this paper is to propose a simple method considered in the frequency domain to solve the design problem of a two-dimensional fractional-order digital differentiator, which utilizes an evolutionary computation called the differential evolution (DE). The DE is a population-based search algorithm and directly uses real-valued manipulations to obtain numerical optimal solutions. We first give a desired frequency response prototype for two-dimensional fractional-order differentiator and a finite impulse response (FIR) digital filter with two dimension structure is employed and designed using the DE algorithm such that its magnitude response may satisfy the desired fractional-order one. This completes the two-dimensional fractional-order digital differentiator design. The rest of this paper can be organized as follows. Section 2 gives a simple problem formulation. Section 3 describes the detailed DE algorithm and presents a DE-based design step for two-dimensional fractional-order digital differentiator. Section 4 illustrates the efficiency of the proposed method with two simulation examples. Finally, a brief conclusion is given in Section 5. 2. Problem formulation In the time-domain analysis, there are two commonly used approximation methods for deriving the fractional-order derivative. One is the Grünwald–Letnikov method and the other is the Riemann–Liouville method. These two approximations are equivalent for a wide class of functions [5]. To perform the fractional-order derivative on a given function f (t ), the definition of the Grünwald–Letnikov method is given by [8] ∞ 1
D α f (t ) = lim
τ →0 t =N τ
α k
=
τα
(−1)k
k=0
α k
f (t − kτ ) and
Γ (α + 1) , Γ (k + 1)Γ (α − k + 1)
(1a)
(1b)
where α ∈ is the arbitrary order, τ is the time increment, Γ (x) is the Gamma function and Γ (k + 1) = k! for a positive integer k. Furthermore, by this definition and assuming f (t ) = 0 for t < 0, the fractional derivative D α f (t ) evaluated at t = nτ can be calculated by the following expression: D α f (nτ ) =
n 1
τα
(−1)k
k=0
α k
f (nτ − kτ ).
(2)
From a signal processing point of view, this definition seems to be intuitive and is considerably suitable for a discrete-time implementation [5]. Another approximation for deriving a fractional derivative is the Riemann–Liouville method defined by D α f (t ) =
1
dn
Γ (n − α ) dt n
t a
f (h) dh, (t − h)α +1−n
α > 0,
(3)
where a is the lower limit of integral and n − 1 < α < n. Recently, in the literature [3] the author developed an alternative expression which is based on using Taylor series expansion for handling the fractional-order derivative. It is defined by the following expression: D α f (t ) =
∞ Γ (n + 1) n−α 1 (n) f (0) t . n! Γ (n − α + 1)
(4)
n=0
The definitions above are the fundamental of designing digital differentiators and they all are considered in the time domain. On the other hand, a two-dimensional FIR digital filter structure can be generally described by y (n1 , n2 ) =
K2 K1
h(k1 , k2 )x(n1 − k1 , n2 − k2 ),
(5)
k1 =0 k2 =0
where x denotes a two-dimensional digital input signal and h represents a unit impulse response of two-dimensional filter, y is the output of the filter, K 1 and K 2 are the numbers of past inputs required, respectively. Eq. (5) is also referred to as the two-dimensional convolution sum. If the input signal x is assumed to be a complex sinusoid x(n1 , n2 ) = exp( jn1 Ω1 + jn2 Ω2 ),
(6)
where Ω1 and Ω2 are called the horizontal and vertical digital frequencies, respectively, then the corresponding output signal y can be expressed by
662
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
y (n1 , n2 ) =
K2 K1
h(k1 , k2 ) exp j (n1 − k1 )Ω1 + j (n2 − k2 )Ω2
k1 =0 k2 =0
= exp( jn1 Ω1 + jn2 Ω2 )
K2 K1
h(k1 , k2 ) exp(− jk1 Ω1 − jk2 Ω2 )
k1 =0 k2 =0
= x(n1 , n2 ) H (Ω1 , Ω2 ), where
K2 K1
H (Ω1 , Ω2 ) ≡
(7)
h(k1 , k2 ) exp(− jk1 Ω1 − jk2 Ω2 )
(8)
k1 =0 k2 =0
is the frequency response of two-dimensional FIR digital filter. In order for application to the DE algorithm, let
Θ = [θ1 , θ2 , . . . , θ( K 1 +1)×( K 2 +1) ] = h(0, 0), . . . , h(0, K 2 ), h(1, 0), . . . , h(1, K 2 ), . . . , h( K 1 , 0), . . . , h( K 1 , K 2 )
(9)
be a manipulated parameter vector with total number of ( K 1 + 1) × ( K 2 + 1). In other words, the impulse response h(k1 , k2 ) is rearranged from two dimension to one dimension. This defined vector Θ is utilized in the DE algorithm. Traditionally, the notation of partial derivative for a function with two variables, f (t 1 , t 2 ), is expressed by
∂ n1 ∂ n2 f (t 1 , t 2 ), ∂ t n11 ∂ t n22
(10)
where n1 and n2 are integers. However, in this study we will explore the fractional-order partial derivative of two variables. Expanding the integer order in Eq. (10), therefore, to the fractional order yields
∂ α 1 ∂ ε2 f (t 1 , t 2 ), ∂ t 1α1 ∂ t 2α2
(11)
where α1 and α2 are arbitrary orders. In fact, numerical analyses and approximations considered in the time domain for solving this kind of derivative are much complicated and also seem to be rarely studied and developed. Thus this paper is to propose an alternative frequency method to achieve the design of fractional-order digital differentiator of two variables using two-dimensional FIR digital filter structure. The frequency prototype of a desired fractional-order differentiator of two variables can be characterized by D (Ω1 , Ω2 ) = ( j Ω1 )α1 ( j Ω2 )α2 .
(12)
The frequency response H (Ω1 , Ω2 ) of Eq. (8) must be designed to meet the differentiation characteristic. The fractional orders of α1 and α2 are always given by the designer for certain actual specifications. In order to accomplish this goal, the DE algorithm is used and it is clearly introduced in the next section. 3. Differential evolution and design steps In 1997, Storn and Price initially presented the differential evolution (DE) algorithm [13] and it has been proven to be an efficient means in solving engineering optimization problems [14–17]. This algorithm is somewhat similar to the genetic algorithm (GA) because it also contains three evolutionary operations: mutation, crossover, and selection. To tackle the two-dimensional digital differentiator design problem, the cost function for evaluating the performance of a manipulated parameter vector as in Eq. (9) is first defined by
E=
D (Ω1 , Ω2 ) − H (Ω1 , Ω2 ) 2 dΩ2 dΩ1 ,
(13)
Ω1 Ω2
where D (Ω1 , Ω2 ) is the required fractional-order differentiator prototype of two variables as in Eq. (12) and H (Ω1 , Ω2 ) is the frequency response of a two-dimensional FIR digital filter described by Eq. (8). This function significantly guides the algorithm converging to the optimal solution. Minimizing Eq. (13) means that the designed filter magnitude response may more meet the required derivative specification. 3.1. Mutation operation Mutation operation begins with picking three parameter vectors Θα , Θβ , and Θγ at random from the population. These three vectors are used to generate a new parameter vector called the donor vector V = [ v 1 , v 2 , . . . , v (k1 +1)×( K 2 +1) ] by the following formula: V = Θα + F · (Θβ − Θγ ),
(14)
where F is called the mutation constant factor. This resulting donor vector V will further make a crossover with the target vector Θ = [θ1 , θ2 , . . . , θ( K 1 +1)×( K 2 +1) ].
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
663
3.2. Crossover operation Crossover operation simply means exchanging some components between vectors V and Θ . To achieve that, it is required to generate a set of ( K 1 + 1) × ( K 2 + 1) random numbers as {r1 , r2 , . . . , r( K 1 +1)×( K 2 +1) }, where r i is selected from the interval (0, 1), and another set of binary sequences {b1 , b2 , . . . , b( K 1 +1)×( K 2 +1) } is derived by
bi =
1 if r i < C 0
otherwise
for i = 1, 2, . . . , ( K 1 + 1) × ( K 2 + 1),
(15)
where C ∈ (0, 1) represents the crossover rate and is always set to 0.5. Then a trial vector W = [ w 1 , w 2 , . . . , w ( K 1 +1)×( K 2 +1) ] is obtained by
wi =
θi
if b i = 1
vi
if b i = 0
for i = 1, 2, . . . , ( K 1 + 1) × ( K 2 + 1).
(16)
This completes the crossover operation. 3.3. Selection operation Having determined the trial vector W , a selection operation is further performed. The cost function of the resulting trial vector is evaluated and then compared with that of the original target vector. Once the cost function of trial vector is less than that of target vector (i.e., the trial vector is better than target vector), then the trial vector take the place of the target vector; otherwise, this trial vector must be rejected and the target vector still survives in the next generation. 3.4. Stopping criterion In the DE algorithm, a complete execution of the mutation, crossover and selection operations on each parameter vector is referred to as “one generation”. The DE algorithm is terminated when the pre-specified number of generations G has been attained. The steps for two-dimensional fractional-order FIR digital differentiator design using the DE algorithm are listed below. Data: Filter parameters K 1 and K 2 in Eq. (5), frequency response prototype of a desired fractional-order differentiator of two variables in Eq. (12), population size P , mutation constant factor F in Eq. (14), crossover rate C in Eq. (15), and number of generations G. Goal: Minimize the defined cost function in Eq. (13) using the DE algorithm to obtain the two-dimensional FIR digital filter. Step 1. Generate N parameter vectors to form an initial population from the interval [−1, 1] at random. Step 2. If the pre-specified number of generations G is attained, then stop the algorithm. for i = 1 to P Evaluate the cost function of target vector Θi using Eq. (13). Obtain the donor vector V using Eq. (14). Derive a set of binary sequences {b1 , b2 , . . . , b( K 1 +1)×( K 2 +1) } using Eq. (15). Apply the crossover formula of Eq. (16) to obtain the trial vector W . Evaluate the cost function of trial vector W also using Eq. (13). Perform the selection operation on both W and Θi to generation a new offspring Θinew . If E ( W ) < E (Θi ), then Θinew = W , otherwise Θinew = Θi . End for i = 1 to P Θi = Θinew . End Step 3. Go back to Step 2. 4. Illustrative examples We here illustrate the efficiency of the proposed method with two different fractional-order examples. The related values assigned to variables of the DE algorithm are listed in Table 1. Besides, the values of K 1 and K 2 are about the orders of two-dimensional FIR digital filter and are always given by the designer. In general, if we require the better design performance, the larger values for K 1 and K 2 are also needed. Here we choose K 1 = 4 and K 2 = 4 for simulations. The DE-based design algorithm is implemented by using the software of Borland C++ 5.02 under PC environments. In this study, the first example is to set fractional orders being α1 = 0.5 and α2 = 0.5 for the desired frequency response specification, i.e., D (Ω1 , Ω2 ) = ( j Ω1 )0.5 ( j Ω2 )0.5 with 0 Ω1 , Ω2 3. The second example is to choose the fractional orders as α1 = 1.2 and α2 = 1.2 for simulations. Simulation results are given in Figs. 1–6 and in Tables 2 and 3, respectively. For the first example, the desired fractional-order magnitude response of two variables is displayed in Fig. 1 and the magnitude response of
664
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
Table 1 The related values assigned to variables of the DE algorithm. Population size Mutation constant factor Crossover rate Number of generations
P = 60 F = 0.05 C = 0.5 G = 1000
Fig. 1. The desired fractional-order magnitude response of two variables, | D (Ω1 , Ω2 )| = |( j Ω1 )0.5 ( j Ω2 )0.5 |.
Fig. 2. The designed magnitude response of two-dimensional FIR digital differentiator, | H (Ω1 , Ω2 )|, for the first example.
two-dimensional FIR digital differentiator designed by the proposed method is then shown in Fig. 2. It is clearly seen that a better approximation can be achieved. Fig. 3 further demonstrates the convergence of the optimal cost function with respect to number of generations to verify the utility of the DE algorithm. Also, a satisfactory convergence is obtained. Table 2 lists the resulting impulse responses of the two-dimensional FIR digital filter and these values are always located in the interval [−1, 1]. Similarly, Figs. 4–6 and Table 3 then show the relative results to the second example. 5. Conclusions In this paper, we have successfully presented a novel design method for designing the fractional-order two-dimensional FIR digital differentiator. The impulse responses of two-dimensional FIR digital filter are solved by using the DE algorithm such that its corresponding magnitude response may satisfy that of the desired fractional-order differentiator of two variables. This scheme is considered in the frequency domain. Finally, two illustrative examples with different fractional orders
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
665
Fig. 3. Convergence trajectory of cost function for the first example.
Table 2 The designed impulse responses of two-dimensional FIR digital filter for the first example. h(k1 , k2 ) k1 k1 k1 k1 k1
=0 =1 =2 =3 =4
k2 = 0
k2 = 1
k2 = 2
k2 = 3
k2 = 4
0.057463 0.277533 −0.266738 −0.121528 0.074812
−0.080047 0.617100 0.048996 −0.428500 −0.186711
−0.273819 −0.165218 −0.216396 0.527541 0.082124
0.214724 −0.742744 0.776760 −0.145418 −0.082790
0.068500 0.045818 −0.375454 0.2.4639 0.048852
Table 3 The designed impulse responses of two-dimensional FIR digital filter for the second example. h(k1 , k2 ) k1 k1 k1 k1 k1
=0 =1 =2 =3 =4
k2 = 0
k2 = 1
k2 = 2
k2 = 3
k2 = 4
−0.549183 0.397936 −0.051481 −0.072136 0.239871
0.884649 −1.758188 1.220823 −0.963963 0.453469
−0.626167 2.165458 −1.542556 0.490688 −0.490421
0.411015 −1.416074 0.115456 0.899382 0.012883
−0.188830 0.614723 0.209984 −0.423311 −0.191813
Fig. 4. The desired fractional-order magnitude response of two variables, | D (Ω1 , Ω2 )| = |( j Ω1 )1.2 ( j Ω2 )1.2 |.
666
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
Fig. 5. The designed magnitude response of two-dimensional FIR digital differentiator, | H (Ω1 , Ω2 )|, for the second example.
Fig. 6. Convergence trajectory of cost function for the second example.
are provided to examine the validity of the proposed method. It can be concluded from simulation results that the DE-based method is considerably valid and efficient in designing the two-dimensional fractional-order FIR digital differentiator. Acknowledgment This work was partially supported by the National Science Council of Taiwan under Grant NSC 96-2221-E-366-001. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
J.J. Shyu, S.C. Pei, A generalized approach to the design of variable fractional-order FIR digital filters, Signal Process. 88 (2008) 1428–1435. C.C. Tseng, Design of variable and adaptive fractional order FIR differentiators, Signal Process. 86 (2006) 2554–2566. C.C. Tseng, Design of fractional order digital FIR differentiators, IEEE Signal Process. Lett. 8 (2001) 77–79. S. Samadi, M. Omair Ahmad, M.N.S. Swamy, Exact fractional-order differentiators for polynomial signals, IEEE Signal Process. Lett. 11 (2004) 529–532. R.S. Barbosa, J.A. Tenreiro Machado, M.F. Silva, Time domain design of fractional differintegrators using least-squares, Signal Process. 86 (2006) 2567– 2581. W.M. Ahmad, Hyperchaos in fractional order nonlinear systems, Chaos Solitons Fractals 26 (2005) 1459–1465. I. Podlubny, Fractional-order systems and P I λ D μ -controllers, IEEE Trans. Automat. Control 44 (1999) 208–214. C. Hwang, J.F. Leu, S.Y. Tsay, A note on time-domain simulation of feedback fractional-order systems, IEEE Trans. Automat. Control 47 (2002) 625–631. H. Samavati, A. Hajimiri, A. Shahani, G. Nasserbakht, T. Lee, Fractal capacitors, IEEE J. Solid-State Circ. 33 (1998) 2035–2041. D.E. Dudgeon, R.M. Mersereau, Multidimensional Digital Signal Processing, Prentice Hall, New Jersey, 1984. V. Singh, New LMI condition for the nonexistence of overflow oscillations in 2-D state-space digital filters using saturation arithmetic, Digit. Signal Process. 17 (2007) 345–352. J.A. Cadzow, Multi-dimensional recursive linear system synthesis, Digit. Signal Process. 8 (1998) 139–157.
W.-D. Chang / Digital Signal Processing 19 (2009) 660–667
667
[13] R. Storn, K. Price, Differential evolution – A simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim. 11 (1997) 341–359. [14] M. Varadarajan, K.S. Swarup, Network loss minimization with voltage security using differential evolution, Electr. Power Syst. Res. 78 (2008) 815–823. [15] B. Liu, L. Wang, Y.H. Jin, D.X. Huang, F. Tang, Control and synchronization of chaotic systems by differential evolution algorithm, Chaos Solitons Fractals 34 (2007) 412–419. [16] R. Angira, A. Santosh, Optimization of dynamic systems: A trigonometric differential evolution approach, Comput. Chem. Eng. 31 (2007) 1055–1063. [17] W.D. Chang, Parameter identification of Chen and Lü systems: A differential evolution approach, Chaos Solitons Fractals 32 (2007) 1464–1471.
Wei-Der Chang was born in Kaohsiung, Taiwan, in 1970. He received the Ph.D. degree in electrical engineering from National Sun Yat-Sen University, Kaohsiung, Taiwan, in 2002. He is now an associate professor at the Department of Computer and Communication, Shu-Te University, Kaohsiung, Taiwan. His research interests are in digital signal processing, digital filter design, evolutionary computations, chaotic secure communication, and control engineering.