Discrete-Time Fractional Difference Calculus: Origins, Evolutions, and New Formalisms

1. Introduction

All everyday human activities give rise to signals that carry a certain type of information about the systems that generated them. These signals are bounded functions that are collected to be studied, transmitted and manipulated in order to extract the information they carry. Discrete-Time Signal Processing (DTSP), also called Digital Signal Processing, is a set of mathematical and engineering techniques that allow the processing (collection, study, analysis, synthesis, transformation, storage, etc.) of signals performed mainly on digital devices.

Combining ideas, theories, algorithms and technologies from different quadrants, the DTSP has not stopped continuously evolving and increasing its already vast field of applications. This evolution was motivated by the enormous progress of digital technologies that allow the construction of processors, in general, more reliable and robust than analog ones and, above all, more flexible. The on-chip implementation of specialized processors (e.g., FFT) has facilitated the application of mathematical techniques that would be difficult (or impossible) to perform analogically. The DTSP plays an important role in communication systems where its mission is to handle signals, both at transmission and reception, in order to achieve an efficient and reliable flow of information between source and receiver. However, it is not only in communication systems that we find DTSP applications. In fact, its field of action has widened and includes areas such as: Speech, Radar and Sonar, Seismology, Bio-medicine, Economics, Astronomy, etc. In Mathematics, it has been very useful in the study of functions and in solving differential equations. The well-known Newton series is very famous.

Mathematically, the DTSP relied on several important tools such as real and complex analysis, difference equations, discrete-time Fourier and Z transforms, algebra, etc. It benefited from the enormous development of signal theory in the 2nd half of the 20th century when signal processing techniques reached a sufficiently high degree of development. However, its origins are much earlier.

In general, we can "date" the beginning of the study of signals to the discovery of periodic phenomena that led to the introduction of the notions of year, week, day, hour, etc. With an equal degree of importance, we can consider the theory and representation of music made by the Pythagoreans as the first spectral analysis. It is important to note that they actually made a discrete time-frequency formulation. More recently, we refer the discovery and study of the spectrum of sunlight by Newton (1666) and the works of mathematicians such as Euler, Gauss (who devised the first algorithm for the fast Fourier transform – 1805), Fourier (who created the basis for spectral analysis), Sturm and Liouville. These works had direct implications on the way of studying signals in the frequency domain, which did not cease to evolve and gain importance from the 1940s, thanks to the works of the theoretical field of stochastic processes (Wiener and Kolmogorov): correlation, adapted filter, Wiener filter, etc [1,2], notions that would become the basis of modern developments in Spectral Analysis (Tukey, Parzen, Akaike, Papoulis, and Burg). It was also Tukey who, with Cooley, rediscovered the algorithm that allowed the implementation of the FFT in 1965, which was a milestone in signal analysis.

The difference equations, taking the form of the ARMA (autoregressive-moving average) model, had a rapid increase in importance due to the works of Box, Jenkins, Oppenheim, Kailath, Kalman, Robinson, Rabiner, and many others in the 1980s of the 20th century [3,4,5,6,7,8]. We can place here the real beginning and affirmation of DTSP. Nevertheless, the discovery of computers was perhaps the biggest impulse given to the DTSP, by the possibility of discrete implementation of processor devices, previously made exclusively with analog technology and to perform simulations that allow to predict, with great accuracy, the behavior of a given system. This led to an autonomization of the theory of “Discrete Signals and Systems” that became an independent branch, leading to alternative technological solutions based on digital design and realization devices [3,4,8,9,10,11,12,13,14,15]. Although the main developments were based on difference equations, the true origins were not forgotten and motivated some attempts to model and identify systems based on the delta difference [16,17,18,19,20].

The emergence of fractional tools has opened new doors to the modeling and estimation of everyday systems that were known to be best described by fractional systems. However, this does not mean that there was a coherent theory of fractional systems in discrete time. Probably the first attempt was made in [21], but the systems described are not really fractional, although they use fractional delays. In the last 20 years, many texts have been published on fractional differences and derivatives in discrete time, leading to different views of what fractional systems in discrete time are and how they are characterized [22,23,24,25,26,27,28]. The purpose of this paper is exactly to describe the mathematical basis underlying the main formulations. We introduce differences through a system approach to highlight the fact that the required definition must be valid for any function regardless of its support. This allows for a broader scope. On the other hand, it is important to make a clear distinction between time flow from left to right (causal case) or the other way around (anti-causal). Under normal conditions, they should not be mixed. To this end, we define "time sequence" as an alternative to "time scale", avoiding the confusion that the latter might introduce. We will proceed with the definitions of nabla (causal) and delta (anti-causal) differences and enumerate their main properties. We proceed with the introduction of other formulations, such as discrete-time, bilateral, tempered differences, and the completely new bilinear differences. These differences are “invariant under translations”. We propose new “scale-invariant differences” that are connected to Hadamard derivatives. For all the presented differences, ARMA type difference equations are proposed.

Given the importance of discrete signals inherent in this work, we review the classical sampling theorem valid for the case of shift-invariant systems [14,29,30,31,32] and another one suitable for scale-invariant systems, but different from similar existing in the literature [33,34].

The paper outlines as follows. In Section 2 we present several mathematical tools useful in the paper and clarify some notions. The sampling theorems are introduced here. Section 3 is used to make a historical overview of the difference evolutions, both continuous and discrete -time. The different approaches are described. The problems created by some definitions are criticized in Section 4. The main contributions in this paper are presented in Section 5 where several shift-invariant differencers and accumulators, say: nabla, delta, two-sided, tempered, and bilinear differences. For all the definitions, continuous-time and discrete-time versions are presented. The scale-invariant differences are introduced and studied in Section 6. All the described differences are suitable for defining ARMA type linear systems. This is done and exemplified in Section 7. Finally we present a brief discusion.

2. Preliminaries

2.1. Glossary and Assumptions

In the evolution of the DTSP several notions have been introduced without retaining a clear meaning. In fractional calculus there is considerable confusion in the terminology adopted having, in some cases, the same name for different operators. Here we try to clarify the meaning of some terms in order to avoid confusion. Therefore, we start with some fundamental terms. Later we will introduce others, which are needed in the rest of the document [35].

• Anti-causal
  An anti-causal system is causal under reverse time flow. A system is anti-causal if the output at any instant depends only on values of the input and/or output at the present and future time instants. The delta derivative is an example of anti-causal system.
• Anti-difference
  The operator that is simultaneously the left and right inverse of the difference will be called anti-difference.
• Backward
  Reverse time flow—from future to past.
• Causal operator or system
  A system is causal if the output at any instant depends only on values of the input and/or output at the present and past instants. The nabla derivative is an example of causal system.
• Forward
  Normal time flow—from past to future.
• Fractional
  Fractional will have the meaning of non integer real number.
• Scale-invariant system
  A system is scale-invariant if a stretching or shrinking in the input produces the same stretching/shrinking in the output. It is described by the Mellin convolution [33,36]

  $$y\left(\tau \right)=x\left(\tau \right)★g\left(\tau \right)={\int }_0^{\infty }x\left(\frac{\tau }{\eta }\right)g\left(\eta \right)\frac{d\eta }{\eta }.$$
• Signal
  Bounded function that conveys some kind of information.
• Shift-invariant system
  A system is shift-invariant if a delay or lead in the input produces the same delay/lead in the output. It is described by the usual convolution [37]

  $$y\left(t\right)=x\left(t\right)\ast g\left(t\right)={\int }_{-\infty }^{+\infty }x(t-\eta )g\left(\eta \right)\mathrm{d}\eta .$$
• System
  Any operator that transforms signals into signals. We will often use the terms system and operator interchangeably.

2.2. Some mathematical tools

Traditionally, the delta symbol is used for the so-called “forward difference” and it comes from long time ago, while nabla uses to be attached to the “backward difference” [38], in contradiction with the time flow. However, they are generalized for any order through the same way: the binomial theorem [39]. Let $\alpha \in \mathbb{R}-{\mathbb{Z}}^{-}$. Then,

$${(1-z)}^{\alpha }=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)z^k,\left|z\right|<1.$$

(1)

We can extend it for negative integer values of $\alpha$, through the Pochhammer symbol.

The binomial coefficients

$$\left(\genfrac{}{}{0pt}{}{\alpha }{\beta }\right)={\displaystyle \frac{\Gamma (\alpha +1)}{\Gamma (\alpha -\beta +1)\Gamma (\beta +1)}}$$

assume a central role in the theory we will develop in the following. They enjoy interesting properties [40,41]:

$$\left|\left(\genfrac{}{}{0pt}{}{\alpha }{\beta }\right)\right|\le {\displaystyle \frac{A}{{\beta }^{\alpha +1}}},\mathrm{for}\beta \to \infty ,$$

with $A>0$.

$$\left(\genfrac{}{}{0pt}{}{\alpha }{n}\right)={\displaystyle \frac{\Gamma (\alpha +1)}{\Gamma (\alpha -n+1)n!}}={(-1)}^n{\displaystyle \frac{{(-\alpha )}_n}{n!}},n\in \mathbb{N},$$

where ${\left(a\right)}_n$ is the Pochhammer symbol for the rising factorial

$${\left(a\right)}_n=\prod _{k=0}^{n-1}(a+k)={\displaystyle \frac{\Gamma (a+n)}{\Gamma \left(a\right)}},\mathrm{with}{\left(a\right)}_0=1,$$

generalized as

$${\left(a\right)}_{\beta }={\displaystyle \frac{\Gamma (a+\beta )}{\Gamma \left(a\right)}},\mathrm{with}{\left(a\right)}_0=1.$$

$$\left(\genfrac{}{}{0pt}{}{\alpha }{n}\right)={(-1)}^n\left(\genfrac{}{}{0pt}{}{n-\alpha -1}{n}\right).$$

$$\left(\genfrac{}{}{0pt}{}{\alpha +n}{n}\right)=\left(\genfrac{}{}{0pt}{}{\alpha +n}{\alpha }\right).$$

$$\sum _{m=0}^{\infty }\left(\genfrac{}{}{0pt}{}{\alpha }{m}\right)\left(\genfrac{}{}{0pt}{}{\beta }{n-m}\right)=\left(\genfrac{}{}{0pt}{}{\alpha +\beta }{n}\right).$$

The falling factorial is represented by [38]

$${\left(a\right)}^{\left(n\right)}=\prod _{k=0}^{n-1}(a-k)={\displaystyle \frac{\Gamma (a+1)}{\Gamma (a+1-n)}},\mathrm{with}{\left(a\right)}_0=1,$$

so that

$$\left(\genfrac{}{}{0pt}{}{\alpha }{n}\right)={\displaystyle \frac{{\left(\alpha \right)}^{\left(n\right)}}{n!}},n\in \mathbb{N},$$

