# A Computational Fractional Signal Derivative Method

Authors: Matías Salinas; Rodrigo Salas; Diego Mellado; Antonio Glaría; Carolina Saavedra
Journal: Modelling and Simulation in Engineering
 (2018)

---

## Abstract

We propose an efficient computational method to obtain the fractional derivative of a digital signal. The proposal consists of a new interpretation of the Grünwald–Letnikov differintegral operator where we have introduced a finite Cauchy convolution with the Grünwald–Letnikov dynamic kernel. The method can be applied to any signal without knowing its analytical form. In the experiments, we have compared the proposed Grünwald–Letnikov computational fractional derivative method with the Riemman–Louville fractional derivative approach for two well-known functions. The simulations exhibit similar results for both methods; however, the Grünwald–Letnikov method outperforms the other approach in execution time. Finally, we show an application of how our proposal can be useful to find the fractional relationship between two well-known biomedical signals.


---

## Body

## 1. Introduction

Signal processing has been a powerful tool for system analysis when the system’s analytical model is unknown. In order to analyse these kinds of systems, it is necessary to apply high-level transformation operations to understand the signal. Within this group of operations, fractional calculus has become very useful for feature extraction, systems control, and modelling speech and more complex applications such as finding the relationship between different biomedical signals [1–5].The fractional calculus has emerged as a nonlocal theory described with operators of fractional nature [6]. Fractional calculus was born as a natural generalization of the traditional calculus (Leibniz 1695, Euler 1730, Fourier 1822, Abel 1823, among others [6–8]); however, until recently, this mathematical theory is playing an active role in disciplines such as physics and control theory [9].In the last decade, several applications have emerged due to the fractional nature of the phenomena. For instance, in physics, fractional calculus has been applied to thermodynamics, materials, and waves [9, 10]. In a fractional optimal control, either the performance index or the differential equations governing the dynamics of the system contains a term with a fractional derivative [11]. Recently, Tapia et al. [1] and Salinas et al. [12] have applied fractional calculus to discover a fractional relationship between two different biomedical signal modalities.Fractional calculus is a terminology that refers to the integration and differentiation of arbitrary order [6, 8]; in other words, the meaning of k-th derivative dky/dxk and k-th iterated integral ∫…∫⋅dx are extended by considering a fractional α∈ℝ+ parameter instead of the integer k∈ℕ parameter.On the other hand, the numerical computation of the derivatives of the signal has many uses in analytical signal processing. The first derivative can be interpreted as the slope of the tangent to the signal at each point, and the second derivative is a measure of the curvature of the signal (the rate of change of the slope). The derivatives have been used to smoothen the signals by filtering the noise, for peak detection, trace analysis, and feature extraction, among others [13–16].The aim of this article is to propose an efficient computational method of the fractional calculus applications to digital signal processing. Thus, we introduce a numerical formulation of the fractional derivatives for continuous and discrete signals, and also, we introduce the concepts ofCauchy finite convolution and the Grünwald–Letnikov (GL) dynamic kernel. The fractional derivatives and integrals are obtained by computing these convolutions with the GL kernel and the signal of interest. Our proposal can be viewed as a dynamic filter applied to a causal realization of a stochastic process.This article is organized as follows. In Section2, we introduce the fundamentals concepts of fractional calculus. Afterwards, in Section 3, we formulate our proposal of the fractional derivative of a signal as a finite convolution. In Section 4, simulations results are obtained in order to numerically quantify the approximation error. Section 5 shows the mathematical relationship between two biosignals by applying our fractional derivative approach. Concluding remarks and further works are given in the last section.


## 2. Basic Principles of Fractional Calculus

