A Computational Fractional Signal Derivative Method

. We propose an eﬃcient computational method to obtain the fractional derivative of a digital signal. The proposal consists of a new interpretation of the Gr¨unwald–Letnikov diﬀerintegral operator where we have introduced a ﬁnite Cauchy convolution with the Gr¨unwald–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¨unwald–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¨unwald–Letnikov method outperforms the other approach in execution time. Finally, we show an application of how our proposal can be useful to ﬁnd the fractional relationship between two well-known biomedical signals.


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][2][3][4][5].
e 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][7][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 d k y/dx k and k-th iterated integral . . . ·dx are extended by considering a fractional α ∈ R + parameter instead of the integer k ∈ N parameter.
On the other hand, the numerical computation of the derivatives of the signal has many uses in analytical signal processing. e 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). e derivatives have been used to smoothen the signals by filtering the noise, for peak detection, trace analysis, and feature extraction, among others [13][14][15][16]. e aim of this article is to propose an efficient computational method of the fractional calculus applications to digital signal processing.
us, we introduce a numerical formulation of the fractional derivatives for continuous and discrete signals, and also, we introduce the concepts of Cauchy finite convolution and the Grünwald-Letnikov (GL) dynamic kernel. e 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. is article is organized as follows. In Section 2, 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.

Basic Principles of Fractional Calculus
e 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 of f(x) with respect to x is the function f ′ (x), defined as e derivative of a function represents an infinitesimal change in the function with respect to one of its variables. For higher-order derivatives, let f (k) be the k-th derivative, and the k-th derivative is defined by taking the derivative of the previous function f (k−1) as follows: If we apply the derivative recursively from the function f(x), we obtain the following expression: where k m 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: 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: In this case, the binomial coefficient for a negative value is computed as follows: In Sections 2.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.

Grünwald-Letnikov Differintegral Operator.
e Grünwald-Letnikov differintegral is a combined differentiation/integration operator. e Grünwald-Letnikov e GL operator is based on a classical concept for which integer order derivatives can be represented as limits of finite differences (3). e binomial coefficient of (4) can be generalized to real arguments by means of the Euler's gamma function: where the following properties hold: Γ(1) � 1, Γ(x + 1) � xΓ(x) and if k is an integer number, then Γ(k + 1) � k!. e previous equation converges for all x ∈ C, (R(x) > 0). If we replace the k argument with a real positive parameter α ∈ R + in (4), then we obtain the following property: Grünwald-Letnikov introduced (8) in (3), obtaining the following definition for the fractional α-th derivative: 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 α ∈ R + , we consider the −α for the GL integral definition to obtain the following expression: where the binomial coefficient is expanded as follows:

Riemann-Liouville Differintegral
Operator. On the other hand, Riemann and Liouville have extended traditional integral in order to obtain its fractional model. e notation for the α − integral is J α f, and the α − differintegral being And, the fractional derivative of the RL is given by (13) where x is the evaluation point, α is the fractional order, and x 0 is the reference point of integration.
Several studies have shown important properties and conditions for the existence of the fractional models [17][18][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].

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 X(t) { } t∈T indexed by T is called a stochastic continuous-time process. When T � N, the stochastic discrete-time process X[n] { } n∈T is obtained. e signals x(t) and x[n] are defined as a realization of a realvalued continuous or discrete random process, respectively. In this work, we will consider only causal signals, and the value of a signal x(t) equals to 0 when t < 0.
On a continuous state, a random variable can take any value between a time interval t ∈ T. e fractional derivative for a continuous signal can be expressed as follows: where mh is the spacing between points on a grid and directly related to the approximation error O(h α ). 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 every T s seconds. e sampled signal x[n] is related to the continuous signal through the equation t � nT s , where T s is called the sampling interval or sampling period. e sampling frequency f s is defined as the number of samples taken per second, thus f s � 1/T s . In order to capture all the information from the continuous signal, it is necessary to fulfill the Nyquist sampling theorem [21]. e 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. e discrete signal x[n], n � [0, . . . , N − 1], can be expressed as a vector: where N is the number of sample points. e fractional derivative of a discrete signal x[n] is given by In order to optimize the definition above, it is possible to define the α-th Grünwald-Letnikov kernel function as follows: where the value of constant C α depends only on α, as follows: where the coefficients G α (m), m � 0, . . . , n − 1 are given by (17). e n-th Cauchy finite discrete convolution is defined as follows: e Grünwald-Letnikov fractional derivative can be rewritten by computing the Cauchy finite discrete convolution, given by (20), between the signal x[n] and the Grünwald-Letnikov (GL) dynamic kernel G α n (19) as follows: In Figure 1, 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.
Notice that, for the integral case, the α-th Grünwald-Letnikov kernel function G α changes to the following expression: and the constant C α changes to e Grünwald-Letnikov fractional integral can be rewritten by computing the Cauchy finite discrete convolution as follows: To numerically compute the fractional derivative of a signal, we apply the following property to the gamma function of (17): where the following well-known property was used Γ(α + 1) αΓ(α). We have optimized the numerical computation of the kernel (17) with a recursive equation: is 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.
Algorithm 1 shows the implementation of the Grünwald-Letnikov fractional derivative as the Cauchynite 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 de ned α parameter of the derivative, respectively. Line 5 through 7 computes several constant terms: the Γ(α + 1), the rst 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 rst component of the fractional derivative. Between lines 13-15, the Cauchy nite discrete convolution between the signal and the dynamic kernel is computed. Line 17 gives the result of the fractional derivative of a signal.
It 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]. e 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 rst derivative; nevertheless, the fractional derivative captures the nonlinear behaviour and long-term memory due to its inherent de nition [23].

Simulation Study
For this study, we implemented two di erentiable signals in order to calculate their fractional derivatives with di erent α 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: with c being the damping coe cient and ω 2πf, with f being its frequency in Hz. For the experiment, we determined ω 5 Hz and the damping coe cient, c 1. If we apply the Hermmann rules [6] to (27), we are able to obtain the following α-th fractional derivative equation for this signal: e second signal, named by us as "sine-cosine function," is modelled with the following equation: with ω 2π. In order to make the equation continuous in R, 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 sinecosine function, as presented on Figure 2, and both signals were sampled at 100 Hz. We can see that both methods have the same behaviour for different α in both signals. e 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.
is comparison was done in order to study the effects of the sampling rate when calculating the fractional derivative on discrete signals with both methods. Figure 3 shows the increase on the RMSE when doing higher fractional derivatives. e Nyquist frequency is 10 Hz for the oscillator, At this sample rate, the RMSE is the highest for almost all derivatives. Table 1 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.
Furthermore, in Figure 4, 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. is RMSE value decreases with higher sample rates, being more stable for sample rates higher than 1 kHz.
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 ( Figure 5). In Figure 5(a), we compared the damped oscillator signal, with both approaches at 10, 50, and 100 Hz and with fractional In Table 1, 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). Table 2 shows the execution times for both GL and RL algorithms. e proposed computational method is highly optimized in execution time compared to the RL method.

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: 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. Figure 6 shows a sequence of derivatives ranging from the PPG signal from α � 0 until the optimal value of α, where the ABP signal is similar. e fractional derivatives are able to capture the characteristics and properties of the relations between two biosignals. e 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. e signals used in this section were extracted from the "Training set for non-Invasive and minimally-Intrusive (nImI) blood pressure estimates" from the website http://nimi.uv.cl [25].

Conclusion
In this paper, we have applied the Grünwald-Letkikov (GL) differintegral operator as a computational method for signal processing.
e 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 a ected 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. e fractional calculation could be useful for extracting information from signals that are related to each other.
Data Availability e training set for non-Invasive and minimally-Intrusive (nImI) blood pressure estimates can be found at http://nimi. uv.cl; for more information, contact the corresponding author.

Conflicts of Interest
e authors declare that they have no con icts of interest.