and

$${\left(a\right)}^{\left(n\right)}={(-1)}^n{(-a)}_n.$$

It is generalized by

$${\left(a\right)}^{\left(\beta \right)}={\displaystyle \frac{\Gamma (a+1)}{\Gamma (a+1-\beta )}},\mathrm{with}{\left(a\right)}_0=1,$$

and

$${\left(a\right)}^{\left(\beta \right)}={(a-\beta +1)}_{\beta }.$$

The (bilateral) Laplace transform (LT) is given [42]:

$$\mathcal{L}\left[h\left(t\right)\right]=H\left(s\right)=\underset{-\infty }{\overset{\infty }{\int }}h\left(t\right)e^{-st}\mathrm{d}t,s\in \mathbb{C},$$

(2)

that is assumed to converge in some non void region (region of convergence – ROC) which may degenerate into the imaginary axis, givig rise to the Fourier transform (with $s=i\omega$.) We define the inverse LT by the Bromwich integral

$$h\left(t\right)={\mathcal{L}}^{-1}F\left(s\right)=\frac{1}{2\pi i}\underset{a-i\infty }{\overset{a+i\infty }{\int }}H\left(s\right)e^{st}\mathrm{d}s,t\in \mathbb{R},$$

(3)

where $a\in \mathbb{R}$ is called abscissa of convergence. Frequently, we denote by $\gamma$ the integration path.

In a similar way, we define the Mellin transform (MT) by [36]

$$\mathcal{M}\left[h\left(t\right)\right]=H\left(v\right)=\underset{0}{\overset{\infty }{\int }}h\left(t\right)t^{-v-1}\mathrm{d}t,v\in \mathbb{C},$$

(4)

with inverse similar to (3)

$$h\left(t\right)={\mathcal{M}}^{-1}H\left(v\right)=\frac{1}{2\pi i}\underset{\gamma }{\int }H\left(v\right)t^v\mathrm{d}v,t\in {\mathbb{R}}^{+}.$$

(5)

For the discrete-time case, we define the Z transform [14,41] by

$$\mathcal{Z}\left[f\left(n\right)\right]=F\left(z\right)=\sum _{n=-\infty }^{+\infty }f\left(n\right)z^{-n},z\in \mathbb{C},$$

(6)

with the inverse given by the Cauchy integral

$$f\left(n\right)=\frac{1}{2\pi i}{\oint }_{c}F\left(z\right)z^{n-1}\mathrm{d}z,$$

(7)

where c is the unit circle. With the change of variable $z=e^{i\omega },-\pi <\omega \le \pi$ we obtain the discrete-time Fourier transform.

2.3. On Time Sequences

A powerful approach into the continuous/discrete unification and generalization was introduced by Aulbach and Hilger through the calculus on time scales [43,44]. These are nonempty closed subsets $\mathbb{T}$ of the set $\mathbb{R}$ of real numbers. Let t be the current instant. Using the language of the time scale calculus, the previous next instant is denoted by $\rho \left(t\right)$. Similarly, the next following point on the time scale $\mathbb{T}$ is denoted by $\sigma \left(t\right)$. One has

$$\rho \left(t\right)=t-\nu \left(t\right),\sigma \left(t\right)=t+\mu \left(t\right),$$

where $\nu \left(t\right)$ and $\mu \left(t\right)$ are called the graininess functions. These functions can be used to construct any time sequences. However, we will no continue this way.

Remark 1. 

The designation “time scale” is misleading, since the word “scale” is associated to a notion of stretching or shrinking, frequently having a relation to a speed or a rate. For example. consider a function $f\left(t\right)$ defined on $\mathbb{R}$ and a parameter $a>0$ which allow us to define a new function $g\left(t\right)=f\left(at\right)$. We modified the way how the flux of information is delivered. An interesting example is given by the classic turntables where we are able to switch from a rotation speed to another one. This parameter is usually called scale. Therefore, we propose here the use of the designation time sequence.

In this work, we will consider time sequences $\mathbb{T}$ defined by a set of discrete instants $t_n$, $n\in \mathbb{Z}$, and by the corresponding graininess functions. We define a direct graininess [45,46]

$$t_n=t_{n-1}+{\nu }_n,n\in \mathbb{Z},$$

and reverse graininess

$$t_n=t_{n+1}-{\mu }_n,n\in \mathbb{Z},$$

where we avoid representing any reference instant $t_0$. These definitions of “irregular” sequences are suitable for dealing with some of the most interesting time sequences we find in practice. However, we have some difficulties in doing some kind of manipulations that are also very common. Let us consider a time sequence defined on $\mathbb{R}$ and unbounded when $t\to \pm \infty$. For example, a time sequence defined by

$$t_n=nT+{\tau }_n,n\in \mathbb{Z},T>0,\left|{\tau }_n\right|<\frac{T}{2},$$

which we can call “almost linear sequence” [47]. However, in the most interesting engineering applications, we consider regular (uniform) sequences

$$\mathbb{T}=h\mathbb{Z}=\left\{\dots ,-3h,-2h,-h,0,h,2h,3h,\dots \right\},$$

with $h\in {\mathbb{R}}^{+}$.

Remark 2. 

We can consider a slided time sequence by a given value, $a+h\mathbb{Z},a<h,$ but this corresponds to introducing another parameter that we cannot determine, due to the relativistic character of any time sequence. In other words: we need another time sequence not depending on a to fix it. However, this may be an acceptable procedure for studying continuous-time functions.

Now, consider a power transformation of a time sequence:

$$\theta =q^t.$$

We generate a new (scale)–sequence which is in ${\mathbb{R}}^{+}$. In particular, we will obtain sequences like

$${\theta }_n={\theta }^n,n\in \mathbb{Z},\theta >0,$$

or

$${\theta }_n={\theta }_{n-1}·{\tau }_n,n\in \mathbb{Z},{\tau }_n>0.$$

We will use these sequences when dealing with scale-invariant differences.

2.4. On the Sampling

Let $f\left(t\right),t\in \mathbb{R}$ be a continuous-time bounded function. The discrete-time function obtained from $f\left(t\right)$ by retaining the values at a set of pre-specified instants is the sampled function, $f\left(t_n\right),n\in \mathbb{Z}$. The procedure for obtaining such function is called ideal sampling. From an operator (system) point of view this is obtained with help of the comb distribution. Although we can consider irregular combs, we will not do it here [47]. The objectives in mind lead us to the usual comb. This is a periodic repetition of the Dirac delta function [48,49,50,51]. Here, we state it in the following format.

$$c\left(t\right)=\sum _{n=-\infty }^{+\infty }\delta (\frac{t}{T}-n),$$

(8)

where T is the sampling interval. The Fourier transform (FT) of this function is also a periodic comb.

$$FT\left[\sum _{n=-\infty }^{+\infty }\delta (t-nT)\right]=\frac{2\pi }{T}\sum _{m=-\infty }^{+\infty }\delta (\omega -m\frac{2\pi }{T}).$$

(9)

The comb is called the ideal sampler because, when multiplying a given function, $x\left(t\right)$, by a comb, $c\left(t\right)$, it retains the samples of the original function, giving rise to a modulated comb:

$$x_s\left(t\right)=x\left(t\right)·c\left(t\right)=\sum _{n=-\infty }^{+\infty }x\left(t\right)\delta (\frac{t}{T}-n)=T\sum _{n=-\infty }^{+\infty }x\left(nT\right)\delta (t-nT).$$

(10)

If $x\left(t\right)$ has a jump at $t=n_{0}T$, we use the half sum of the lateral limits $x\left(n_{0}T\right)=\frac{x\left(n_0{T}^{+}\right)+x\left(n_0{T}^{-}\right)}{2}$, in agreement with the inverse Laplace and Fourier integrals. Let $X\left(s\right)$ be the LT of $x\left(t\right)$. Then, the LT of $x_s\left(t\right)$ is given by [28]

$$X_s\left(s\right)=FT\left[x\left(t\right)·c\left(t\right)\right]=\sum _{m=-\infty }^{+\infty }X(s-im\frac{2\pi }{T}),$$

(11)

stating the well-known phenomenon: sampling in a given domain implies a repetition in the transform domain, meaning that the sampling operation produces a repetition of the transform in parallel to the real axis in strips of width $\frac{2\pi }{T}$. We observe here the reason for including T in (11): the term corresponding to $m=0$ is $X\left(s\right)$. The study we performed was based on the Laplace transform, but it can be performed also with the Fourier transform.

The choice of T depends on the objectives of the work at hand. In general we can choose any value, excepting if we have in mind the recovery of $X\left(s\right)$ from $X_s\left(s\right)$. In such case, we can impose that $X\left(s\right)$ and $X_s\left(s\right)$ have the same poles in the strip defined by $\left|Im\left(s\right)\right|<\frac{\pi }{T}$. However, the most known approach is the Whittaker–Kotel’nikov–Shannon sampling theorem [12,14,29,30,31,32]

Theorem 1. 

If a function is bandlimited to $[-\frac{\pi }{T},\frac{\pi }{T}],$ it is completely determined by giving its ordinates at a series of points spaced T seconds apart [29,30,31]:

$$f\left(t\right)=\sum _{k=-\infty }^{\infty }f\left(kT\right)sinc(\frac{t}{T}-k),$$

(12)

where

$$sinc\left(t\right)=\frac{sin\left(\pi t\right)}{\pi t},$$

with $sinc\left(0\right)=1$ is the so-called sinc function that is the impulse response of the ideal lowpass filter [14].

Consider the instant $a+nT,n\in \mathbb{Z},\gamma =a/T<1$ and denote $f\left(nT\right)$ by $f_n$. We obtain

$$f(a+nT)=\sum _{k=-\infty }^{\infty }f\left(kT\right)sinc(\frac{a}{T}+n-k)$$

and

$$f_{n+\gamma }=\sum _{k=-\infty }^{\infty }f\left(kT\right)sinc(\gamma +n-k),$$

that states what we can call fractional translation which allows expressing unknown intermediate values in terms of the uniformly spaced samples. The reverse can also be done [21].

The Whittaker–Kotel’nikov–Shannon sampling theorem is based on the usual Fourier transform that has relation with the shift-invariant systems, since it is defined in terms of the eigenfunctions of such systems: the exponentials. To obtain similar theorem for scale-invariant systems, we start from the Fourier-Mellin transform that we define by:

$$F\left(i\nu \right)={\int }_0^{\infty }f\left(\tau \right){\tau }^{-i\nu -1}d\tau ,$$

(13)

with inverse

