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.

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 [15].

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 [68]); 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 and k-th iterated integral are extended by considering a fractional parameter instead of the integer 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 [1316].

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 of Cauchy 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 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.

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 of with respect to x is the function , defined as

The derivative of a function represents an infinitesimal change in the function with respect to one of its variables. For higher-order derivatives, let be the k-th derivative, and the k-th derivative is defined by taking the derivative of the previous function as follows:

If we apply the derivative recursively from the function , we obtain the following expression:where is the binomial coefficient indexed by a pair of integers, and , . If the notation 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.

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 . 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:where the following properties hold: , and if is an integer number, then . The previous equation converges for all . If we replace the argument with a real positive parameter 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 and 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:where the binomial coefficient is expanded as follows:

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 is , and the being . Thus, the RL fractional integral is defined as

And, the fractional derivative of the RL is given bywhere is the evaluation point, α is the fractional order, and is the reference point of integration.

Several studies have shown important properties and conditions for the existence of the fractional models [1719]. 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 be a subset of . A set of stochastic or random variables indexed by is called a stochastic continuous-time process. When , the stochastic discrete-time process is obtained. The signals and 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 equals to 0 when .

On a continuous state, a random variable can take any value between a time interval . The fractional derivative for a continuous signal can be expressed as follows:where is the spacing between points on a grid and directly related to the approximation error . Moreover, is the upper limit of differentiation, and in this case, is considered.

A continuous signal can be converted into a discrete signal to a sequence of samples by sampling the signal every seconds. The sampled signal is related to the continuous signal through the equation , where is called the sampling interval or sampling period. The sampling frequency is defined as the number of samples taken per second, thus . 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 signal , , can be expressed as a vector:where is the number of sample points. The fractional derivative of a discrete signal 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 depends only on , as follows:

Let us define the Grünwald–Letnikov (GL) dynamic kernel of size n aswhere the coefficients , are given by (17).

The n-th Cauchy finite discrete convolution is defined as follows:

The Grünwald–Letnikov fractional derivative can be rewritten by computing the Cauchy finite discrete convolution, given by (20), between the signal and the Grünwald–Letnikov (GL) dynamic kernel (19) as follows:

In Figure 1, the values of the coefficients of the kernel , for and , are shown. Notice that, when 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 changes to the following expression:and the constant changes to

The 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 . We have optimized the numerical computation of the kernel (17) with a recursive equation:

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.

Algorithm 1 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 , the first term of the dynamic kernel , and the constant . 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.

Discrete signal
Sampling rate
Number of samples in the signal
user defined value
for to do
  
end for
for to do
   
for to n do
   
  end for
end for
Return the array D of the fractional derivative

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]. The fractional derivative gives insight about the functional relationship for modelling complex systems. For instance, the application of a fractional derivative with 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:with being the damping coefficient and , with being its frequency in Hz. For the experiment, we determined  Hz and the damping coefficient, . If we apply the Hermmann rules [6] to (27), we are able to obtain the following -th fractional derivative equation for this signal:

The second signal, named by us as “sine-cosine function,” is modelled with the following equation:with . 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 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.

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. Figure 3 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.

Table 1 shows the values obtained for Figure 3. For , 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. This 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 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.

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. The proposed computational method is highly optimized in execution time compared to the RL method.

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:

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 until the optimal value of α, where the ABP signal is similar.

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 website http://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.

Data Availability

The 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

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors acknowledge the support of CONICYT + PAI/CONCURSO NACIONAL INSERCIÓN EN LA ACADEMIA, CONVOCATORIA 2014 + Folio (79140057), FONDEF ID16I10322, and REDI170367 from CONICYT.