The fractional calculus can be obtained either from the generalization of the definition of the derivative or the definition of the integral. Depending on the type of approach selected, the resulting mathematical equation is different. If the extension is build from the derivative, then the Grünwald–Letkikov method is obtained. On the other hand, the Riemman–Louville (RL) method is achieved when the extension is made from the integral.It is well known that the simple derivative offx with respect to x is the function f′x, defined as(1)f′x=ddxfx=limh→0fx−fx−hh.The derivative of a function represents an infinitesimal change in the function with respect to one of its variables. For higher-order derivatives, letfk be the k-th derivative, and the k-th derivative is defined by taking the derivative of the previous function fk−1 as follows:(2)fkx=ddxfk−1x=limh→0fk−1x−fk−1x−hh.If we apply the derivative recursively from the functionfx, we obtain the following expression:(3)fkx=dkdxkfx=limh→01hk∑m=0k−1mkmfx−mh,where km is the binomial coefficient indexed by a pair of integers, k and m, 0≤m≤k. If the notation k! corresponds to the k-th factorial, then the value of the coefficient is obtained as follows:(4)km=k!m!k−m!.Equation (3) can be adjusted to define the k-th integral as an extension of the derivative, where a negative sign for the k-th derivative is used. As a result, the following expression for the integral is obtained:(5)f−kx=∫⋯∫fxdx=limh→0hk∑m=0k−1m−kmfx−mh.In this case, the binomial coefficient for a negative value is computed as follows:(6)−km=−1mk+m−1!k−1!m!.In Sections2.1 and 2.2, we introduce the definitions of fractional derivatives and integrals given by both the Grünwald–Letnikov derivative approach and the Riemann–Liouville integrative approach.


### 2.1. Grünwald–Letnikov Differintegral Operator

The Grünwald–Letnikov differintegral is a combined differentiation/integration operator. The Grünwald–Letnikovα-differintegral of function f is denoted by Dαf. The GL operator is based on a classical concept for which integer order derivatives can be represented as limits of finite differences (3).The binomial coefficient of (4) can be generalized to real arguments by means of the Euler’s gamma function:(7)Γx=∫0∞tx−1e−tdtℜx>0,where the following properties hold: Γ1=1, Γx+1=xΓx and if k is an integer number, then Γk+1=k!. The previous equation converges for all x∈ℂ,ℜx>0. If we replace the k argument with a real positive parameter α∈ℝ+ in (4), then we obtain the following property:(8)αm=Γα+1m!Γα−m+1.Grünwald–Letnikov introduced (8) in (3), obtaining the following definition for the fractional α-th derivative:(9)Dαfx=limh→01hα∑m=0p−q/h−1mΓα+1m!Γα−m+1fx−mh,where p and q are the upper and lower limits of differentiation, respectively.If we consider that theα parameter can take negative values, then the fractional integrative is defined. For α∈ℝ+, we consider the −α for the GL integral definition to obtain the following expression:(10)D−αfx=limh→0hα∑m=0p−q/hΓα+mm!Γαfx−mh,where the binomial coefficient is expanded as follows:(11)−αm=−1mΓα+mm!Γα.


### 2.2. Riemann–Liouville Differintegral Operator

On the other hand, Riemann and Liouville have extended traditional integral in order to obtain its fractional model. The notation for theα−integral is Jαf, and the α−differintegral being Dαx0f. Thus, the RL fractional integral is defined as(12)Jαft=1Γα∫0tt−τα−1fτdτ.And, the fractional derivative of the RL is given by(13)Dxαx0fx=1Γm−αdmdxm∫x0xx−tm−α−1ftdt,where x is the evaluation point, α is the fractional order, and x0 is the reference point of integration.Several studies have shown important properties and conditions for the existence of the fractional models [17–19]. Although, the analysis of these properties is out of the scope of this article, a detailed discussion of Grünwald–Letnikov differential and integral operators can be found in the monograph of Samko et al. [20].


## 2.1. Grünwald–Letnikov Differintegral Operator