$$f\left(\tau \right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F\left(i\nu \right){\tau }^{i\nu }d\nu ,$$

(14)

obtained from (4) and (5), by letting $v=i\nu$.

Let $Q\in {\mathbb{R}}^{+}$. As above, consider a bandlimited scale function as the one verifying

$$F\left(i\nu \right)=0,\left|\nu \right|>\frac{\pi }{Q},$$

that has an associated Fourier series

$$F\left(i\nu \right)=\sum _{k=-\infty }^{\infty }{F}_k{e}^{-ikQ\nu },$$

with

$$F_k=\frac{Q}{2\pi }{\int }_{-\frac{\pi }{Q}}^{\frac{\pi }{Q}}F\left(iu\right)e^{ikQu}\mathrm{d}u=Qf\left(e^{kQ}\right),k\in \mathbb{Z}.$$

We have successively

$$\begin{array}{cc}\hfill f\left(\tau \right)& =\frac{1}{2\pi }{\int }_{-\frac{\pi }{Q}}^{\frac{\pi }{Q}}\sum _{k=-\infty }^{\infty }{F}_k{e}^{-ikQ\nu }{\tau }^{i\nu }d\nu ,\hfill \\ & =\sum _{k=-\infty }^{\infty }{F}_k\frac{1}{2\pi }{\int }_{-\frac{\pi }{Q}}^{\frac{\pi }{Q}}{e}^{i\nu (ln(\tau )-kQ)}{\tau }^{i\nu }d\nu ,\hfill \\ & =\sum _{k=-\infty }^{\infty }{F}_k\frac{sin\left[\frac{\pi }{Q}(ln\left(\tau \right)-kQ)\right]}{\pi (ln(\tau )-kQ)},\hfill \\ & =\sum _{k=-\infty }^{\infty }{F}_k\frac{sin\left[\pi \left(\frac{ln\left(\tau \right)}{Q}-k\right)\right]}{Q\pi (\frac{ln\left(\tau \right)}{Q}-k)}.\hfill \end{array}$$

Therefore, the scale-invariant sampling theorem reads

Theorem 2. 

Let $f\left(\tau \right),\tau \in {\mathbb{R}}^{+}$ be a bandlimited scale function. It is completely determined by giving its ordinates at a series of exponentialy spaced points

$$f\left(\tau \right)=\sum _{k=-\infty }^{\infty }f\left(e^{kQ}\right)sinc\left(\frac{ln\left(\tau \right)}{Q}-k\right).$$

(15)

From this and through a comparison with theorem 1, we conclude that the scale-invariant comb is given by:

$$s\left(\tau \right)=\sum _{k=-\infty }^{\infty }\delta \left(\frac{ln\left(\tau \right)}{Q}-k\right),$$

(16)

that represents a sequence of impulses located at the points in the set $\mathcal{Q}:\tau =e^{kQ},k\in \mathbb{Z}$. These results, slightly different from similar existing in literature [33,34], will be used later in the corresponding difference definitions.

These two sampling theorems require that the functions have Fourier transforms with bounded support. However, this may not happen in practice originating what is known by “aliasing”, existing several procedures to alleviate its effect [14,29,30,31]. If we do not want to recover the continuous-time function, the sampling theorems give us a chance to reduce the losses in the continuous-to-discrete conversion.

Remark 3. 

Traditionally, we use to simplify the notation by writing $f_n=f\left(nT\right)$ and $f_{n+\gamma }=f\left((n+\gamma )T\right).$ We will adopt this procedure. For the case stated in (15), we set $q=e^Q$ and $f_{q^n}=f\left(e^{nQ}\right).$

3. Historical Overview

3.1. Euler Procedure

Newton and Leibniz introduced their (different) approaches to infinitesimal calculus in the 17th century. Leibniz’ approach was based on generalizations of sums and differences. In particular, his definition of derivative was formulated in terms of limits of increments

$$f^{\prime }\left(t\right)={\displaystyle \frac{df}{dt}}=\underset{h\to 0}{lim}{\displaystyle \frac{f\left(t\right)-f(t-h)}{h}}{f}^{\prime }\left(t\right)={\displaystyle \frac{df}{dt}}=\underset{h\to 0}{lim}{\displaystyle \frac{f(t+h)-f\left(t\right)}{h}}.$$

(17)

Euler (1768) took these formulae, removed the limit operation and used the “incremental quotients” as approximations of the derivative in solving differential equations on a set, $\mathbb{T}$, of pre-defined values of t, $\mathbb{T}=t_n,n=0,1,2,\cdots$ . With this procedure he was led to discrete functions, $f\left(t_n\right).$ It was the born of the “difference equations” and the procedure is the currently known Euler method. So the difference equations gained the protagonism that belonged to the differential equations. This procedure originated the “numerical analysis” in Mathematics and, in the XXth century, the discrete-time signal processing [13,14,15,37,41,52,53] that is a well established scientific area being responsible for important realizations in our daily life.

Anyway, the former procedure, retaining the incremental ratio, was not completely abandoned; it continued being important as an intermediate step to obtain difference equations [14] and used in some applications [17,18]. The modern approach to differential discrete equations dates back to the Hilger’s works on looking for a continuous/discrete unification [43,44,54].

To give a more precise idea and clarify the nomenclature, consider the differential equation

$${\displaystyle \frac{df\left(t\right)}{dt}}+af\left(t\right)={\displaystyle \frac{dg\left(t\right)}{dt}}+bg\left(t\right),$$

where $f\left(t\right)$ and $g\left(t\right)$ are real functions of real variable and $a,b\in \mathbb{R}.$ Using the first incremental ratio

$$\nabla f\left(t\right)={\displaystyle \frac{f\left(t\right)-f(t-h)}{h}},$$

(18)

in the equation, we obtain

$$\nabla f\left(t\right)+af\left(t\right)=\nabla g\left(t\right)+bg\left(t\right),$$

(19)

that leads to

$$(1+ah)f\left(t\right)-f(t-h)=(1+bh)g\left(t\right)-g(t-h),$$

(20)

that is a difference equation. Note that only present and past values are involved: it represents a causal system. With a similar procedure for the other derivative

$$\Delta f\left(t\right)={\displaystyle \frac{f(t+h)-f\left(t\right)}{h}},$$

(21)

we obtain

$$\Delta f\left(t\right)+af\left(t\right)=\Delta g\left(t\right)+bg\left(t\right),$$

(22)

that gives

$$f(t+h)-(1+ah)f\left(t\right)=g(t+h)-(1+ah)g\left(t\right),$$

(23)

that is again a difference equation, but involving only present and future values: anti-causal.

If we sample the function at a given discrete time grid $\mathbb{T}:t_n,n\in \mathbb{Z},$ (19) and (22) become discrete-time differential equations, while (20) and (23) transform into discrete-time difference equations.

3.2. Differences and Fractional Calculus

A fractional difference was introduced for the first time in 1832 by Liouville [55,56] who used it for defining a fractional derivative. His first formula, derived directly to be a generalization of the delta derivative, assumed the form:

$${\overline{\Delta }}_L^{\alpha }f\left(t\right)=h^{-\alpha }\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t+(\alpha -k)h).$$

(24)

This formula constituted a second step into defining fractional order derivatives. However, Liouville recognized that such a formula was not exactly what he was looking for and introduced also another one, much more interesting, that we can write as

$${\nabla }^{\alpha }f\left(t\right)=h^{-\alpha }\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t-kh).$$

(25)

He observed that, if $f\left(t\right)=e^{st},t\in \mathbb{R}$, the summation was convergent if $Re\left(s\right)>0$. This formula was the base of the fractional derivative definition used by Grünwald [57] and Letnikov [58,59]. It is important to note that Liouville definition was assumed to be valid for functions of exponential type and defined on $\mathbb{R}$ or $\mathbb{C}$, while Grünwald and Letnikov worked on $[t_0,t]\in \mathbb{R}.$

In the study of operators based on his operational methods, Heaviside obtained the same result and made a study of the binomial series which he generalised [39,60,61]. This approach was retaken later by E. Post (1930) [62] and by P. Butzer & U. Westphal (1974) [63]. Also in 1974, J. Diaz and T. Osier obtained a particular case of (24) with $h=1$ [64]. In a later section, we will perform the analysis of such a definition.

It is important to highlight that, while searching for a definition convergent for $Re\left(s\right)<0,$ Liouville arrived at the delta-type difference

$${\Delta }_L^{\alpha }f\left(t\right)={(-h)}^{-\alpha }\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t+kh).$$

(26)

However, this route to fractional calculus never gained the favor of the vast majority of researchers, so it was basically considered as an approximation and used in numerical methods. Very interesting are the integral representations for (25) and (26) [25,64] and the two sided differences [65,66]. The approach and applications introduced by V. Tarasov are also very important [23,24].

Due to its (bad) influence we are going to study the “difference” (24) that states the first attempt Liouville made to define a fractional derivative using an infinite summation. It seems that in an independent way, Diaz and Osler proposed as delta difference the expression [64]

$${\Delta }_{D+O}^{\alpha }f\left(t\right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t+\alpha -k).$$

(27)

As it is easy to observe, this formula, while agreeing with the requests in the positive integer order, fails in the fractional order, since it uses symultaneously past and future values. This fact has implications in the discrete-time differences deduced from it.

To make a fair analysis, let us apply the LT to (24). We obtain

$$F_{{\Delta }_{D+O}^{\alpha }}\left(s\right)=e^{\alpha s}\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)e^{-ks}F\left(s\right).$$

(28)

The series converges only for $Re\left(s\right)>0$, that implies causality, contrarily to what we were expecting, since the delta difference should be anti-causal. Therefore, the formula (24) does not make sense and not good consequences should be expected from it.

3.3. Discrete-Time Differences

A less likely appearance of a discrete difference happened in a different context, the summation of series using Cesàro’s method [39,67,68], where S. Chapman introduced a fractional delta discrete difference as being given by

$${\Delta }_C^{\alpha }{v}_n=\underset{p\to \infty }{lim}\sum _{k=0}^{p-n+1}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)v_{n+k},n\in \mathbb{N},$$

(29)

that was considered by G. Isaacs as alternative to Diaz and Osler’s approach, mainly due to the validity of the associativity of orders [69]. However, in 1962, Isaacs proposed another formulation reading

$${\Delta }_I^N{v}_n=\sum _{k=n}^{\infty }\left(\genfrac{}{}{0pt}{}{k-n-N-1}{k-n}\right)v_k,n\in \mathbb{N}.$$

(30)

It is important to highlight that Isaacs’ formulae express anti-causal differences as we should expect. Among several interesting results, he presented what can be considered as discrete Leibniz rule

$${\Delta }_I^N\left(v_n{w}_n\right)=\sum _{k=N}^{\infty }\left(\genfrac{}{}{0pt}{}{N}{k}\right){\Delta }_I^k{v}_n{\Delta }_I^{N-k}{w}_{n+k}.$$

Meanwhile, discrete signal (time series) analysis, which has benefited from major developments since World War II, received a spectacular boost with the publication of the landmark book ”Time Series Analysis: Forecasting and Control” by George E. P. Box and Gwilym M. Jenkins [4]. Here, the authors presented a coherent and mathematically well founded study of the ARMA models and their evolution, the Autoregressive-Moving Integrated Average models (ARIMA) [10,70]. The “Integrated” factor pointed already to what was going to happen next: the insertion of a “Fractionally Integrated” term, leading to the ARFIMA models [71,72,73,74,75,76,77]. To the formalization of this new model, the concept of fractional differencing was introduced through (1) [71,72]

$${(1-z^{-1})}^{\alpha }=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)z^{-k}=\sum _{k=0}^{\infty }\frac{{(-\alpha )}_k}{k!}{z}^{-k},$$

for $\left|z\right|>1$, meaning that

$${\nabla }^{\alpha }{v}_n=\sum _{k=0}^{\infty }\frac{{(-\alpha )}_k}{k!}{v}_{n-k},n\in \mathbb{Z},$$

(31)

that corresponds to make $h=1$ in the nabla (Liouville–)Grünwald-Letnikov difference. These two formulae were also introduced earlier by Cargo and Shisha in the study of the zeros of polynomials [78]. The fractionally differenced and integrated models have gained roots and continue being used today [10,15,79]. In engineering, namely in signal processing, a fractional generalization of the ARMA model has been proposed [25,27].

A brief glance at (31) may give us an impression of repulsion due to the summation to infinity. However, such an impression is incorrect, because, if the function is null for $n<n_0\in \mathbb{Z}$, the summation becomes finite. For example, if $v_n=0,n<0$,

$${\nabla }^{\alpha }{v}_n=\sum _{k=0}^n{(-1)}^k\frac{{(-\alpha )}_n}{n!}{v}_{n-k},n\in \mathbb{N}.$$

(32)

However, H. Gray and N. Zhang did not notice this fact and, in 1988, proposed an alternative through a new definition of the fractional difference. They had in mind to preserve many properties of fractional derivative definitions, namely the exponent law and the important Leibniz rule. To get it, they started from the summation formula that they repeated using the Cauchy procedure to obtained the n-fold summation

$${\nabla }^{-N}f\left(t\right)={\displaystyle \frac{1}{\Gamma \left(N\right)}}\sum _{k=k_0}^t{(t-k+1)}_{N-1}f\left(k\right).$$

This definition can be extended to $N\in {\mathbb{Z}}_0^{-}.$ It is important to refer that such a formula was already known from the theory of Z transform [9], since it corresponds to the inverse of a multiple pole at the point $z=1$ in the complex plane (discrete integrator or accumulator).

Gray and Zhang were led to the following definition:

Let $\alpha \in \mathbb{R}$ and $f\left(t\right)$ defined over the set $\mathbb{T}=k_0-N,k_0-N+1,\cdots ,0,1,\cdots ,t\in \mathbb{Z}.$ The $\alpha$–order summation over ${\mathbb{T}}_0=k_0,k_0+1,\cdots ,t\in \mathbb{Z}$ is defined by

$${\nabla }^{-\alpha }f\left(t\right)={\displaystyle \frac{{\nabla }^N}{\Gamma (\alpha +N)}}\sum _{k=k_0}^t{(t-k+1)}_{N+\alpha -1}f\left(k\right),$$

(33)

where $N=max0,N_0\in \mathbb{Z}$ such that $0<\alpha +N<1.$ The $\alpha$–order derivative is defined by the substitution $-\alpha \to \alpha .$ This definition is coherent with the usual integer order difference. It is interesting to note that the above definition mimics the Riemann-Liouville definition of fractional derivative [40,41,80]. It can be shown that, setting $N-1<\alpha \le N\in {\mathbb{Z}}_0^{+}$, we obtain

$$\begin{array}{cc}\hfill {\nabla }^{\alpha }f\left(t\right)& ={\displaystyle \frac{{\nabla }^N}{\Gamma (N-\alpha )}}\sum _{k=k_0}^t{(t-k+1)}_{N-\alpha -1}f\left(k\right)\hfill \\ \hfill & ={\displaystyle \frac{{\nabla }^N}{\Gamma (N-\alpha )}}\sum _{k=0}^{t-k_0}{(k+1)}_{N-\alpha -1}f(t-k).\hfill \end{array}$$

(34)

Consider the $N=0$ case and note that

$${(k+1)}_{-\alpha -1}=\frac{\Gamma (-\alpha +k)}{\Gamma (k+1)}={(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)$$

meaning that (33) coincides with (26), provided that we remove the constraint on the domain of $f\left(t\right)$ and let it be defined over $\mathbb{Z}$. On the other hand, we can prove that, if $N>0$,

$${\nabla }^N\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)=\left(\genfrac{}{}{0pt}{}{\alpha +N}{k}\right).$$

In such case, we can write, for any $\alpha$

$$\begin{array}{cc}\hfill {\nabla }^{\alpha }f\left(t\right)& ={\displaystyle \frac{{\nabla }^N}{\Gamma (N-\alpha )}}\sum _{k=0}^{\infty }{(k+1)}_{-\alpha -1}f(t-k)\hfill \\ \hfill & =\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t-k).\hfill \end{array}$$

(35)

This last expression coincides with (26), when $h=1$. We can conclude that the work of H. Gray and N. Zhang shows a different but equivalent to Liouville’s and Chapman’s way to introduce fractional differences. They presented the causal (nabla) version, but it is not very difficult to obtain the corresponding anti-causal (delta). It is important to remark that they introduced, for the first time, the discrete-time Riemann-Liouville type difference. The Caputo type difference was introduced later for both nabla and delta cases [81,82,83]. Meanwhile, K. Miller and B. Ross [84] presented an approach to delta difference that we will consider next. After these publications, it seems there existed a time gap of almost 11 years without any novelty on the subject excepting the revision of the difference concepts having in consideration the notions and unification introduced by Hilger [43,44] and the discrete-time approximations to continuous-time derivatives by I. Podlubny [85]. In [54], Bohner and Peterson revised the concepts of differences and derivatives in terms of the notion of time scale.

4. A Critical View of soMe Aspects Related to Differences

4.1. A “Fractional Delta Difference” That Is Not a Delta Difference

Concerning the fractional difference, only in 2007, F. Atici and P. Eloe returned to the subject, by considering the delta difference first [86,87] and the nabla later [88]. Although they introduce a similar formulation to Gray and Zhang’s, it seems they were unaware of it, at least initially. They intended to follow closely the continuous-time fractional derivative definitions. The starting point was Miller and Ross’ definition [84] they recovered [86,87]. However, this definition was based of the Diaz & Osler continuous-time difference that we studied above. These works were followed by many others in the last 15 years. This had, as consequence, a parallel introduction of several similar but different formalizations, even from the same author, of both nabla and delta differences, as the Riemann-Liouville like and Caputo like formulations, their compositions, and difference equations [22,81,82,86,87,89,90,91,92,93,94,95,96,97,98,99,100].

As referred, the Atici & Eloe’s approach started by using the Miller and Ross’ definition for the fractional discrete delta sum given by

$$\left({\Delta }_a^{-\alpha }f\right)\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}\sum _{u=a}^{t-\alpha }{(t-\sigma \left(u\right))}^{(\alpha -1)}f\left(u\right)$$

(36)

where $\alpha >0,$$\sigma \left(t\right)=t+1,$$f\left(t\right)$ was defined in ${\mathbb{T}}_a=a,a+1,a+2,\cdots$, and with $t=a+\alpha +k,k=0,1,2,\cdots ,N-1$. As shown by D. Mozyrska [96] it can assume the form

$${\Delta }^{-\alpha }f\left(t\right)=\sum _{m=0}^k\left(\genfrac{}{}{0pt}{}{k-m+\alpha -1}{k-m}\right)f(a+m)$$

(37)

For the nabla difference, again in disagreement with Gray & Zhang, they proposed [87,88,101]

$${\nabla }^{-\alpha }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}\sum _{u=a}^t{\displaystyle \frac{\Gamma (t-u+\alpha -1)}{\Gamma (t-u)}}f\left(u\right)$$

As above, set $u=a+m$ and $t=a+1+k$ to obtain

$${\nabla }^{-\alpha }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}\sum _{n=0}^k{\displaystyle \frac{\Gamma (k-m+\alpha -1)}{\Gamma (k-m+1)}}f(a+m)=\frac{1}{\Gamma \left(\alpha \right)}\sum _{n=0}^k{\displaystyle \frac{\Gamma (n+\alpha -1)}{\Gamma (n+1)}}f(a+k-m).$$

For other versions, see [81,89,92,95,101,102].

Remark 4. 

Strangely, the domain of the above delta sum is not ${\mathbb{T}}_a$ but ${\mathbb{T}}_{a+\alpha }$. This means that the value at a given instant, depends on past values, a causal behaviour contrarily to what we where expecting, since it is assumed to be a delta difference, which is anti-causal. Besides, most applications to engineering, economics, or statistics involve discrete-time systems where all the components, expressed in terms of sums and differences are defined on the same time sequence. However, a lot of people accepted the Miller & Ross’ formulation as correct for the delta difference with corresponding results. In fact, similar situations are repeated in other formulations [81,89,92,95,101,103]. This has a very important consequence: it is impossible to define ARMA type equations using these differences.

Example 1. 

To have an idea of the incorrectness of delta difference definition, assume that $a=0$ and $f\left(t\right)=e^{-t},t\in \mathbb{Z}$. Let us compute 3 cases:

• order 1 difference

  $$\Delta f\left(t\right)=e^{-t-1}-e^{-t}=(e^{-1}-1)e^{-t}.$$
  For a given t, the difference depends on the present and future values.
• order $-1$ difference

  $${\Delta }^{-1}f\left(t\right)=-e^{-t}-e^{-t-1}-e^{-t-2}\cdots =-e^{-t}\sum _{n=0}^{\infty }{e}^{-n}={(e^{-1}-1)}^{-1}{e}^{-t}.$$
  Again, for a given t, the $-1$–order difference (sum) depends on the present and future values.
• order $-1/2$ difference
  with $t=1/2+k$, we have

  $$\begin{array}{cc}\hfill {\Delta }^{-1/2}f\left(t\right)& =\sum _{m=0}^k\left(\genfrac{}{}{0pt}{}{k-m-1/2}{k-m}\right)e^{-m}\hfill \\ & =e^{-k}\sum _{m=0}^k\left(\genfrac{}{}{0pt}{}{n-1/2}{n}\right)e^n\hfill \\ & =e^{-t+1/2}\sum _{m=0}^k{(-1)}^m\left(\genfrac{}{}{0pt}{}{-1/2}{m}\right)e^m\hfill \\ & =e^{-(t+1/2)}+1/2e^{-(t-1/2)}-\frac{1}{8}{e}^{-(t-3/2)}\cdots \hfill \end{array}$$
  Contrarily to the above examples, the $-1/2$ order derivative depends on one future and infinite past values.