The Grünwald–Letnikov differintegral is a combined differentiation/integration operator. The Grünwald–Letnikovα-differintegral of function f is denoted by Dαf. The GL operator is based on a classical concept for which integer order derivatives can be represented as limits of finite differences (3).The binomial coefficient of (4) can be generalized to real arguments by means of the Euler’s gamma function:(7)Γx=∫0∞tx−1e−tdtℜx>0,where the following properties hold: Γ1=1, Γx+1=xΓx and if k is an integer number, then Γk+1=k!. The previous equation converges for all x∈ℂ,ℜx>0. If we replace the k argument with a real positive parameter α∈ℝ+ in (4), then we obtain the following property:(8)αm=Γα+1m!Γα−m+1.Grünwald–Letnikov introduced (8) in (3), obtaining the following definition for the fractional α-th derivative:(9)Dαfx=limh→01hα∑m=0p−q/h−1mΓα+1m!Γα−m+1fx−mh,where p and q are the upper and lower limits of differentiation, respectively.If we consider that theα parameter can take negative values, then the fractional integrative is defined. For α∈ℝ+, we consider the −α for the GL integral definition to obtain the following expression:(10)D−αfx=limh→0hα∑m=0p−q/hΓα+mm!Γαfx−mh,where the binomial coefficient is expanded as follows:(11)−αm=−1mΓα+mm!Γα.


## 2.2. Riemann–Liouville Differintegral Operator

On the other hand, Riemann and Liouville have extended traditional integral in order to obtain its fractional model. The notation for theα−integral is Jαf, and the α−differintegral being Dαx0f. Thus, the RL fractional integral is defined as(12)Jαft=1Γα∫0tt−τα−1fτdτ.And, the fractional derivative of the RL is given by(13)Dxαx0fx=1Γm−αdmdxm∫x0xx−tm−α−1ftdt,where x is the evaluation point, α is the fractional order, and x0 is the reference point of integration.Several studies have shown important properties and conditions for the existence of the fractional models [17–19]. Although, the analysis of these properties is out of the scope of this article, a detailed discussion of Grünwald–Letnikov differential and integral operators can be found in the monograph of Samko et al. [20].


## 3. Fractional Calculus for Signal Processing