4.2. One for All or One for Each

In some discrete-time formulations, as well as in the usual continuous-time fractional calculus, it is current to attach to a given difference/derivative definition the support of the function at hand. This means that, for any function, there is a particular difference/derivative definition. It is a one to one case. This was the procedure used by Gray and Zhang. However, it creates difficulties, since we cannot add functions and difference/derivative with different durations. This is strange because it is assumed that the difference/derivative are linear operators.

Alternatively, we can define general differences/derivatives over the whole possible domain ($\mathbb{R}$ or $\mathbb{Z}$) and only particularize at the moment of the computation. To understand the situation, consider two functions

$$f\left(n\right)=\left\{\begin{array}{cc}0\hfill & n>a\hfill \\ n\hfill & a\le n\le b\hfill \\ 0\hfill & n>b.\hfill \end{array}\right.$$

$$g\left(n\right)=n,a\le n\le b.$$

At first glance, it looks like they are the same. However, they express different situations:

• $f\left(n\right)$ expresses a situation where there is a past and a future. It is like some system that exists, is in stand-by first, acts for some time, and returns to the previous state. It is the situation corresponding to many physical, biological, and social systems.
• $g\left(n\right)$, on contrary, has no past and will have no future: something is born, lives for some time and disappears.

Assume that we want to define and compute a Z transform of such functions. In the first, we use the usual definition

$$F\left(z\right)=\sum _{n=-\infty }^{\infty }f\left(n\right)z^{-n},$$

(38)

that gives

$$F\left(z\right)=\sum _{n=a}^{b}n{z}^{-n}.$$

On the other hand, for $g\left(n\right)$ we have to make a new definition

$${ }_a{G}_b\left(z\right)=\sum _{n=a}^{b}n{z}^{-n},$$

(39)

valid only for all the functions defined in the interval $a\le n\le b.$ However, $F\left(z\right)$ and ${ }_a{G}_b\left(z\right)$ are complex variable functions analytic in the same region and to which corresponds the same inversion integral that leads to the same inverse. In fact, the Cauchy integral

$$h\left(n\right)=\frac{1}{2\pi i}{\oint }_{c}H\left(z\right)z^{n-1}dz$$

defines a function $h\left(n\right)$ on $\mathbb{Z}$. This is interesting: the inverse transform imposes an extrapolation of $g\left(t\right)$ so that its domain is $\mathbb{Z}$ although its support is $a\le n\le b.$ This means that we should use (38) as definition of Z transform and consider functions of the type $f\left(n\right)$, even if its support is finite. Therefore, we do not need to specify the support in defining differences, derivatives, and transforms.

4.3. The Riemann-Liouvile and Caputo Like Procedures

In many papers and due to the influence of the traditional fractional calculus, it became current the use of Riemann-Liouvile and Caputo like procedures. However and while being mathematically correct, it is not a suitable way of doing computations, since

• we are increasing unnecessarily the number of operations; this is very important in computational implementations, because we are increasing the numerical error [14,104];
• we throw most of the computational burden on negative order binomial coefficients that behave asymptotically like $\frac{1}{n^{-\alpha +1}}$, so decreasing very slowly or even increasing.

In continuous fractional calculus, the Riemann-Liouville integral becomes singular when the order is positive (derivative case), as Liouville recognized in his first paper [105], having proposed both procedures that transfer the singular behavior to the derivative of integer order. However, such difficulty does not appear in discrete-time formulations, since the computational burden falls heavily on the binomial coefficients, directly or otherwise, which can be computed through the Pochhammer symbol that is always non singular. Therefore, there is neither particular reason nor advantage in using those procedures.

5. Shift-Invariant Differencers and Accumulators

5.1. Causal

Consider two bounded piecewise continuous functions $f\left(t\right),g\left(t\right),t\in \mathbb{R}$ with $f(-\infty )=g(-\infty )=0.$ For simplicity, assume they are of exponential order, so that we assure they have Laplace transforms, $F\left(s\right),G\left(s\right),s\in \mathbb{C}$, analytic over suitable ROC.

Definition 1. 

Let us define nabla or causal differencer as a linear system whose output, at a given instant, is the difference between the input at that instant and at a previous one:

$$g\left(t\right)=f\left(t\right)-f(t-h),$$

(40)

where $h>0$ is the delay.

The output will be called nabla difference and represented by $\nabla f\left(t\right)$. From this definition we can draw some conclusions:

1.

It is a moving-average system, sometimes called “feedforward” system;

2.

Its impulse response is given by:

$$\varphi \left(t\right)=\delta \left(t\right)-\delta (t-h),$$

so that

$$g\left(t\right)=\varphi \left(t\right)\ast f\left(t\right),$$

implying that

$$G\left(s\right)=(1-e^{-sh})F\left(s\right).$$

(41)

3.

The transfer function is

$$H\left(s\right)=1-e^{-sh},$$

(42)

having $\mathbb{C}$ as ROC.

4.

We can associate in series as many systems as we can in such a way that the output of the $(n-1)$–th system is the input of the next one, n–th

$$f_n\left(t\right)=g_{n-1}\left(t\right),g_0\left(t\right)=f\left(t\right)$$

and

$$g_n\left(t\right)=f_n\left(t\right)-f_n(t-h)=g_{n-1}\left(t\right)-g_{n-1}(t-h).$$

The transfer function of the association is given by

$$G\left(s\right)={\left[1-e^{-sh}\right]}^{n}F\left(s\right),$$

that inverted gives the n–th order nabla difference

$${\nabla }^{n}f\left(t\right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)f(t-kh).$$

(43)

Now, let us return back to (40) and invert the role of the functions: assume that the input is $g\left(t\right)$ and the output is $f\left(t\right)$ (feedback system)

$$f\left(t\right)=f(t-h)+g\left(t\right)$$

that can be reused to give:

$$f\left(t\right)=g\left(t\right)+g(t-h)+f(t-2h)=g\left(t\right)+g(t-h)+g(t-2h)+f(t-3h)=\cdots$$

It is not hard to see that

$$f\left(t\right)={\nabla }^{-1}g\left(t\right)=\sum _{n=0}^{\infty }g(t-nh),$$

(44)

with LT

$$F\left(s\right)=G\left(s\right)\sum _{n=0}^{\infty }{e}^{-nsh}={\left[1-e^{-sh}\right]}^{-1}G\left(s\right),Re\left(s\right)>0.$$

(45)

Relation (44) shows why this operator is called accumulator or sum. The series association of n equal accumulators gives the n–th order nabla sum:

$${\nabla }^{-n}g\left(t\right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{-n}{k}\right)g(t-kh).$$

(46)

Definition 2. 

This result together with (43) suggest that the α–order nabla differencer/accumulator be given by

$${\nabla }^{\alpha }f\left(t\right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t-kh),$$

(47)

where α is any real number.

The corresponding LT is

$$\mathcal{L}\left[{\nabla }^{\alpha }f\left(t\right)\right]={\left[1-e^{-sh}\right]}^{\alpha }F\left(s\right),Re\left(s\right)>0.$$

(48)

Therefore, there are some facts that we must emphasize:

1.

A difference/sum is the output of a system: differencer/accumulator;

2.

The system structure is independent of the inputs;

3.

If the order is not a positive integer, even if the input function has finite support, the output has infinite support; in particular, if $f\left(t\right)$ has support $[a,b]\subset \mathbb{R}$, $g\left(t\right)$ is not identically null above any real value: the support is $[a,\infty [$. This is a very important fact, frequently forgotten or dismissed.

4.

If the input is a right hand function, so is the output; in particular, if $f\left(t\right)=0,t<0,$ then ${\nabla }^{\alpha }f\left(t\right)=0,$ for negative t and for $t\ge 0$ we have

$${\nabla }^{\alpha }f\left(t\right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t-kh),$$

(49)

where $n=\lfloor t/h\rfloor$ is the the great integer such that $n\le t/h$.

5.

The ROC of the transfer function is defined by $Re\left(s\right)>0$, as expected, since we are dealing with a causal system.

5.2. Anti-Causal

Consider the conditions of the previous subsection, but now with with $f(\infty )=g(\infty )=0.$

Definition 3. 

We define delta or anti-causal differencer as a linear system whose output, at a given instant, is the difference between the input at that instant and at a future one:

$$g\left(t\right)=f\left(t\right)-f(t+h),$$

(50)

where $h>0$ is the lead.

The output will be called delta difference and represented by $\Delta f\left(t\right)$. Note that it is the symmetric of current definition (21). Traditionally, the symmetric of (50) is considered: $g\left(t\right)=f(t+h)-f\left(t\right).$ To simplify and highlight the similarity to the nabla difference, we found preferable to omitt the − sign.

From this definition and the previous results we conclude that:

1.

It is also a moving-average system;

2.

The LT gives

$$G\left(s\right)=(1-e^{+sh})F\left(s\right).$$

(51)

3.

The transfer function is

$$H\left(s\right)=1-e^{+sh},$$

(52)

having $\mathbb{C}$ as ROC.

4.

The association in series of n systems as above has transfer function given by

$$G\left(s\right)={\left[1-e^{sh}\right]}^{n}F\left(s\right),$$

that inverted gives the n–th order delta difference

$${\Delta }^{n}f\left(t\right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)f(t+kh).$$

(53)

5.

In a similar way, the delta accumulator is

$${\Delta }^{-1}g\left(t\right)=\sum _{n=0}^{\infty }g(t+nh).$$

(54)

The series association of n accumulators is expressed by

$${\nabla }^{-n}g\left(t\right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{-n}{k}\right)g(t-kh)$$

(55)

and has LT

$$\mathcal{L}\left[{\Delta }^{-n}g\left(t\right)\right]={\left[1-e^{sh}\right]}^{-n}G\left(s\right),Re\left(s\right)<0.$$

(56)

Definition 4. 

This result, together with (43), suggests that the α–order delta differencer/accumulator be defined by

$${\Delta }^{\alpha }f\left(t\right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(t+kh),$$

(57)

where α is any real number.

The corresponding LT is

$$\mathcal{L}\left[{\Delta }^{\alpha }f\left(t\right)\right]={\left[1-e^{sh}\right]}^{\alpha }F\left(s\right),Re\left(s\right)<0.$$

(58)

This and the nabla differences have similar characteristics. Therefore and in particular, we can say relatively to the delta difference:

1.

If the order is not a positive integer, even if the input function has finite support, the output has infinite support; in particular, if $f\left(t\right)$ has support $[a,b]\subset \mathbb{R}$, $g\left(t\right)$ is not identically null below any real value; the support is $]-\infty ,b].$

2.

If the input is a left hand function, so is the output; in particular, if $f\left(t\right)=0,t>0,$ then ${\Delta }^{\alpha }f\left(t\right)=0,$ for positive t and for $t\le 0$ we have

$${\Delta }^{\alpha }f\left(t\right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f(kh-|t\left|\right),$$

(59)

where $n=\lfloor -t/h\rfloor$ is the the great integer such that $n\le -t/h$.

3.

The ROC of the transfer function is defined by $Re\left(s\right)<0$, as expected, since we are dealing with an anti-causal system.

4.

We can account for the $(-)$ sign we removed above, by inserting the factor ${(-1)}^{\alpha }$ into (57).

5.3. Properties

The nabla and delta differencers have similar properties that we will describe for the first. The most important are [25]

• Additivity and Commutativity of the orders
  ${\nabla }^{\alpha }\left[{\nabla }^{\beta }f\left(t\right)\right]={\nabla }^{\beta }\left[{\nabla }^{\alpha }f\left(t\right)\right]={\nabla }^{\alpha +\beta }f\left(t\right).$
• Neutral element
  This comes from the last property by putting $\beta =-\alpha$, ${\nabla }^{\alpha }\left[{\nabla }^{-\alpha }f\left(t\right)\right]={\nabla }^{0}f\left(t\right)=f\left(t\right)$. This is very important because it states the existence of inverse, in coherence with the previous sub-sections.
• Inverse element
  From the last result we conclude that there is always an inverse element: for every $\alpha$ order difference, there is always a $-\alpha$ order difference given by the same formula.
• Associativity of the orders

  $${\nabla }^{\gamma }\left[{\nabla }^{\alpha }{\nabla }^{\beta }\right]f\left(t\right)={\nabla }^{\gamma +\alpha +\beta }f\left(t\right)={\nabla }^{\alpha +\beta +\gamma }f\left(t\right)={\nabla }^{\alpha }\left[{\nabla }^{\beta +\gamma }\right]f\left(t\right).$$
  It is a consequence of the additivity.
• Derivative of the product

  $${\nabla }^{\alpha }\left[f\left(t\right)g\left(t\right)\right]=\sum _{i=0}^{\infty }\left(\genfrac{}{}{0pt}{}{\alpha }{i}\right){\nabla }^{i}f\left(t\right){\nabla }^{\alpha -i}g(t-ih).$$

  (60)
  The delta case is slightly different as expected

  $${\Delta }^{\alpha }\left[f\left(t\right)g\left(t\right)\right]=\sum _{i=0}^{\infty }\left(\genfrac{}{}{0pt}{}{\alpha }{i}\right){\Delta }^{i}f\left(t\right){\Delta }^{\alpha -i}g(t+ih).$$

5.4. Discrete-Time Differences

In the Introduction, we referred the importance of discrete-time signals. In sub-Section 2.4 we showed how to obtain them by sampling continuous-time signals. However, we must highlight an important fact: discrete-time signals exist by themselves, without there needing to be a continuous-time signal from which they resulted. There are many signals that are inherently discrete. This means that, in each case, there is necessarily a clock that defines the time sequence in which we work. Therefore, the delay/lead, $h,$ must be equal or a multiple of the sampling interval used to obtain the discrete-time formulation for differences. For simplicity, we choose the sampling interval equal to h, so that, $\mathbb{T}=\left\{t=nh,n\in \mathbb{Z}\right\},$ giving

$${\nabla }^{\alpha }{f}_k=\sum _{n=0}^{\infty }{\displaystyle \frac{{(-\alpha )}_n}{n!}}{f}_{k-n},k\in \mathbb{Z},$$

(61)

and

$${\Delta }^{\alpha }{f}_k=\sum _{n=0}^{\infty }{\displaystyle \frac{{(-\alpha )}_n}{n!}}{f}_{k+n}.k\in \mathbb{Z},$$

(62)

that can assume different forms:

$${\nabla }^{\alpha }{f}_k=\sum _{n=-\infty }^k{\displaystyle \frac{{(-\alpha )}_{k-n}}{(k-n)!}}{f}_n,k\in \mathbb{Z}$$

(63)

and

$${\Delta }^{\alpha }{f}_k=\sum _{n=k}^{\infty }{\displaystyle \frac{{(-\alpha )}_{n-k}}{(n-k)!}}{f}_n,k\in \mathbb{Z}.$$

(64)

As it is obvious, the formulae (61) to (64) express input/output relations that are discrete-time convolutions between the input $f_n$ and the binomial coefficients (impulse responses) to produce the outputs that are the differences. Such relations are mediated by the called impulse responses that are the outputs, $h_n$, when the input is the Kroneckker delta defined by

$${\delta }_n=\left\{\begin{array}{cc}1\hfill & n=0\hfill \\ 0\hfill & n\ne 0.\hfill \end{array}\right.$$

In passing, we define the discrete-time step by

$${\epsilon }_n=\left\{\begin{array}{cc}1\hfill & n\ge 0\hfill \\ 0\hfill & n<0.\hfill \end{array}\right.$$

Therefore, the impulse response of the nabla and delta differences are respectively

$$h_n={\displaystyle \frac{{(-\alpha )}_n}{n!}}{\epsilon }_n$$

and

$$g_n={\displaystyle \frac{{(-\alpha )}_{-n}}{(-n)!}}{\epsilon }_{-n},$$

where $n\in \mathbb{Z}$, in agreement with the causality. It is important to remark that

1.

Both responses have finite duration if $\alpha \in \mathbb{N}$ and the systems are called FIR (finite impulse systems) [14];

2.

If $\alpha \notin \mathbb{N}$, both responses extend to infinite and the corresponding systems are IIR (infinite impulse response).

3.

If $f_k=0,k<0,$ then ${\nabla }^{\alpha }{f}_k=0,$ for negative k and for $k\ge 0$ we have

$${\nabla }^{\alpha }{f}_k=\sum _{n=0}^k{\displaystyle \frac{{(-\alpha )}_n}{n!}}{f}_{k-n},k\in \mathbb{Z}$$

(65)

or,

$${\nabla }^{\alpha }{f}_k=\sum _{n=0}^k{\displaystyle \frac{{(-\alpha )}_{k-n}}{(k-n)!}}{f}_n,k\in \mathbb{Z}.$$

(66)

4.

Similarly, if $f_k=0,k>0,$ then ${\Delta }^{\alpha }{f}_k=0,$ for positive k and we obtain for $k\le 0$

$${\Delta }^{\alpha }{f}_k=\sum _{n=0}^{\left|k\right|}{\displaystyle \frac{{(-\alpha )}_n}{n!}}{f}_{n-\left|k\right|}.k\in \mathbb{Z}$$

(67)

or

$${\Delta }^{\alpha }{f}_k=\sum _{n=0}^{\left|k\right|}{\displaystyle \frac{{(-\alpha )}_{n-k}}{(n-k)!}}{f}_n,k\in \mathbb{Z}.$$

(68)

5.

It is a simple task to obtain formulae for functions with other supports.

6.

The Z transforms of the above discrete-time differences can be obtained from the corresponding LT by setting $z=e^{sh}.$ For example, the Z transform of the nabla difference (61) is

$$\mathcal{Z}\left[{\nabla }^{\alpha }{f}_k\right]={(1-z^{-1})}^{\alpha }F\left(z\right),$$

(69)

in the suitable ROC in the region defined by $\left|z\right|>1$.

7.

If, in any particular application, a time sequence with the form $\mathbb{T}=a+nh,;n\in \mathbb{Z},a\in \mathbb{R}$ is used, we can make a substitution $f(a+nh)$ for $f\left(nh\right)$.

5.5. Two-Sided Differences

There is another differencer, two-sided (bilateral), given by

$${\Theta }_0^{1}f\left(t\right)=f(t+h/2)-f(t-h/2),$$

(70)

that originates two bilateral fractional differences, but we will not study both here [65,66,106]. We can easily obtain one of them by combining a delta with a nabla difference.

Definition 5. 

Let, $t=kh,k\in \mathbb{Z}$. We define a bilateral differencer through

$${\Theta }_{\theta }^{\beta }{f}_k={\nabla }^a{\Delta }_{ }^b{f}_k=\sum _{n=-\infty }^{+\infty }\frac{{(-1)}^n\Gamma (\beta +1)}{\Gamma \left(\frac{\beta +\theta }{2}-n+1\right)\Gamma \left(\frac{\beta -\theta }{2}+n+1\right)}{f}_{k-n},$$

(71)

with $\beta =a+b$ (difference order) and $\theta =a-b$ (asymmetry parameter).

For $\beta \in {\mathbb{Z}}^{-},$ we obtain singular cases that were treated in [66]. Suitable choices of these paramters allow us to recover the above introduced differences:

• With $\theta =\beta$ we obtain (47);
• Similarly, with $\theta =-\beta$ and variable change we get (57);

Some particular cases of (71) are very interesting:

1.

Riesz type difference, $\theta =0,$

$${\Theta }_0^{\beta }{f}_k=\sum _{n=-\infty }^{+\infty }\frac{{(-1)}^n\Gamma (\beta +1)}{\Gamma (\frac{\beta }{2}-n+1)\Gamma (\frac{\beta }{2}+n+1)}{f}_{k-n}.$$

(72)

2.

Feller type difference, $\theta =\pm 1,$

$${\Theta }_1^{\beta }{f}_k=\sum _{n=-\infty }^{+\infty }\frac{{(-1)}^n\Gamma (\beta +1)}{\Gamma (\frac{\beta +1}{2}-n+1)\Gamma (\frac{\beta -1}{2}+n+1)}{f}_{k-n}.$$

(73)

3.

two-sided discrete Hilbert transform, $\beta =0,$

$${\Theta }_{\theta }^0{f}_k=\sum _{n=-\infty }^{+\infty }{(-1)}^n·\frac{1}{\Gamma \left(\frac{\theta }{2}-n+1\right)\Gamma \left(\frac{-\theta }{2}+n+1\right)}{f}_{k-n}.$$

(74)

With $\theta =1$ we obtain the usual discrete-time formulation of the Hilbert transform [14]

Remark 5. 

It must be noticed that, in the general case stated in (71), the difference ${\Theta }_{\theta }^{\beta }{f}_k$ does not have Z transform, since the ROC degenerates in the unit circle ($z=e^{i\omega }$), meaning that it has discrete-time Fourier transform, $F\left(e^{i\omega }\right)$ [65,106]:

$$\mathcal{F}\left[{\Theta }_{\theta }^{\beta }{f}_k\right]={|2sin(\omega /2)|}^{\beta }{\mathrm{e}}^{-i\omega \theta }{\mathrm{e}}^{isgn\left(\omega \right)\theta \pi /2}F\left(e^{i\omega }\right),\left|\omega \right|<\pi ,$$

(75)

where $sgn(.)$ is the signum function.

5.6. The Tempered Differences

We introduced above the three standard definitions for the differences of order one (40), (50), and (70). With a slight modification we obtain their tempered versions. We only need to make adaptations of the results introduced in [107].

Definition 6. 

Let $\lambda \in \mathbb{R}$. We can define three tempered differences (TD):

• Nabla TD

  $${\nabla }_{\lambda ,f}f\left(t\right)=f\left(t\right)-e^{-\lambda h}f(t-h),$$

  (76)
  that has LT

  $$\mathcal{L}\left[{\nabla }_{\lambda }f\left(t\right)\right]=\left[1-e^{-\lambda h}{e}^{-sh}\right]F\left(s\right),$$

  (77)
  for any $s\in \mathbb{C}.$
• Delta TD

  $${\Delta }_{\lambda }f\left(t\right)=f\left(t\right)-e^{\lambda h}f(t+h).$$

  (78)
  Its LT is valid for any $s\in \mathbb{C}$ and given by

  $$\mathcal{L}\left[{\Delta }_{\lambda }f\left(t\right)\right]=\left[1-e^{\lambda h}{e}^{sh}\right]F\left(s\right),$$

  (79)
  As above, we removed a (−) sign.
• Two-sidedTD

  $${\Theta }_{\theta ,\lambda }f\left(t\right)=\left[e^{\lambda \frac{h}{2}}f(t+\frac{h}{2})-e^{-\lambda \frac{h}{2}}f(t-\frac{h}{2}).\right]$$

  (80)

It is straightforward to invert the relations (76) and (78), and so we obtain the first order anti-differences

$${\nabla }_{\lambda }^{-1}f\left(t\right)=\sum _{n=0}^{\infty }{e}^{-n\lambda h}f(t-nh)$$

(81)

and

$${\Delta }_{\lambda }^{-1}f\left(t\right)=\sum _{n=0}^{\infty }{e}^{n\lambda h}f(t+nh)$$

(82)

that can be generalized for any real order.

Definition 7. 

For $\alpha \in \mathbb{R}$ we can write [107]

$${\nabla }_{\lambda }^{\alpha }f\left(t\right)=\sum _{n=0}^{\infty }\frac{{(-\alpha )}_n}{n!}{e}^{-n\lambda h}f(t-nh),$$

(83)

and

$${\Delta }_{\lambda }^{\alpha }f\left(t\right)=\sum _{n=0}^{\infty }\frac{{(-\alpha )}_n}{n!}{e}^{n\lambda h}f(t+nh).$$

(84)

Definition 8. 

To obtain the discrete-time versions we set $t=nh$ and $\mu =e^{\lambda h}$ so that we can write [107]

$${\nabla }_{\lambda }^{\alpha }{f}_k=\sum _{n=0}^{\infty }\frac{{(-\alpha )}_n}{n!}{\mu }^{-n}{f}_{k-n},$$

(85)

and

$${\Delta }_{\lambda }^{\alpha }{f}_k=\sum _{n=0}^{\infty }\frac{{(-\alpha )}_n}{n!}{\mu }^n{f}_{k+n},$$

(86)

valid for any $\alpha \in \mathbb{R}$.

The bilateral tempered fractional differences are somehow more involved. They can be obtained from the results introduced in [108]. We will not do it here.

5.7. Bilinear Differences

One of the most interesting methods of obtaining discrete-time systems from continuous-time templates is through the conformal transformations, so that we pass from the LT context, (s), to the one of Z transform, (z), [28]. The traditional procedures can be described by

$$s=1-z^{-1},$$

for the nabla case and

$$s=z-1,$$

for the delta case. Basically, $z^{-1}$ and z correspond to delay and lead, respectively. This procedure is equivalent to the one described in sub-Section 5.4. It has a great drawback: the imaginary axis in the s plane is transformed into the Hilger circles $|z-1|=1$ and $|z+1|=1$ [25]. This fact brings some inconvenients [28].

Another alternative to obtain continuous to discrete conversion is through the bilinear (Tustin) transformation [28,109]

$$s=\frac{2}{h}{\displaystyle \frac{1-z^{-1}}{1+z^{-1}}},$$

which leads to discrete-time formulations similar to the GL-like above obtained. In [110] new discrete-time fractional derivatives based on the bilinear transformation were introduced and study. In agreemennt with such theory we propose new differences here.

Definition 9. 

Let us define bilinear nabla differencer as a linear system whose output, at a given instant, is given by

$$g\left(t\right)+g(t-h)=f\left(t\right)-f(t-h)$$

(87)

and the bilinear delta differencer by

$$g\left(t\right)+g(t+h)=f\left(t\right)-f(t+h).$$

(88)

In terms of the LT, we can write, respectively

$$G\left(s\right)={\displaystyle \frac{1-e^{-sh}}{1+e^{-sh}}}F\left(s\right)$$

(89)

and

$$G\left(s\right)={\displaystyle \frac{1-e^{sh}}{1+e^{sh}}}F\left(s\right).$$

(90)

It is obvious that we can formulate symmetric bilinear differences from the LT

$$G\left(s\right)={\displaystyle \frac{e^{sh/2}-e^{-sh/2}}{e^{sh/2}+e^{-sh/2}}}F\left(s\right).$$

(91)

To avoid repetition of concepts, we will consider the nabla case, only. The above formulae suggest immediately how we make a continuous-to-discrete conversion by using the substitution $e^{sh}\to z$.

Definition 10. 

De define the bilinear nabla fractional difference through

$$\mathcal{Z}\left[{\nabla }_b^{\alpha }{f}_k\right]={\left[{\displaystyle \frac{1-z^{-1}}{1+z^{-1}}}\right]}^{\alpha }F\left(z\right),\left|z\right|>1.$$

(92)

Attending to the results in [110], we can write

$${\nabla }_b^{\alpha }{f}_n=\sum _{k=0}^{\infty }{\psi }_k^{\alpha }{f}_{n-k},n\in \mathbb{Z},$$

(93)

where [110]

$${\psi }_k^{\alpha }=\sum _{m=0}^k\frac{{(-\alpha )}_m}{m!}\frac{{(-1)}^{k-m}{\left(\alpha \right)}_{k-m}}{(k-m)!},k\in {\mathbb{Z}}_0^{+}$$

(94)

In particular, for $\alpha =\pm 1$ we have

$$\begin{array}{c}•{\psi }_k^1=\left\{\begin{array}{cc}0\hfill & k<0\hfill \\ 1\hfill & k=0\hfill \\ 2{(-1)}^k\hfill & k>0\hfill \end{array}\right.\hfill \\ •{\psi }_k^{-1}=\left\{\begin{array}{cc}0\hfill & k<0\hfill \\ 1\hfill & k=0\hfill \\ 2\hfill & k>0\hfill \end{array}\right.\hfill \end{array}$$

Although these new differences seem to be rather involved, they are easily implemented with the help of the fast Fourier transform (FFT) [110].

We can also obtain two-sided bilinear differences, but the procedure is rather involved and not very useful. However, the frequency domain representation is easily obtained from (91).

6. Scale–Invariant Differences

In previous sections, we delt preferably with the shift-invariant differences. Here, we will consider other that have deep relations with the scale-invariant derivatives [36].

Consider two bounded piecewise continuous functions $f\left(t\right),g\left(t\right),t\in {\mathbb{R}}^{+}$ with $f\left(0\right),g\left(0\right)=0.$ For simplicity, assume they are of polynomial order, so that we assure they have Mellin transforms, $F\left(v\right),G\left(v\right),v\in \mathbb{C}$, analytic over suitable ROC.

Definition 11. 

Let $q>1$. We define stretching differencer as a linear system whose output, at a given scale, is the difference between the input at different scales:

$$g\left(\tau \right)=f\left(\tau \right)-f\left(\tau q^{-1}\right),\tau \in {\mathbb{R}}^{+},$$

(95)

where q is the scale constant.

The output will be called stretching difference and represented by $⊲f\left(t\right)$. Letting $q<1$, we obtain the shrinking difference, $⊳f\left(t\right)$. Therefore, their properties are similar. We will study the first only, because the other is easy to obtain. From this definition we can draw some conclusions:

1.

Its impulse response is given by:

$$\varphi \left(\tau \right)=\delta (\tau -1)-\delta (\tau /q-1),$$

so that

$$g\left(\tau \right)=\varphi \left(\tau \right)★f\left(\tau \right),$$

implying that

$$G\left(v\right)=(1-q^{-v})F\left(v\right).$$

(96)

2.

The transfer function is

$$H\left(v\right)=1-q^{-v},$$

(97)

having $\mathbb{C}$ as ROC.

3.

As in the shift-invariant case, if associate in series n systems the resulting system defines the n–th order stretching difference that has a transfer function given by

$$H\left(v\right)={\left[1-q^{-v}\right]}^n,$$

from where we obtain the n–th order stretching difference

$${⊲}^{n}f\left(\tau \right)=\sum _{k=0}^n{(-1)}^k\left(\genfrac{}{}{0pt}{}{n}{k}\right)f\left(\tau q^{-k}\right).$$

(98)

Now, let us return back to (40) and invert the role of the functions: assume that the input is $g\left(\tau \right)$ and the output is $f\left(\tau \right)$

$$f\left(\tau \right)=f\left(\tau q^{-1}\right)+g\left(\tau \right),$$

that can be reused to give:

$$f\left(\tau \right)=g\left(\tau \right)+g\left(\tau q^{-1}\right)+f\left(\tau q^{-2}\right)=g\left(\tau \right)+g\left(\tau q^{-1}\right)+g\left(\tau q^{-2}\right)+f(\tau q^{-3}=\cdots$$

It is not hard to see that

$$f\left(\tau \right)={⊲}^{-1}g\left(\tau \right)=\sum _{n=0}^{\infty }g\left(\tau q^{-n}\right),$$

(99)

with MT

$$F\left(v\right)=G\left(v\right)\sum _{n=0}^{\infty }{q}^{-nv}={\left[1-q^{-v}\right]}^{-1}G\left(s\right),\mathit{for}Re\left(v\right)>0.$$

(100)

Relation (99) shows why this operator is again an accumulator. The series association of n equal accumulators gives the n–th order stretching sum:

$${⊲}^{-n}g\left(\tau \right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{-n}{k}\right)g\left(\tau q^{-k}\right).$$

(101)

Definition 12. 

This result together with (98) suggest that the α–order stretching differencer/accumulator must be given by

$${⊲}^{\alpha }f\left(\tau \right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f\left(\tau q^{-k}\right),$$

(102)

where α is any real number.

As we observe, this difference uses an exponential domain in agreement with our considerations above (Section 2.3). The corresponding MT is

$$\mathcal{M}\left[{⊲}^{\alpha }f\left(\tau \right)\right]={\left[1-q^{-v}\right]}^{\alpha }F\left(v\right),\mathit{for}Re\left(v\right)>0.$$

(103)

Therefore, the shrinking difference is defined by

$${⊳}^{\alpha }f\left(\tau \right)=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f\left(\tau q^k\right),$$

(104)

having MT

$$\mathcal{M}\left[{⊳}^{\alpha }f\left(\tau \right)\right]={\left[1-q^v\right]}^{\alpha }F\left(v\right),\mathit{for}Re\left(v\right)<0.$$

(105)

Remark 6. 

The relations (100) and (105) show that there is a scale-invariant system that produces the difference as output.

To obtain the discrete-scale versions we only need to make a sampling in agreement with theorem 2. Let $\tau =q^n$. For the stretching difference we obtain

$${⊲}^{\alpha }{f}_{q^n}=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f_{q^{n-k}},$$

(106)

while for the shrinkning one it comes

$${⊳}^{\alpha }{f}_{q^n}=\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f_{q^{n+k}}.$$

(107)

As we can see, they are formally similar to the nabla and delta differences; they only differ in the sampling sequence used: linear or exponential. Thus, from a purely discrete point of view, we have no way of making any distinction between linearly or exponentially sampled functions. This means that if we want to define stretching and shrinking differences in discrete time, we have to break the continuous connection: we have to work exclusively with sequences defined in $\mathbb{Z}$. So, having a sequence $f_n,n\in \mathbb{Z}$, we wonder how to define stretched and shrunk sequences so that we can introduce differences. In [111] ways to produce stretched and shrunk sequences were presented. In principle, we could use them to define differences, but this procedure has difficulties, since the operations of stretching and shrinking involve all knowledge of the underlying sequence and the scale transformation system is two-sided.

However, we can use the traditional “decimation” operation used in signal processing to define a stretching difference [5,8,14].

Definition 13. 

Let N be a positive integer (decimation parameter). We define a stretching difference through

$${⊲}_N{f}_n=f_n-f_{Nn}.$$

(108)

We can show immediately that, if M is a positive integer, then

$${⊲}_N^M{f}_n=\sum _{m=0}^M{(-1)}^m\left(\genfrac{}{}{0pt}{}{M}{m}\right)f_{N^{m}n}.$$

(109)

If $f_n=n^a,n\in \mathbb{N},a\in \mathbb{R}$, then

$${⊲}_N^M{n}^a=\sum _{m=0}^M{(-1)}^m\left(\genfrac{}{}{0pt}{}{M}{m}\right)N^{am}{n}^a={\left[1-N^a\right]}^M{n}^a.$$

(110)

Proceeding as above, we obtain

$${⊲}_N^{-M}{f}_n=\sum _{m=0}^M\frac{{\left(M\right)}_m}{m!}{f}_{N^{m}n},$$

(111)

that allows us to write

$${⊲}_N^M{n}^a=\sum _{m=0}^M{(-1)}^m\frac{{\left(M\right)}_m}{m!}{N}^{am}{n}^a={\left[1-N^a\right]}^{-M}{n}^a.$$

(112)

Definition 14. 

These relations suggest we write

$${⊲}_N^{\alpha }{f}_n=\sum _{m=0}^{\infty }\frac{{(-\alpha )}_m}{m!}{f}_{N^{m}n},$$

(113)

so that

$${⊲}_N^{\alpha }{n}^a=\sum _{m=0}^{\infty }\frac{{(-\alpha )}_m}{m!}{N}^{am}{n}^a={\left[1-N^a\right]}^{\alpha }{n}^a.$$

(114)

This last relation seems pointing to a definition of a “discrete Mellin transform,” as

$$F\left(v\right)=\sum _{k=1}^{\infty }{f}_k{k}^v,v\in \mathbb{C},$$

which is different from other proposed in the past [33,34]. We do not go further in this way.

The properties of the stretching discrete difference just proposed are readily deduced from the results in Section 5.3. From the such a definition, it is evident the reason for not defining a shrinking difference.

7. The ARMA Type Difference Linear Systems

Definition 15. 

We define an ARMA type difference linear system through the following equation [25,28]

$$\sum _{k=0}^N{a}_k{\mathcal{D}}^{{\alpha }_k}y\left(t\right)=\sum _{k=0}^M{b}_k{\mathcal{D}}^{{\beta }_k}x\left(t\right),$$

(115)

where $a_k$ and $b_k$ ($k=0,1,\cdots ,$) with $a_N=1$ are real numbers. The operator $\mathcal{D}$ is any difference defined previously. The orders N and M are any positive integers. The positive real numbers ${\alpha }_k$ and ${\beta }_k$ with $k=0,1,\cdots ,$ form strictly increasing sequences.

The most interesting systems are those with commensurate orders:

$$\sum _{k=0}^N{a}_k{\mathcal{D}}^{k\alpha }y\left(t\right)=\sum _{k=0}^M{b}_k{\mathcal{D}}^{k\alpha }x\left(t\right).$$

(116)

In the shift-invariant cases, the exponentials are the eigenfunctions and the eigenvalues are the transfer functions [27,41,112]. In the scale-invariant, such role is played by the powers [36]. The corresponding (Laplace, Z, Mellin) transforms give the impulse response or Green function of the system.

Example 2. 

Consider the ARMA(1,1) system:

$${\nabla }^{\alpha }y\left(t\right)+ay\left(t\right)=b_1{\nabla }^{\alpha }x\left(t\right)+b_{0}x\left(t\right)$$

and let $x\left(t\right)=e^{st},t\in \mathbb{R}$. The output will be $y\left(t\right)=H\left(s\right)e^{st},t\in \mathbb{R},c\in \mathbb{C},$ with

$$H\left(s\right)={\displaystyle \frac{b_1{(1-e^{-sh})}^{\alpha }+b_0}{{(1-e^{-sh})}^{\alpha }+a}},Re\left(s\right)>0.$$

The impulse response is not easily obtained: From,

$$\begin{array}{cc}\hfill H\left(s\right)& ={\displaystyle \frac{b_1{(1-e^{-sh})}^{\alpha }+b_0+b_{1}a-b_{1}a}{{(1-e^{-sh})}^{\alpha }+a}}\hfill \\ & =b_1+{\displaystyle \frac{b_0-b_{1}a}{{(1-e^{-sh})}^{\alpha }+a}}=b_1+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m{(1-e^{-sh})}^{-(m+1)\alpha },\hfill \\ & =b_1+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m\sum _{k=0}^{\infty }\frac{{\left((m+1)\alpha \right)}_k}{k!}{e}^{-ksh}.\hfill \end{array}$$

we get

$$h\left(t\right)=b_1\delta \left(t\right)+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m\sum _{k=0}^{\infty }\frac{{\left((m+1)\alpha \right)}_k}{k!}\delta (t-kh).$$

The output for any $x\left(t\right)$ is given by the usual convolution

$$y\left(t\right)=b_{1}x\left(t\right)+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m\sum _{k=0}^{\infty }\frac{{\left((m+1)\alpha \right)}_k}{k!}x(t-kh).$$

In the discrete-time case, we set $t=nh,n\in \mathbb{Z}$ to give

$${\nabla }^{\alpha }{y}_n+a{y}_n=b_1{\nabla }^{\alpha }{x}_n+b_0{x}_n.$$

If the input is the exponential $z^n,n\in \mathbb{Z},z\in \mathbb{C}$, the output is $y_n=H\left(z\right)z^n,n\in \mathbb{Z},z\in \mathbb{C}$, where

$$H\left(z\right)={\displaystyle \frac{b_1{(1-z^{-1})}^{\alpha }+b_0}{{(1-z^{-1})}^{\alpha }+a}},\left|z\right|>1.$$

We can obtain the Z transform inverse of this function approximately through the FFT. However, proceeding as above, we get

$$h_n=b_1{\delta }_n+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m\sum _{k=0}^{\infty }\frac{{\left((m+1)\alpha \right)}_k}{k!}{\delta }_{n-k}.$$

As in the continuous-time case, the output corresponding to a given input is obtained through the discrete convolution

$$y_n=b_1{x}_n+(b_0-b_{1}a)\sum _{m=0}^{\infty }{a}^m\sum _{k=0}^{\infty }\frac{{\left((m+1)\alpha \right)}_k}{k!}{x}_{n-k}.$$

Example 3. 

Consider the scale-invariant AR(1) difference system defined by the equation

$${⊳}_3^{1/2}{y}_n={⊳}_3^{1/2}{x}_n-x_n$$

and let $x_n=n^{-2},n\in {\mathbb{Z}}_0^{-}$.

Accordingly to what we wrote in the last section, the solution is, for $n\in \mathbb{N}$,

$$y_n=x_n-\sum _{m=0}^{\infty }\frac{{(1/2)}_m}{m!}{x}_{3^{m}n}.$$

So,

$$y_n=n^{-2}-\sum _{m=0}^{\infty }\frac{{(1/2)}_m}{m!}{3}^{-2m}{n}^{-2}=n^{-2}\left[1-{(1-1/9)}^{-1/2}\right]=\frac{2\sqrt{2}-3}{2\sqrt{2}}{n}^{-2}.$$

8. Discussion

Differences are basic building blocks for defining derivatives, but they can be used in many applications, to solve differential equations and model many systems. In most situations, shift-invariant differences are used, although scale-invariant versions are also useful. Here they have been studied separately. General continuous-time cases have been introduced, although the main interest has been placed on the discrete versions. These are fundamental in computational applications. We adopted a system point of view to emphasize that differences are outputs of linear systems, which implies that they are defined independently of the inputs, notably their duration. Moreover, outputs are generally of infinite duration, even if the input support is finite. In addition to the classic differences, we introduce new ones such as bilateral and tempered differences.

The option by system approach for introducing differences allowed us to define ARMA type linear systems enlarging the classic procedure used in time series analysis and processing which supported many important applications in Engineering, Economics, Statistics, and so on. It is important to remark that many interesting functions we find in applications are acquired under a discrete-time form, without having any analytical formula. This implies that we have to deal with functions (signals) defined in the set of integers. Anyway, implicit in any application, there is a time sequence imposed by an underlying clock which imposes a working domain that we cannot change. This aspect was frequently dismissed in the past, a fact that led to some “abnormalities” like the lost of (anti-)causality. This happened, for example, with the assumed “delta difference” that is not really a delta difference, since it should be anti-causal, but it is bilateral. This fact is expected to make a review of some associated concepts and tools.