In this section, we introduce our proposal of a computational method to obtain theα fractional derivative of the signal. Consider T be a subset of 0,∞. A set of stochastic or random variables Xtt∈T indexed by T is called a stochastic continuous-time process. When T=ℕ, the stochastic discrete-time process Xnn∈T is obtained. The signals xt and xn are defined as a realization of a real-valued continuous or discrete random process, respectively. In this work, we will consider only causal signals, and the value of a signal xt equals to 0 when t<0.On a continuous state, a random variable can take any value between a time intervalt∈T. The fractional derivative for a continuous signal can be expressed as follows:(14)Dαxt=limh→01hα∑m=0p/h−1mΓα+1m!Γα−m+1xt−mh,where mh is the spacing between points on a grid and directly related to the approximation error Ohα. Moreover, p is the upper limit of differentiation, and in this case, q=0 is considered.A continuous signal can be converted into a discrete signal to a sequence of samples by sampling the signal everyTs seconds. The sampled signal xn is related to the continuous signal through the equation t=nTs, where Ts is called the sampling interval or sampling period. The sampling frequency fs is defined as the number of samples taken per second, thus fs=1/Ts. In order to capture all the information from the continuous signal, it is necessary to fulfill the Nyquist sampling theorem [21]. The Nyquist rate is the minimum sampling rate that satisfies the sampling theorem and corresponds to twice the maximum component frequency of the function being sampled.The discrete signalxn, n=0,…,N−1, can be expressed as a vector:(15)x=x0,x1,x2,…,xN−1,where N is the number of sample points. The fractional derivative of a discrete signal xn is given by(16)Dαxn=1hα∑m=0n−1mΓα+1m!Γα−m+1xn−m.In order to optimize the definition above, it is possible to define theα-th Grünwald–Letnikov kernel function as follows:(17)Gαm=Cα−1mm!Γα−m+1,where the value of constant Cα depends only on α, as follows:(18)Cα=Γα+1hα.Let us define theGrünwald–Letnikov (GL) dynamic kernel Gnα of size n as(19)Gnα=Gα0,Gα1,Gα2,…,Gαn−1,where the coefficients Gαm, m=0,…,n−1 are given by (17).Then-th Cauchy finite discrete convolution is defined as follows:(20)f⊛ngn=∑m=0nfmgn−m.The Grünwald–Letnikov fractional derivative can be rewritten by computing the Cauchy finite discrete convolution, given by (20), between the signal xn and the Grünwald–Letnikov (GL) dynamic kernel Gnα (19) as follows:(21)Dαxn=Gnα⊛nxn=∑m=0nGαmxn−m.In Figure1, the values of the coefficients of the kernel Gαm, for m=0,1,2,…,10 and α=0.3,0.5,0.7,0.9, are shown. Notice that, when m increases the kernel rapidly converges to 0; however, simulation results have shown that even the small terms significantly contribute to the computation of the fractional derivative.Figure 1
Grünwald–Letnikov kernel (Gαm) under different alpha values.Notice that, for the integral case, theα-th Grünwald–Letnikov kernel function Gα changes to the following expression:(22)G−αm=C−α−1mΓα+mm!,and the constant Cα changes to(23)C−α=hαΓα.The Grünwald–Letnikov fractional integral can be rewritten by computing the Cauchy finite discrete convolution as follows:(24)D−αxn=Gn−α⊛nxn=∑m=0nG−αmxn−m.To numerically compute the fractional derivative of a signal, we apply the following property to the gamma function of (17):(25)Γα−m+1=Γα−m+2α−m+1,where the following well-known property was used Γα+1=αΓα. We have optimized the numerical computation of the kernel (17) with a recursive equation:(26)Gαm=−Gαm−1α−m+1m.This property will optimize the computing time by recursively computing the gamma function. Moreover, it will avoid numerical problems due to the computation of the factorial term for large numbers.Algorithm1 shows the implementation of the Grünwald–Letnikov fractional derivative as the Cauchy finite discrete convolution between the signal and the dynamic kernel. Line 1 through 4 asks for the input signal x, sampling rate h, number of samples points N, and the user defined α parameter of the derivative, respectively. Line 5 through 7 computes several constant terms: the Γα+1, the first term of the dynamic kernel Gα0, and the constant Cα. Between lines 8–10, the GL dynamic kernel is computed. In line 9, we have been applied the property of (26). Between lines 11–17, we compute the fractional derivative. In line 12, we initialize the first component of the fractional derivative. Between lines 13–15, the Cauchy finite discrete convolution between the signal and the dynamic kernel is computed. Line 17 gives the result of the fractional derivative of a signal.Algorithm 1: Algorithm of the fractional derivative of the signal as the Cauchy finite discrete convolution with the GL dynamic kernel.
x
← Discrete signal
h
← Sampling rate
N
← Number of samples in the signal
α
← user defined value
g
←
Γ
α
+
1
G
0
←
1
/
g
C
α
←
g
/
h
α
for m=1 to N−1 do
Gm=−Gm−1∗α−m+1/m
end for
for n=0 to N−1 do
Dn←0
for m=0 to ndo
Dn=Dn+Cα⋅Gm⋅xn−m
end for
end for
Return the arrayD of the fractional derivativeIt is not simple to give an interpretation or meaning to the fractional derivative in signals without knowing the context; however, the fractional derivatives accurately describes natural complex phenomena that occur in physics, engineering, and biological systems [22]. The fractional derivative gives insight about the functional relationship for modelling complex systems. For instance, the application of a fractional derivative with α∈0,1 is similar to a linear convex combination between the signal and its first derivative; nevertheless, the fractional derivative captures the nonlinear behaviour and long-term memory due to its inherent definition [23].


## 4. Simulation Study

For this study, we implemented two differentiable signals in order to calculate their fractional derivatives with differentα values between 0.1 and 0.9 at steps of 0.1. First, an underdamped harmonic oscillator was implemented, and its equation corresponds to the following expression:(27)xt=e−γtsinωt,with γ being the damping coefficient and ω=2πf, with f being its frequency in Hz. For the experiment, we determined ω=5 Hz and the damping coefficient, γ=1. If we apply the Hermmann rules [6] to (27), we are able to obtain the following α-th fractional derivative equation for this signal:(28)∂αxt=−γαe−γtsinωt+e−γtωαsinωt+π2α.The second signal, named by us as “sine-cosine function,” is modelled with the following equation:(29)xt=sinωtcosωt2+1.1,with ω=2π. In order to make the equation continuous in ℝ, we added to the denominator a constant value of 1.1 as shown in (29).We implemented the fractional derivative of the signal as the Cauchy finite discrete convolution with the GL dynamic kernel on Python 3.6, using Numpy and Scipy libraries in order to process continuous and discrete signals. For the RL approach, we used SymPy’s implementation [24]. It is important to mention that this toolbox is able to compute the fractional derivative only for continuous signals.First most, we simulated different fractional derivatives, up to the first derivative for both the harmonic oscillator and sine-cosine function, as presented on Figure2, and both signals were sampled at 100 Hz. We can see that both methods have the same behaviour for different α in both signals.Figure 2
Fractional derivatives for the underdamped harmonic oscillator using (a) Grünwald–Letnikov and (b) Riemann–Liouville (b) approaches. And fractional derivatives for the “sine-cosine function” using (c) Grünwald–Letnikov and (d) Riemann–Liouville approaches.
(a)
(b)
(c)
(d)The second simulation calculates the root mean-squared error (RMSE) for both signals, using SymPy’s RL approach as ground-truth, at sample frequencies of 10, 50, and 100 Hz. This comparison was done in order to study the effects of the sampling rate when calculating the fractional derivative on discrete signals with both methods. Figure3 shows the increase on the RMSE when doing higher fractional derivatives. The Nyquist frequency is 10 Hz for the oscillator, At this sample rate, the RMSE is the highest for almost all derivatives.Figure 3
Root mean-squared error (RMSE) between fractional derivatives calculated with both GL and RL methods, with different sample rates. (a) RMSE of oscillator signal. (b) RMSE of the sine-cosine signal.
(a)
(b)Table1 shows the values obtained for Figure 3. For α<1, the RMSE between both methods is low, and starts to exponentially increase at higher α values, at all sample rates for both signals.Table 1
Root mean-squared error for different fractional derivatives.
Item
Oscillator
Sine-cosine
Alpha
10 Hz
50 Hz
100 Hz
10 Hz
50 Hz
100 Hz
0.1
0.922
0.882
0.142
5.129
4.355
4.320
0.2
1.298
1.231
0.394
6.641
5.570
5.510
0.3
1.829
1.718
0.813
8.643
7.228
7.127
0.4
2.580
2.399
1.479
11.243
9.486
9.327
0.5
3.644
3.351
2.499
14.554
12.549
12.327
0.6
5.152
4.681
4.017
18.674
16.679
16.434
0.7
7.295
6.543
6.214
23.678
22.215
22.098
0.8
10.344
9.153
9.315
29.574
29.589
29.904
0.9
14.691
12.820
13.580
36.347
39.350
40.804
μ
5.376
4.753
4.272
17.154
16.335
16.427
σ
4.445
3.826
4.357
10.202
11.272
11.678Furthermore, in Figure4, the RMSE between both methods is shown, while increasing the sample rate for the oscillator signals, at α values of 0.3, 0.5, and 0.7. This RMSE value decreases with higher sample rates, being more stable for sample rates higher than 1 kHz.Figure 4
Root mean-squared error of different fractional derivatives with different sample rates.Finally, we extended this analysis by doing a comparison between both GL and RL approaches at different fractional derivatives and sample rates for both signals (Figure5). In Figure 5(a), we compared the damped oscillator signal, with both approaches at 10, 50, and 100 Hz and with fractional derivatives at 0.3, 0.5, and 0.7. For 10 Hz (left column), our method has performance issues due to the signal being sampled at the Nyquist frequency. SymPy’s RL method does not have this limitation because it requires the original equation and the evaluation time (t), ignoring any sampling.Figure 5
Comparison between fractional derivatives atα of 0.3, 0.5, and 0.7 (top to bottom rows) and for sample rates of 10, 50, and 100 Hz (left to right columns) for both (a) and (b) sine-cosine signals. (a) Comparison for oscillator signal. (b) Comparison for sine-cosine signal.
(a)
(b)In Table1, the mean RMSE between α values of 0.1 and 0.9 is shown. When the sample rate increases, our GL implementation for the oscillator has a mean RSME in this experiment of 4.753 for the derivatives at 50 Hz and 4.272 at 100 Hz compared to the RL approach. For the sine-cosine signal, we have a similar performance, as shown on Figure 5(b).Table2 shows the execution times for both GL and RL algorithms. The proposed computational method is highly optimized in execution time compared to the RL method.Table 2
Execution times, in seconds, for GL and RL algorithms usingα=0.3,0.5,and0.7 for both oscillator and sine-cosine signals at different samples rates.
α
Oscillator (seconds)
Sine-cosine (seconds)
10 Hz
50 Hz
100 Hz
10 Hz
50 Hz
100 Hz
0.3 GL
0.0012
0.0010
0.0019
0.0009
0.0010
0.0015
0.3 RL
26.6902
132.2715
255.0312
107.1661
512.4475
1006.4040
0.5 GL
0.0003
0.0010
0.0014
0.0004
0.0010
0.0015
0.5 RL
26.0732
129.7112
246.8918
105.5369
507.6593
975.6912
0.7 GL
0.0003
0.0010
0.0015
0.0003
0.0010
0.0019
0.7 RL
26.2456
130.8888
250.3541
106.9648
510.5240
990.7355


## 5. Fractional Derivative Applied to Biosignals

Studies in the literature have quantified that there is a relationship between the photoplethismography (PPG) and the arterial blood pressure (ABP) signals. In [1], we have shown the mathematical relation between these two signals. To accomplish this, we propose to apply the computational method of the fractional derivative in order to find a relationship between these biosignals of interest.We assume that theα fractional derivative of the PPG signal resembles the arterial blood pressure signal as follows:(30)ABP=DαPPGn.To fit the fractional alpha derivative of PPG to the ABP signal, we need to optimize the alpha parameter by minimizing the least square error. Figure6 shows a sequence of derivatives ranging from the PPG signal from α=0 until the optimal value of α, where the ABP signal is similar.Figure 6
Fitting from PPG to BP by fractional derivatives, the parameterα controls the degree of the derivative.The fractional derivatives are able to capture the characteristics and properties of the relations between two biosignals. The fractional approach opens new possibilities in the field of medicine, where the signals may have a symbiotic relationship between them. In other words, one signal may contain information from another, and this leads us to think of methods for extracting that hidden information.The signals used in this section were extracted from the “Training set for non-Invasive and minimally-Intrusive (nImI) blood pressure estimates” from the websitehttp://nimi.uv.cl [25].


## 6. Conclusion

In this paper, we have applied the Grünwald–Letkikov (GL) differintegral operator as a computational method for signal processing. The experiments show that the Cauchy finite convolution with the GL dynamic kernel can be applied to any signal without knowing its analytical form. Furthermore, when the signal’s analytical form is known, we empirically showed that Grünwald–Letkikov and Riemann–Liouville (RL) methods obtain similar results, but the execution time of the GL method outperforms the RL method with several orders of magnitude. It is important to note that the GL computational method is affected when the signal sampling rate is close to the Nyquist rate bound; meanwhile, the RL method still works. However, the RL method requires the analytical equation of the signal that it is hardly ever available in real applications.Finally, the proposed computational method of the GL derivative of the signal can be applied for feature extraction and processing of any stochastic signals.The fractional calculation could be useful for extracting information from signals that are related to each other.



---

*Source: 7280306-2018-08-01.xml*