The Fractional Orthogonal Difference with Applications

Abstract

This paper is a follow-up of a previous paper of the author published in Mathematics journal in 2015, which treats the so-called continuous fractional orthogonal derivative. In this paper, we treat the discrete case using the fractional orthogonal difference. The theory is illustrated with an application of a fractional differentiating filter. In particular, graphs are presented of the absolutel value of the modulus of the frequency response. These make clear that for a good insight into the behavior of a fractional differentiating filter, one has to look for the modulus of its frequency response in a log-log plot, rather than for plots in the time domain.

1. Introduction

In [1], a comprehensive survey is given about the orthogonal derivative. This derivative has a long, but almost forgotten history. The fractional derivative has also a long history. In [2], these two subjects are combined with the fractional orthogonal derivative. This derivative make use of orthogonal polynomials, both continuous (e.g., Jacobi polynomials) and discrete (Hahn polynomials). In this paper, the theory is applied to functions given for a discrete set of points. We use the fractional difference. Grünwald [3] and Letnikov [4] developed a theory of discrete fractional calculus for causal functions (f (x) = 0 for x < 0). We use the definition defined by, e.g., Kuttner [5] or Isaacs [6]. We show that the definition of Grünwald–Letnikov for the fractional difference is a special case of the definition of the fractional orthogonal difference for the Hahn polynomials for causal functions. Bhrawy et al. [7] use orthogonal polynomials to solve fractional differential equations. However, there formulation of the fractional derivative is entirely different from ours.

In Section 2, we give an overview of some results for the fractional derivative following [2] (Section 2).

In Section 3, the theory of the fractional derivative will be applied with a discrete function, and so, the fractional difference is used as a base for deriving a formula for the fractional orthogonal difference.

In Section 4, the theory of the fractional difference is applied with the discrete Hahn polynomials.

In Section 5, we derive the frequency response for the fractional difference with the discrete Hahn polynomials.

In Section 6, the theory will be applied to a fractional differentiating filter. Analog, as well discrete filters are treated. The theory is illustrated by the graphs of the modulus of the frequency responses.

2. Overview of the Results of the Fractional Derivative

In [2], we find the following formulas for the fractional integral: the Weyl and the Riemann–Liouville transform. In terms of these, we can define a fractional derivative.

Let µ ∈ ℂ with $ℛ$e(µ) > 0. Let f be a function on (−∞,b), which is integrable on bounded subintervals. Then, the Riemann–Liouville integral of order µ is defined as:

$$R^{-\mu }\left[f\right](x)=f^{(-\mu )}(x):=\frac{1}{\mathrm{\Gamma }(\mu )}{\displaystyle {\int }_{-\infty }^{x}f(y){(x-y)}^{\mu -1}dy}$$

(1)

A sufficient condition for the absolute convergence of this integral is f (−x) = O (x^−µ−ϵ),ϵ>0, x→∞

Similarly, for Re(µ) > 0 and f, a locally-integrable function on (a, ∞), such that f (x) = O (x^−µ−ϵ), ϵ>0, x→∞, the Weyl integral of order µ is defined as:

$$W^{-\mu }[f](x)=f^{(-\mu )}(x):=\frac{1}{\mathrm{\Gamma }(x)}{\displaystyle {\int }_x^{\infty }f(y){(y-x)}^{\mu -1}}dy$$

Clearly:

$$R^{-\mu }\left[f({-}_{·})\right](-x)=W^{-\mu }\left[f\right](x)$$

The Riemann–Liouville integral is often given for a function f on [0, b) in the form:

$$R^{-\mu }\left[f\right](x)=f^{(-\mu )}(x):=\frac{1}{\mathrm{\Gamma }(\mu )}{\displaystyle {\int }_0^{x}f(y){(x-y)}^{\mu -1}}dy(x>0)$$

This can be obtained from Equation (1) if we assume that f (x) = 0 for x ≤ 0 (then, f is called a causal function).

The fractional derivatives are defined by the following formulas with Re(ν) < n,n a positive integer and $D=\frac{d}{dx}$:

$$R^v[f](x)=D^n\left[R^{v-n}[f]\right](x)$$

(2)

$$W^v[f](x)={(-1)}^n{D}^n\left[W^{v-n}[f]\right](x)$$

(3)

Put W⁰ = id = R⁰. Under the assumption of sufficient differentiability and convergence, we have W^μW^V = W^μ+v, R^μR^v = R^μ+v for all µ, ν ∈ ℂ and Rⁿ = Dⁿ, Wⁿ = (−1)ⁿ Dⁿ.

For Re(ν) < 0, formula Equations (2) and (3) are easily derived, and for Re(ν) > 0, they are definitions. This is under the condition that all derivatives of order less then n + 1 of the function f (x) should exist.

Because we need the Fourier transform of Formulas (2) and (3), we use the following theorem. The proof is in [8] (Chapter 2, Section 7) for the Riemann–Liouville derivative, as well as for the Weyl derivative.

Theorem 1. Let f and g be functions on ℝ for which the Fourier transforms exist and are given by:

$$F(w)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{iwx}f(x)dx}$$

$$G(w)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{iwx}g(x)dx}$$

If g = W^v [f], then:

$$\frac{G(w)}{F(w)}={(iw)}^v$$

Here and in the rest of this paper, z^v = e^vInz with − π < arg (z) < π. Therefore, (−iw)^v =(−iw+0)^v = e^−πiv/2(w+i0)^v. If z is on the cut (−∞, 0] and z ≠ 0, then we distinguish between (z+i0)^v=e^iπv(−z)^v and (z−i0)^v=e^−iπv(−z)^v

The quotient $H(w)=\frac{G(w)}{F(w)}$ will be called the frequency response, where we follow the terminology of filter theory.

In [1], the following theorem is proven:

Theorem 2. Let n be a positive integer. Let p_n be an orthogonal polynomial of degree n with respect to the orthogonality measure µ for which all moments exist. Let x ∈ ℝ. Let I be a closed interval, such that, for some ε > 0, x + δu ∈ I, if 0 ≤ δ ≤ ε and u ∈ supp(µ). Let f be a continuous function on I, such that its derivatives of order 1, 2,…, n at x exist. In addition, if I is unbounded, assume that f is of at most polynomial growth on I. Then:

$$f^{(n)}(x)=\underset{\delta \downarrow 0}{\lim }{D}_{\delta }^n[f](x)$$

where:

$$D_{\delta }^n[f](x)=\frac{k_{n}n!}{h_n}\frac{1}{{\delta }^n}{\displaystyle {\int }_{ℝ}f(x+\delta \xi )p_n(\xi )d\mu (\xi )}$$

From this formula, with f absolutely integrable on ℝ and $g=D_{\delta }^n[f]$, we immediately derive for the frequency response $H_{\delta }^n(w)$:

$$H_{\delta }^n(w)=\frac{G(w)}{F(w)}=\frac{k_{n}n!}{h_n{\delta }^n}{\displaystyle {\int }_{ℝ}{p}_n(\xi )e^{-iw\delta \xi }d\mu (\xi )}$$

(4)

3. The Fractional Orthogonal Difference

In this section, we treat the fractional difference. If we want to compute the fractional derivative of a function given at a set of equidistant points, we can use instead of the fractional derivative, the fractional difference. Then, we have to work in the formula for the orthogonal derivative with discrete orthogonal polynomials.

As an analog of the Riemann–Liouville fractional integral Equation (1), we can define a fractional summation (see, for example, [5] or [6]) as:

$$I_{\delta }^{-\mu }[f](x):={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}f(x-k\delta )\delta >0}$$

(5)

In contrast to the Riemann–Liouville fractional integral, this summation has no a priori singularity for µ < 0. If f is a causal function, then the upper bound for the summation should be [x/δ], and the sum is certainly well-defined. For µ ≠ 0,− 1,− 2,…, we can write:

$$\frac{{(\mu )}_k}{k!}=\frac{1}{\mathrm{\Gamma }(\mu )}\frac{\mathrm{\Gamma }(\mu +k)}{\mathrm{\Gamma }(1+k)}\sim \frac{k^{\mu -1}}{\mathrm{\Gamma }(\mu )}$$

as k → ∞. Hence, the infinite sum converges if $f(x)=O\left({\left|x\right|}^{-\lambda }\right)$ as x → ∞ with λ > Re(µ). For µ = −n(n a nonnegative integer), we get:

$$I_{\delta }^n[f](x)={\delta }^{-n}{\displaystyle \sum _{k=0}^{\infty }{(-1)}^k}\left({}_k^n\right)f(x-k\delta )={(D_{\delta }^{-})}^n[f](x)$$

where:

$$D_{\delta }^{-}[f](x):=\frac{f(x)-f(x-\delta )}{\delta }$$

(6)

Note that:

$$D_{\delta }^{-}\left[I_{\delta }^{-\mu }[f]\right](x)=I_{\delta }^{-\mu +1}[f](x)$$

More generally, we have:

$$I_{\delta }^{-\mu }\left[I_{\delta }^{-v}[f]\right](x)=I_{\delta }^{-(\mu +v)}[f](x)$$

For µ > 0, the fractional summation $I_{\delta }^{-\mu }[f]$ formally approximates the Riemann–Liouville integral R−^µ [f]. Indeed, putting y = kδ, we can rewrite the formula for $I_{\delta }^{-\mu }[f](x)$ as:

$$I_{\delta }^{-\mu }[f](x)=\frac{{\delta }^{\mu }}{\mathrm{\Gamma }(\mu )}{\displaystyle \sum _{y\in \delta \mathbb{Z}\ge 0}\frac{\mathrm{\Gamma }\left(\mu +\frac{y}{\delta }\right)}{\mathrm{\Gamma }\left(1+\frac{y}{\delta }\right)}f(x-y)=\frac{\delta }{\mathrm{\Gamma }(\mu )}{\displaystyle \sum _{y\in \delta \mathbb{Z}\ge 0}\frac{y^{\mu -1}\mathrm{\Gamma }\left(\mu +\frac{y}{\delta }\right)}{{\left(\frac{y}{\delta }\right)}^{\mu -1}\mathrm{\Gamma }\left(1+\frac{y}{\delta }\right)}}}f(x-y)$$

(7)

The work in [9] (5.11.12) gives, if δ ↓ 0:

$$\underset{\delta \downarrow 0}{\lim }\frac{\mathrm{\Gamma }\left(\frac{y}{\delta }+\mu \right)}{\mathrm{\Gamma }\left(\frac{y}{\delta }+1\right){\left(\frac{y}{\delta }\right)}^{\mu -1}}=1$$

(8)

Formally, we have:

$$\underset{\delta \downarrow 0}{\lim }\delta {\displaystyle \sum _{y\in \delta \mathbb{Z}\ge 0}{g}_{\delta }(y)={\displaystyle {\int }_0^{\infty }g(y)dy}}$$

(9)

with g_δ (y) → g (y) as δ ↓ 0. Application of Equations (8) and (9) to Equation (7) yields:

$$\underset{\delta \downarrow 0}{\lim }{I}_{\delta }^{-\mu }[f](x)=\frac{1}{\mathrm{\Gamma }(\mu )}{\displaystyle {\int }_0^{\infty }{y}^{\mu -1}f(x-y)dy}$$

(10)

This is the fractional integral of Riemann–Liouville. Therefore, it is shown that the fractional difference formally tends to the fractional integral in the limit as δ ↓ 0.

Remark 1. Grünwald [3] and Letnikov [4] developed a theory of discrete fractional calculus for causal functions (f (x) = 0 for x < 0). For a good description of their theory, see [10]. They started with the backwards difference:

$${\mathrm{\Delta }}^{-}{}_{\delta }[f](x)=\frac{f(x)-f(x-\delta )}{\delta }$$

and next, they define the fractional derivative as:

$$\frac{d^v}{d{x}^v}f(x):=\frac{1}{\mathrm{\Gamma }(-v)}\underset{\delta \to 0,N\to \infty }{\lim }\frac{1}{{\delta }^v}{\displaystyle \sum _{m=0}^{N-1}\frac{\mathrm{\Gamma }(m-v)}{\mathrm{\Gamma }(m+1)}f(x-k\delta )}$$

(11)

where$N=\left[\frac{x}{\delta }\right]$. This definition has the advantage that ν is fully arbitrary. No use is made of derivatives or integrals. The disadvantage is that computation of the limit is often very difficult. However, it can be shown that Equation (11) is the same as the Riemann–Liouville derivative for causal functions for Re(ν) > 0 and f is n times differentiable with n > Re(ν). At least formally, this follows also our approach if we write$I_{\delta }^v[f](x)={(D_{\delta }^{-})}^n\left[I_{\delta }^{v-n}[f]\right](x)$ and then use Equation (10) and the fact that Equation (6) approaches the derivative as δ ↓ 0.

Comparing Equation (11) with Equation (5), the differences are the upper bound and the limit. For practical calculations, Formula (11) is often used without the limit. In the next section, we will develop the fractional orthogonal difference for the Hahn polynomials. In Section 5, we will see that the definition of Grünwald–Letnikov for the fractional difference is a special case of the definition of the fractional orthogonal difference for the Hahn polynomials for causal functions.

For obtaining the formula for the fractional orthogonal difference, we let ${\left\{p_n\right\}}_{n=0}$ be a system of orthogonal polynomials with respect to weights w (j) on points j (j = 0, 1,…, N). Then, the approximate orthogonal derivative [2] (3.1b) takes the form:

$$D_{\delta }^n[f](x)=\frac{k_{n}n!}{h_n}\frac{1}{{\delta }^n}{\displaystyle \sum _{j=0}^{N}f(x+j\delta )w(j)p_n(j)n=0,1,\dots ,N}$$

(12)

Let Re(ν) < n with n a positive integer. The fractional derivative $R^v[f](x)=D^n\left[R^{v-n}[f]\right](x)$ can be formally approximated by:

$$I_{\delta ,n,w}^v[f](x):=D_{\delta }^n\left[I_{\delta }^{v-n}[f]\right](x)$$

(13)

as δ ↓ 0. This is a motivation to compute $I_{\delta ,n,w}^v[f]$ more explicitly, in particular for special choices of the weights w. Substitution of Equations (5) and (12) in Equation (13) gives:

$$I_{\delta ,n.w}^v[f](x)=\frac{k_{n}n!}{h_n}\frac{1}{{\delta }^v}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{j=0}^N\frac{{(n-v)}_k}{k!}f(x+(j-k)\delta )w(j)p_n(j)}}$$

(14)

For the double sum, we can write:

$$\begin{array}{r}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{j=0}^{N}g(k,j)={\displaystyle \sum _{j=0}^N{\displaystyle \sum _{k=0}^{\infty }g(k,j)={\displaystyle \sum _{j=0}^N{\displaystyle \sum _{m=-j}^{\infty }g(j+m,j)=}}}}}}\\ ={\displaystyle \sum _{j=0}^N{\displaystyle \sum _{m=1}^{\infty }g(j+m,j)+{\displaystyle \sum _{j=0}^N{\displaystyle \sum _{m=0}^{j}g(j-m,j)=}}}}\\ ={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{j=0}^{N}g(j+m,j)+{\displaystyle \sum _{m=0}^N{\displaystyle \sum _{j=m}^{N}g(j-m,j)}}}}\end{array}$$

Application of this formula to Equation (14) yields:

$$\begin{array}{r}{I}_{\delta ,n,w}^v[f]=\frac{k_n}{h_n}\mathrm{\Gamma }(n+1)\frac{1}{{\delta }^v}{\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )\left[{\displaystyle \sum _{j=0}^N\frac{{(n-v)}_{j+m}}{(j+m)!}{p}_n(j)w(j)}\right]+}\\ +\frac{k_n}{h_n}\mathrm{\Gamma }(n+1)\frac{1}{{\delta }^v}{\displaystyle \sum _{m=0}^{N}f(x+m\delta )\left[{\displaystyle \sum _{j=m}^N\frac{{(n-v)}_{j-m}}{(j-m)!}{p}_n(j)w(j)}\right]}\end{array}$$

The summations inside the square brackets are in principle known after a choice of the discrete orthogonal polynomials p_n (x). Then, this formula can be used for approximating the fractional orthogonal derivative.

Remark 2. When taking the Fourier transform of the fractional summation, one can show that for δ ↓ 0, this Fourier transform formally tends to the Fourier transform of the fractional integral. For this purpose, we write Equation (5) as:

$$I_{\delta }^{-\mu }[f](n\delta )={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}}f((n-k)\delta )$$

Taking the discrete Fourier transform gives:

$${\displaystyle \sum _{n=-\infty }^{\infty }{I}_{\delta }^{-\mu }[f](n\delta )e^{in\delta w}={\delta }^{\mu }}{\displaystyle \sum _{n=-\infty }^{\infty }{e}^{in\delta w}}{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}f((n-k)\delta )}$$

while working formally, we interchange the summations on the right. This gives:

$$\begin{array}{l}{\displaystyle \sum _{n=-\infty }^{\infty }{I}_{\delta }^{-\mu }[f](n\delta )e^{in\delta w}={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}{e}^{ik\delta w}{\displaystyle \sum _{n=-\infty }^{\infty }f((n-k)\delta )e^{i(n-k)\delta w}=}}}\\ ={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}{e}^{ik\delta w}}{\displaystyle \sum _{n=-\infty }^{\infty }f(n\delta )e^{in\delta w}}\end{array}$$

The first summation is well known. We obtain:

$$\begin{array}{l}{{\displaystyle \sum _{n=-\infty }^{\infty }{I}_{\delta }^{-\mu }[f](n\delta )e^{in\delta w}=\left(\frac{1-e^{i\delta w}}{\delta }\right)}}^{-\mu }{\displaystyle \sum _{n=-\infty }^{\infty }f(n\delta )e^{in\delta w}}=\\ =w^{-\mu }{\left(\frac{\sin (\delta w/2)}{\delta w/2}\right)}^{-\mu }{e}^{-i\mu (\delta w-\pi )/2}{\displaystyle \sum _{n=-\infty }^{\infty }f(n\delta )e^{in\delta w}}\end{array}$$

Multiply both sides with δ, and take the formal limit for δ ↓ 0. We obtain the well-known formula:

$${\displaystyle {\int }_{-\infty }^{\infty }{R}^{-\mu }[f](x)e^{ixw}dx=e^{i\pi \mu /2}{w}^{-\mu }{\displaystyle {\int }_{-\infty }^{\infty }f(x)e^{ixw}}dx}$$

4. The Fractional Orthogonal Difference for the Discrete Hahn Polynomials

In this section, we substitute the discrete Hahn polynomials Q_n(x;α,β,N) in Formula (14) to calculate the fractional difference. We do this analogously to Jacobi polynomials in [2]. For the Hahn polynomials, we use the definition in [11]. They give:

$$p_n(j)=Q_n(j;\alpha ,\beta ,N):{=}_3{F}_2\left(\begin{array}{cc}\begin{array}{c}-j,-n,n+\alpha +\beta +1\\ \alpha +1,-N\end{array}& ;1\end{array}\right)$$

(15)

For the Hahn polynomials, there are the following properties with n = 0,1, ,…, N:

$$\begin{array}{c}w(j;\alpha ,\beta ,N)=\frac{{(\alpha +1)}_j{(B+1)}_{N-j}}{j!(N-j)!}\\ \frac{k_n}{h_n}={(-1)}^n\frac{(2_n+\alpha +\beta +1)}{{(\beta +1)}_n}\frac{{(n+\alpha +\beta +1)}_n}{{(n+\alpha +\beta +1)}_{N+1}}\frac{N!}{n!}\end{array}$$

(16)

The work in [11] (9.5.10) gives:

$$Q_n(x;\alpha ,\beta ,N)w(x)=\frac{{(-1)}^n(\beta +1)}{{(-N)}_n}{\nabla }_x^n\left(\frac{{(\alpha +n+1)}_x}{x!}\frac{{(\beta +n+1)}_{N-n-x}}{(N-n-x)!}\right)$$

(17)

where:

$${\nabla }_x(f(x)):=f(x)-f(x-1)$$

Substitution in Equation (14) gives:

$$I_{\delta ,n,w}^v[f]=H_1(n)\left[{\displaystyle \sum _{m-1}^{\infty }f(x-m\delta )J_1(m)+{\displaystyle \sum _{m=0}^{N}f(x+m\delta )J_2(m)}}\right]$$

(18)

with:

$$J_1(m)={\displaystyle \sum _{j=0}^N\frac{{(n-v)}_{j+m}}{(j+m)!}(w{Q}_n)(j;\alpha ,\beta ,N)}$$

(19)

$$J_2(m)={\displaystyle \sum _{j=m}^N\frac{{(n-v)}_{j-m}}{(j-m)!}}(w{Q}_n)(j;\alpha ,\beta ,N)={\displaystyle \sum _{j=0}^{N-m}\frac{{(n-v)}_j}{j!}}(w{Q}_n)(j+m;\alpha ,\beta ,N)$$

(20)

$$H_1(n)={(-1)}^n\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)\mathrm{\Gamma }(\beta +1)\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+n+\alpha +\beta +2)\mathrm{\Gamma }(n+\beta +1)}\frac{1}{{\delta }^v}$$

(21)

To compute the sums in Equations (19) and (20), we use partial summation. The general formula for partial summation is:

$${\displaystyle \sum _{j=0}^N(\nabla f)(j)g(j)=f(N)g(N)-f(-1)g(0)-{\displaystyle \sum _{j=0}^{N-1}f(j)(\nabla g)(j+1)}}$$

(22)

In particular, if f (N) = f (−1) = 0, then:

$${\displaystyle \sum _{i=0}^N(\nabla f)(j)g(j)=-{\displaystyle \sum _{i=0}^{N-1}f(j)(\nabla g)(j+1)}}$$

Applying this formula n times and remarking that f (N − n) = f (−1) = 0 gives:

$${\displaystyle \sum _{j=0}^N{(\nabla f)}^n(j)g(j)={(-1)}^n{\displaystyle \sum _{j=0}^{N-n}f(j){(\nabla g)}^n(j+n)}}$$

(23)

Applying Formula (17) for the summation in Equation (19), we get:

$$J_1(m)=\frac{{(-1)}^n{(\beta +1)}_n}{{(-N)}_n}{\displaystyle \sum _{j=0}^N\frac{{(n-v)}_{j+m}}{(j+m)!}}{\nabla }_j^n\left(\frac{{(\alpha +n+1)}_j}{j!}\frac{{(\beta +n+1)}_{N-n-j}}{(N-n-j)!}\right)$$

Suppose:

$$\begin{array}{l}f(j):=\frac{{(\alpha +n+1)}_j}{j!}\frac{{(\beta +n+1)}_{N-n-j}}{(N-n-j)!}\\ g(j):=\frac{{(n-v)}_{j+m}}{(j+m)!}=\frac{{(j+m+1)}_{n-v-1}}{\mathrm{\Gamma }(n-v)}\end{array}$$

and seeing that f (N − n) = f (−1) = 0, we use Equation (23) to get:

$$J_1(m)=\frac{{(\beta +1)}_n}{\mathrm{\Gamma }(n-v)(-N)}{\displaystyle \sum _{j=0}^{N-n}\frac{{(\alpha +n+1)}_j}{j!}}\frac{{(\beta +n+1)}_{N-n-j}}{(N-n-j)!}{\nabla }_j^n({(j+m+n+1)}_{n-v-1})$$

For the second fraction inside the summation, we use [12] (Appendix I.(30)) ):

$$\mathrm{\Gamma }(a-j)=\frac{\mathrm{\Gamma }(a){(-1)}^j}{{(1-a)}_j}$$

Then, it follows:

$$\frac{{(\beta +n+1)}_{N-n-j}}{(N-n-j)!}=\frac{{(\beta +n+1)}_{N-n}}{(N-n)!}\frac{{(n-N)}_j}{{(-\beta -N)}_j}$$

For the third fraction inside the summation, we use the formula:

$${\nabla }^n({(j+a)}_b)={(a)}_{b-n}{(b-n+1)}_n\frac{{(a+b+n)}_j}{{(a)}_j}$$

(24)

Then, we obtain:

$${\nabla }_j^n({(j+m+n+1)}_{n-v-1})={(-v)}_n\frac{\mathrm{\Gamma }(m+n-v)}{\mathrm{\Gamma }(m+n+1)}\frac{{(m+n-v)}_j}{{(m+n+1)}_j}$$

After substitution in J₁ (m), we get:

$$J_1(m)=\frac{1}{\mathrm{\Gamma }(-v){(-N)}_n}\frac{{(1+\beta )}_N}{(N-n)!}\frac{\mathrm{\Gamma }(m+n-v)}{\mathrm{\Gamma }(m+n+1)}{\displaystyle \sum _{j=0}^{N-n}\frac{{(\alpha +n+1)}_j{(m+n-v)}_j{(n-N)}_j}{{(m+n+1)}_j{(-\beta -N)}_j}}\frac{1}{j!}$$

Then, J₁ (m) can be written as a hypergeometric function:

$$J_1(m)=\frac{1}{{(-N)}_n}\frac{{(1+\beta )}_N}{(N-n)!}\frac{{(-v)}_{m+n}}{(m+n)!}{}_3{F}_2\left(\begin{array}{c}n-N,n+\alpha +1,m+n-v\\ -N-\beta ,m+n+1\end{array};1\right)$$

(25)

Observing the hypergeometric function, we will see that this function represents a series with a finite number of terms (if β is not an integer with − N − β ≠ 0). With α = β = 0 and Re (α) < m, it follows:

$$J_1(m)=\frac{{(N+1-n)}_n}{{(-N)}_n}\frac{{(-v)}_{m+n}}{(m+n)!}{}_3{F}_2\left(\begin{array}{c}n-N,n+1,m+n-v\\ -N,m+n+1\end{array};1\right)$$

Observing the hypergeometric function, we will see that this function represents an infinite series. To compute J₂ (m), we write Formula (20) with Equation (17) as:

$$J_2(m)=\frac{{(-1)}^n{(\beta +1)}_n}{{(-N)}_n}{\displaystyle \sum _{j=0}^{N-m}\frac{{(n-v)}_j}{j!}}{\nabla }_j^n\left(\frac{{(\alpha +n+1)}_{j+m}}{(j+m)!}\frac{{(\beta +n+1)}_{N-n-m-j}}{(N-n-m-j)!}\right)$$

(26)

Suppose:

$$\begin{array}{c}f(j)=\frac{{(\alpha +n+1)}_{j+m}}{(j+m)!}\frac{{(\beta +n+1)}_{N-n-m-j}}{(N-n-m-j)!}\\ g(j)=\frac{{(n-v)}_j}{j!}\end{array}$$

Then, we rewrite the formula for the partial summation Equation (22) as:

$${\displaystyle \sum _{j=-1}^{N-m}(\nabla f)(j)g(j)=f(N-m)g(N-m)-f(-2)g(-1)-{\displaystyle \sum _{j=-1}^{N-m-1}f(j){\nabla }_j^n(j+1)}}$$

In particular, if f (N — m) = g (−1) = 0, then:

$${\displaystyle \sum _{j=0}^{N-m}(\nabla f)(j)g(j)=-{\displaystyle \sum _{j=-1}^{N-m-1}f(j)(\nabla g)(j+1)=-{\displaystyle \sum _{j=0}^{N-m}f(j-1)(\nabla g)(j)}}}$$

Applying this formula n times gives:

$${\displaystyle \sum _{j=0}^{N-m}{(\nabla f)}^n(j)g(j)={(-1)}^n{\displaystyle \sum _{j=0}^{N-m}f(j-n){(\nabla g)}^n(j)}}$$

Substitution in Equation (26) gives:

$$J_2(m)=\frac{{(\beta +1)}_n}{\mathrm{\Gamma }(n-v){(-N)}_n}{\displaystyle \sum _{j=0}^{N-m}\frac{{(\alpha +n+1)}_{j+m-n}}{(j-n+m)}}\frac{{(\beta +n+1)}_{N-m-j}}{(N-m-j)!}{\nabla }^n({(j+1)}_{n-v-1})$$

For the ∇ function, we use Formula (24). We obtain:

$${\nabla }^n({(j+1)}_{n-v-1})=\mathrm{\Gamma }(n-v)\frac{{(-v)}_j}{j!}$$

After substitution in J₂ (m), we get:

$$J_2(m)=\frac{{(\beta +1)}_n}{\mathrm{\Gamma }(n-\nu ){(-N)}_n}{\displaystyle \sum _{j=0}^{N-m}\frac{{(\alpha +n+1)}_{j+m-n}}{(j-n+m)!}}\frac{{(\beta +n+1)}_{N-m-j}}{(N-m-j)!}\frac{{(-\nu )}_j}{j!}$$

For j = 0, the factor (m − n + 1) could be less than or equal to zero. Then, the terms concerning this condition should not be taken. In this case, the summation should be reversed. We set j = N − m − i. This gives:

$$J_2(m)=\frac{{(\beta +1)}_n}{\mathrm{\Gamma }(n-\nu ){(-N)}_n}{\displaystyle \sum _{i=0}^{\min (N-m,N-n)}\frac{{(\alpha +n+1)}_{N-n-i}}{(N-n+i)!}}\frac{{(\beta +n+1)}_i}{i!}\frac{{(-\nu )}_{N-m-i}}{(N-m-j)!}$$

After some manipulations with the Pochhammer symbols, J₂ (m) can be written as a hypergeometric function.

$$J_2(m)=\frac{{(\beta +1)}_n}{{(-N)}_n}\frac{{(\alpha +n+1)}_{N-n}}{(N-n)!}\frac{{(-\nu )}_{N-m}}{(N-m)!}F3{}_2\left(\begin{array}{c}n-N,\beta +n+1,-N+m\\ -\alpha -N,\nu -N+m+1\end{array};1\right)$$

(27)

With α = β = 0, it follows again, after some manipulations of the Pochhammer symbols:

$$J_2(m)=\frac{{(-1)}^n}{\mathrm{\Gamma }(-\nu )}\frac{\mathrm{\Gamma }(N-m-\nu )}{(N-m)!}F3{}_2\left(\begin{array}{c}n-N,n+1,-N+m\\ -N,\nu -N+m+1\end{array};1\right)$$

It can be shown that the hypergeometric function represents a finite series. Then, with [9] (16.4.11), we can rewrite the hypergeometric function. This results in:

$$J_2(m)=\frac{\mathrm{\Gamma }(\nu -2n)}{\mathrm{\Gamma }(\nu +1-n)\mathrm{\Gamma }(-\nu )}\frac{\mathrm{\Gamma }(N-m-\nu +n+1)}{(N-m)!}F3{}_2\left(\begin{array}{c}-n,n+1,-m\\ -N,\nu +1-n\end{array};1\right)$$

For the approximate fractional difference, it results from Equations (18)–(21), (25) and (27):

$$\begin{array}{l}{I}_{\delta ,n}^{\nu }[f]=\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{\mathrm{\Gamma }(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(N+1+\beta )}{\mathrm{\Gamma }(n+1+\beta )}\frac{1}{{\delta }^{\nu }}\\ {\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )\frac{{(-\nu )}_{m+n}}{(m+n)!}}F3{}_2\left(\begin{array}{c}n-N,n+\alpha ,m+n-\nu \\ -N-\beta ,m+n+1\end{array};1\right)+\\ +\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{\mathrm{\Gamma }(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(\alpha +N+1)}{\mathrm{\Gamma }(\alpha +n+1)}\frac{1}{{\delta }^{\nu }}\\ {\displaystyle \sum _{m=0}^{N}f(x+m\delta )\frac{{(-\nu )}_{N-m}}{(N-m)!}F3{}_2\left(\begin{array}{c}n-N,\beta +n+1,-N+m\\ -\alpha -N,\nu -N,m+1\end{array};1\right)}\end{array}$$

Remark 3. To see what conditions are necessary for the given function f (x), so that the first summation is convergent, we split up this summation into two summations:

$$\begin{array}{l}\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(N+1+\beta )}{\mathrm{\Gamma }(N+1+\beta )}\\ \frac{1}{{\delta }^v}{\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )\frac{{(-\nu )}_{m+n}}{(m+n)!}}F3{}_2\left(\begin{array}{c}n-N,n+\alpha +1,m+n-\nu \\ -N-\beta ,m+n+1\end{array};1\right)=\\ =\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(N+1+\beta )}{\mathrm{\Gamma }(n+1+\beta )}\frac{1}{{\delta }^{\nu }}\\ \frac{1}{{\delta }^{\nu }}{\displaystyle \sum _{m=1}^{M}f(x-m\delta )\frac{{(-\nu )}_{m+n}}{(m+n)!}}F3{}_2\left(\begin{array}{c}n-N,n+\alpha +1,m+n-\nu \\ -N-\beta ,m+n+1\end{array};1\right)+\frac{1}{{\delta }^{\nu }}R(x)\end{array}$$

with:

$$\begin{array}{l}R(x)=\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(N+1+\beta )}{\mathrm{\Gamma }(n+1+\beta )}\\ \frac{1}{{\delta }^v}{\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )\frac{{(-\nu )}_{m+n}}{(m+n)!}}F3{}_2\left(\begin{array}{c}n-N,n+\alpha +1,m+n-\nu \\ -N-\beta ,m+n+1\end{array};1\right)\end{array}$$

If M goes to infinity, then for R (x), there remains (with [9] (5.11.12)):

$$\begin{array}{l}R(x)=\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{\mathrm{\Gamma }(2n+\alpha +\beta +2)}{(N+n+\alpha +\beta +2)}\frac{\mathrm{\Gamma }(N+1+\beta )}{\mathrm{\Gamma }(n+1+\beta )}\\ \frac{1}{{\delta }^v}F\left(\begin{array}{c}n-N,n+\alpha +1\\ -N-\beta \end{array};1\right)\frac{1}{\mathrm{\Gamma }(-\nu )}{\displaystyle \sum _{m=M+1}^{\infty }f(x-m\delta )\frac{1}{m^{\nu +1}}}\end{array}$$

The hypergeometric function can be summed. There remains:

$$\begin{array}{l}R(x)=\frac{\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(N+1)}\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{1}{{\delta }^{\nu }}{\displaystyle \sum _{m=M+1}^{\infty }f(x-m\delta )\frac{1}{m^{\nu +1}}}\\ \end{array}$$

(28)

so for convergence of approximate fractional difference, the summation in Equation (28) should converge.

For α = β = 0, there remains:

$$\begin{array}{l}{I}_{\delta ,n}^{\nu }[f]=\frac{\mathrm{\Gamma }(2n+2)}{\mathrm{\Gamma }(N+n+2)}\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{1}{{\delta }^{\nu }}\\ {\displaystyle \sum _{m=M+1}^{\infty }f(x-m\delta )\frac{\mathrm{\Gamma }(m+n-\nu )}{(m+n)}}F3{}_2\left(\begin{array}{c}n-N,n+1,m+n-\nu \\ -N,m+n+1\end{array};1\right)+\\ +{(-1)}^n\frac{\mathrm{\Gamma }(2n+2)}{\mathrm{\Gamma }(N+n+2)}\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{\mathrm{\Gamma }(\nu -2n)}{\mathrm{\Gamma }(\nu +1-n)}\frac{1}{{\delta }^{\nu }}\\ {\displaystyle \sum _{m=0}^{N}f(x+m}\delta )\frac{\mathrm{\Gamma }(N-m-\nu +n+1)}{(N-m)!}F3{}_2\left(\begin{array}{c}-n,n+1,-m\\ -N,\nu +1-n\end{array};1\right)\end{array}$$

For n = 1, it follows:

$$\begin{array}{l}{I}_{\delta ,n}^{\nu }[f]=\frac{6}{\mathrm{\Gamma }(N+3)}\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{1}{{\delta }^{\nu }}{\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )\frac{\mathrm{\Gamma }(m+1-\nu )}{\mathrm{\Gamma }(m+2)}}F3{}_2\left(\begin{array}{c}n-N,2,m+1,-\nu \\ -N,m+2\end{array};1\right)-\\ -\frac{6}{\mathrm{\Gamma }(N+3)}\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{\mathrm{\Gamma }(\nu -2)}{\mathrm{\Gamma }(\nu )}\frac{1}{{\delta }^{\nu }}{\displaystyle \sum _{m=0}^{\infty }f(x-m\delta )\frac{\mathrm{\Gamma }(N-m-\nu +2)}{\mathrm{\Gamma }(N-m+1)}}F3{}_2\left(\begin{array}{c}-1,2,-m\\ -N,\nu \end{array};1\right)\\ \end{array}$$

(29)

The first hypergeometric function can be written as:

$$F3{}_2\left(\begin{array}{c}1-N,2,m+1-\nu \\ -N,m+2\end{array};1\right)=\frac{\mathrm{\Gamma }(m+2)}{(-N)\mathrm{\Gamma }(m+1-\nu )}{\displaystyle \sum _{k=0}^{N-1}(k-N)(k+1)\frac{\mathrm{\Gamma }(k+m+1-\nu )}{\mathrm{\Gamma }(k+m+2)}}$$

Making use of:

$$(k-N)(k+1)=(k+m+1)(k+m)-(N+2m)(k+m+1)+m(N+m+1)$$

and of [13] (2.5(16)):

$${\displaystyle \sum _{k=0}^{N-1}\frac{\mathrm{\Gamma }(k+a)}{\mathrm{\Gamma }(k+b)}}=\frac{1}{(a-b+1)}\left[\frac{\mathrm{\Gamma }(N+a)}{\mathrm{\Gamma }(N+b-1)}-\frac{\mathrm{\Gamma }(a)}{\mathrm{\Gamma }(b-1)}\right](a+1\ne b,N>0)$$

there remains:

$$\begin{array}{l}F3{}_2\left(\begin{array}{c}1-N,2,m+1-\nu \\ -N,m+2\end{array};1\right)=\frac{\mathrm{\Gamma }(m+2)\mathrm{\Gamma }(-\nu )}{\mathrm{\Gamma }(m+1-\nu )\mathrm{\Gamma }(3-\nu )}\\ \left[(2N+2m-2\nu -N\nu +2\right]\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}(2m-N\nu )\frac{\mathrm{\Gamma }(N+m-\nu +2)}{\mathrm{\Gamma }(N+m+1)}\end{array}$$

For the second hypergeometric functions, it is yielded:

$$F3{}_2\left(\begin{array}{c}-1,2,-m\\ -N,\nu \end{array};1\right)=\frac{(N\nu -2m)}{N\nu }$$

After substitution in Equation (29), there remains:

$$\begin{array}{l}{I}_{\delta ,1}^{\nu }[f](x){\left[\frac{6}{N\mathrm{\Gamma }(N+3)\mathrm{\Gamma }(3-\nu )}\frac{1}{{\delta }^{\nu }}\right]}^{-1}=\\ ={\displaystyle \sum _{m=1}^{\infty }f(x-m\delta )}\left[(2N+2m-2\nu -N\nu +2)\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}-(2m+N\nu )\frac{\mathrm{\Gamma }(N+m-\nu +2)}{\mathrm{\Gamma }(N+m+1)}\right]+\\ +{\displaystyle \sum _{m=0}^{\infty }f(x-m\delta )}\frac{\mathrm{\Gamma }(N-m-\nu +2)}{\mathrm{\Gamma }(N-m+1)}(2m-N\nu )\end{array}$$

(30)

5. The Frequency Response for the Approximate Fractional Orthogonal Difference

Just as for the fractional orthogonal derivative, it is possible to compute the frequency response for the fractional orthogonal difference. In that case, we can apply Theorem 1, where we take another function for the function g. In this case, we use Equation (5):

$$g=I_{\delta }^{-\mu }[f](x)={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}}f(x-k\delta )\delta >0$$

Fourier transform of this function gives:

$$G(\omega )=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{i\omega x}}g(x)dx={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}}\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }{e}^{i\omega x}}f(x-k\delta )dx={\delta }^{\mu }{\displaystyle \sum _{k=0}^{\infty }\frac{{(\mu )}_k}{k!}}{e}^{i\omega k\delta }f(\omega )$$

The summation is known, and it follows:

$$G(\omega )={\left(\frac{1-e^{iw\delta }}{\delta }\right)}^{-\mu }F(\omega )$$

(31)

From Equations (4) and (31), we obtain that if $g=I_{\delta ,n,\omega }^v[f]$, then:

$$H(\omega )=\frac{G(\omega )}{F(\omega )}={\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{v-n}\frac{k_{n}n!}{h_n{\delta }^n}{\displaystyle \sum _{x=0}^N{p}_n(x)w(x)e^{-ix\delta \omega }}$$

For the Hahn polynomials, this becomes:

$$H(\omega )={\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{v-n}\frac{k_{n}n!}{h_n{\delta }^n}S$$

(32)

For the summation S, [2] (5.28) gives:

$$S=\frac{{(\beta +1)}_N}{N!}{(1-e^{-i\delta \omega })}^{n}F\left(\begin{array}{c}-N+n,\alpha +n+1\\ -\beta -N\end{array};e^{-i\delta \omega }\right)$$

Substitution in Equation (32) gives:

$$H(\omega )={\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }\frac{k_{n}n!}{h_n}\frac{{(\beta +1)}_N}{N!}F\left(\begin{array}{c}-N+n,\alpha +n+1\\ -\beta -N\end{array};e^{-i\delta \omega }\right)$$

Substitution of Equation (16) yields:

$$H(\omega )={(1)}^n{\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }\frac{\mathrm{\Gamma }(N+\beta +1)\mathrm{\Gamma }(2n+\alpha +\beta +2)}{\mathrm{\Gamma }(N+\beta +1)\mathrm{\Gamma }(N+n+\alpha +\beta +2)}F\left(\begin{array}{c}-N+n,\alpha +n+1\\ -\beta -N\end{array};e^{-i\delta \varpi }\right)$$

(33)

Comparing this formula with [2] (5.28), we will see that the only difference between these formulas is the factor ${(i\omega )}^{\nu -n}$. In the next section, we will discuss this fact.

With α = β = 0 and n = 1, there remains:

$$\begin{array}{l}H(\omega )=-{\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }\frac{6\mathrm{\Gamma }(N+1)}{\mathrm{\Gamma }(N+3)}F\left(\begin{array}{c}1-N,2\\ -N\end{array};e^{-i\delta \omega }\right)=\\ =-{\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }\frac{6\mathrm{\Gamma }(N)}{\mathrm{\Gamma }(N+3)}{\displaystyle \sum _{k=0}^{N-1}(k-N)(k+1)e^{-i\delta k\omega }}=\\ =\frac{{(1-e^{i\omega \delta })}^{\nu }}{{\delta }^{\nu }{(e^{-i\delta \omega }-1)}^3}\frac{6\mathrm{\Gamma }(N)}{\mathrm{\Gamma }(N+3)}(e^{-i\delta \omega (N+1)}(N-N{e}^{-i\delta \omega }+2)-(2e^{-i\delta \omega }+N{e}^{-i\delta \omega }-N))\end{array}$$

(34)

For N = 1, there remains:

$$H(\omega )=-{\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }$$

Remark 4. Note that the frequency response in Equation (33) is derived with a summation to infinity. For practical reasons, one should take the summation to a finite number, say M. Then, we have to use the frequency response corresponding to Equation (30). This is important, because Equation (30) is the formula to compute the approximate fractional Hahn derivative. Then, the frequency response will differ from the ideal one. Taking the Fourier transform, the result for the frequency response is:

$$\begin{array}{l}H(\omega ){\left[\frac{6}{\mathrm{\Gamma }(N+3)}\frac{1}{{\delta }^{\nu }}\right]}^{-1}=\\ =\frac{1}{\mathrm{\Gamma }(2-\nu )}{e}^{-i(N+1)\omega \delta }{\displaystyle \sum _{m=1}^{M+N+1}{e}^{im\omega \delta }}\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}+\frac{1}{\mathrm{\Gamma }(2-\nu )}{\displaystyle \sum _{m=1}^M{e}^{im\omega \delta }\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}}-\\ -\frac{2}{N}\frac{1}{\mathrm{\Gamma }(3-\nu )}{e}^{-iN\omega \delta }{\displaystyle \sum _{m=1}^{M+N}{e}^{im\omega \delta }}\frac{\mathrm{\Gamma }(m-\nu +2)}{\mathrm{\Gamma }(m)}+\frac{2}{N}\frac{1}{\mathrm{\Gamma }(3-\nu )}{\displaystyle \sum _{m=1}^M{e}^{im\omega \delta }\frac{\mathrm{\Gamma }(m-\nu +2)}{\mathrm{\Gamma }(m)}}\end{array}$$

(35)

For M →∞, this function becomes the same as Equation (34).

If f (x) is causal, then Equation (33) can be derived from a finite sum, and we do not cut offthe summation. Then, we can use the frequency response Equation (33) instead of Equation (35).

Remark 5. For the filter for the fractional differentiation of Grünwald and Letnikov, the approximate frequency response (with N finite) can be simply derived from Equation (11) by taking the Fourier transform. We first take the Fourier transform of Equation (11) without the limits. This gives: with N = [x/δ]

$$F[G{L}^{\nu }[f](x)(\omega )=\frac{1}{\mathrm{\Gamma }(-\nu )}{\displaystyle {\int }_0^{\infty }{e}^{i\omega x}}\frac{1}{\delta \nu }{\displaystyle \sum _{k=0}^{[x/\delta ]}\frac{\mathrm{\Gamma }(k-\nu )}{\mathrm{\Gamma }(K+1)}}f(x-k\delta )dx$$

Because f (x) is a causal function, it follows that x ≥ kδ. Substitution gives:

$$\begin{array}{l}H(\omega )=\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{1}{\delta \nu }{\displaystyle \sum _{k=0}^{\infty }\frac{\mathrm{\Gamma }(k-\nu )}{\mathrm{\Gamma }(k+1)}}{\displaystyle {\int }_{k\delta }^{\infty }{e}^{i\omega x}}f(x-k\delta )dx\\ =\frac{1}{\mathrm{\Gamma }(-\nu )}\frac{1}{\delta \nu }{\displaystyle \sum _{k=0}^{\infty }\frac{\mathrm{\Gamma }(k-\nu )}{\mathrm{\Gamma }(k+1)}}{e}^{i\omega k\delta }F(\omega )={\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }F(\omega )\end{array}$$

(36)

For δ → 0, the frequency response becomes:

$$H(\omega )=\underset{\delta \to 0}{\lim }{\left(\frac{1-e^{i\omega \delta }}{\delta }\right)}^{\nu }={(-i\omega )}^{\nu }$$

Because of Section 3.11 from [1], the Grünwald-Letnikov derivative follows by taking w (j) = 1 and N = n in Equation (12) or α = β = 0 and N = n in Equation (15). Therefore, Equation (36) becomes a special case of Equation (33).

6. Application of the Theory to a Fractional Differentiating Filter

This section treats the application of the fractional derivative in linear filter theory. In that theory, a filter is described by three properties, e.g., the input signal, the so-called impulse response function of the filter and the output signal. These properties are described in the time domain. The output signal is the convolution integral of the input signal and the impulse response function of the filter. The linearity means that the output signal depends linearly on the input signal.

One distinguishes between continuous and discrete filters, corresponding to continuity and discreteness, respectively, of the signals involved. In the discrete case, the output signal can be computed with a discrete filter equation.

In the usage of filters, there are two important items for consideration. The first item is the computation of the output signal given the input signal. The second item is a qualification of the working of the filter.

In our opinion, this second task should be preferably done in the frequency domain, where one can see the spectra of the signal and eventually the noise. In this domain, one obtains the Fourier transform of the output signal by multiplication of the Fourier transform of the input signal with the frequency response.

In this section, we will give graphs of the absolute value of the frequency response associated with various approximate fractional derivatives introduced in this paper. From these graphs, we will get insight into how these approximations will work in practice. We think this kind of analysis, in particular when using log-log plots, is preferable to the analysis of filters in the time domain. As an example of the latter type of analysis, actually for the fractional Jacobi derivative, see a paper of Tenreiro Machado [14].

For the second item, the frequency response of the filter should be computed. We treat here the approximate fractional Jacobi derivative and the approximate fractional Hahn derivative. For practical computations, these become the approximate fractional Legendre derivative and the approximate fractional Gram derivative.

For the first item, we can use [2] (4.18) for the Jacobi derivative and [2] (4.39) for the Hahn derivative. For the second item, we had to go to the frequency domain.

In the frequency domain, we suppose that the Fourier transforms of the input signal x (t) and the output signal y (t) are X (ω), and Y (ω) and does exist. The Fourier transformation of the impulse response function h (t) of the filter is the frequency response H (ω) of the filter. For this frequency response, there is the definition:

$$H(\omega ):=\frac{Y(\omega )}{X(\omega )}$$

The function H (ω) is a complex function. For a graphic display of this function, one uses the modulus and the phase of the frequency response. In our case, we want to see the property of the differentiation, and then, the graph of the modulus of the frequency response is preferred.

For an n-th order differentiator, one can show:

$$H(\omega )={(-i\omega )}^n$$

(37)

In the frequency domain, the graph of the modulus of the frequency response of this filter is a straight line with a slope that depends on n. See Figure 1.

We see that for ω→∞, the modulus of the frequency response |H (ω) | goes to ∞. This means that (because every system has some high-frequency noise) this filter is unstable. That is the reason why a differentiating filter should always have an integrating factor. This factor will be added to the filter, so that for ω→∞, the modulus of the frequency response goes to zero. The filter for the orthogonal derivatives has this property. These derivatives appear from an averaging process (least squares), and such a process will always give an integrating factor.

The formula for the frequency response of the approximate Legendre derivative is [1] (Section 5):

$$H(\omega )=\frac{\mathrm{\Gamma }(2n+2)}{2^n\mathrm{\Gamma }(n+1)}\frac{1}{{\delta }^n}{j}_n(\omega \delta )$$

where the functions j_n are the spherical Bessel functions [9] (10.49.3). The graph of the modulus of the first order approximate Legendre derivative with δ = 1 is given in Figure 2.

It is clear that for large ω, the modulus of the frequency response goes to zero. For small ω, the graph is a straight line with slope one. In the case that the order of the differentiating filter is not an integer, the formula of the frequency response should be analogous to Equation (37):

$$H(\omega )={(-i\omega )}^{\nu }$$

where υ can be an integer, a fractional or even a complex number. The graph of the modulus of the frequency response is given in Figure 3.

Furthermore, in this case, there is an instability. To prevent this instability, one can use an approximate fractional Jacobi differentiating filter. For the approximate Jacobi derivative, the frequency response is already computed in [2] (Section 5). We repeat the formula [2] (5.5).

$$H(\omega )={(-i\omega )}^v{e}^{-i\omega \delta }M(n+\alpha +1,2n+\alpha +\beta +2;2i\omega \delta )$$

With δ = 1, the squared absolute value of the frequency response of the approximate fractional Jacobi derivative can be written as a series expansion as follows:

$${\left|H(\omega )\right|}^2={\displaystyle \sum _{m=0}^{\infty }{(-4{\omega }^2)}^m}{\displaystyle \sum _{K=0}^{2m}\frac{{(n+\alpha +\beta +1)}_k{(n+\alpha +\beta +1)}_{2m-k}}{{(n+\alpha +\beta +1)}_k{(n+\alpha +\beta +1)}_{2m-k}}}\frac{{(-1)}^k}{k!(2m-k)!}$$

For small ω, a first approximation of the modulus of the frequency response gives:

$$|H(\omega )|~{\omega }^v\sqrt{1-\frac{4(n+\alpha +1)(n+\beta +1)}{{(2n+\alpha +\beta +2)}^2(2n+\alpha +\beta +3)}{\omega }^2}$$

(38)

It is clear that the formula is symmetric in α and β. The choices of α and β detect the cut-off frequency For simplicity, we look at the frequency (and not at the exact cut-off frequency) where the frequency response has a first maximum. From Formula (38), the maximum frequency is:

$${\omega }_{\max }=\frac{(2n+\alpha +\beta +2)}{2}\sqrt{\frac{\nu (2n+\alpha +\beta +3)}{(n+\alpha +1)(n+\beta +1)(1+\nu )}}$$

To simplify this formula, set α = β. Then, there remains:

$${\omega }_{\max }=\sqrt{\frac{\nu (2n+2\alpha +3)}{(1+\nu )}}$$

The shape of the curves does not change in principle. Therefore, if one wants the cut-off frequency as high as possible, then α and β should be chosen as high as possible. For the case α = β = 0, the formula for the frequency response of the approximate fractional Jacobi derivative simplifies to the approximate fractional Legendre derivative:

$$H(\omega )={\omega }^{\nu }\frac{\mathrm{\Gamma }(2n+2)}{2^n\mathrm{\Gamma }(n+1)}\frac{1}{{(\omega \delta )}^n}{j}_n(\omega \delta )$$

where the functions j_n are the spherical Bessel functions with n − 1 ≤ υ ≤ n. For n = 1 and δ =1, there remains:

$$H(\omega )=3{\omega }^{\nu -3}(\sin \omega -\omega \cos \omega )$$

The graph of the modulus of this frequency response is given in Figure 4 for different values of υ.

Next, we treat the discrete filter.

When we apply a discrete filter, the input signal and the output signal are sampled with a sample frequency that is higher than twice the maximum frequency of the input signal. Therefore, there is always a maximum frequency for the frequency response of the filter. What the frequency response will do above this maximum frequency is not important, presupposing that for ω →∞, the frequency response will go to zero. The input signal is known over N points. Therefore, the input signal is an approximation of the true input signal.

The output signal can be computed for the N points of the input signal. This is the reason that the output signal is always an approximation of the filtered true input signal. If the filter is a differentiator, the output signal is an approximation of the derivative of the input signal. This approximation is very dependent on N.

For the fractional orthogonal derivative, we use the frequency response of the approximate derivative of the discrete Hahn difference as derived in Equation (33). For α = β = 0, this frequency response becomes that of the fractional Gram derivative. In the following figure, we take with 0 < υ ≤ 1 and n = 1 in Formula (34). For N = 1, the modulus of the frequency response of the ideal filter has a maximum at ω = π. For different values of N, the modulus of the frequency response with δ = 1 is given in Figure 5.

In practice (taking a finite number of points), we use the frequency response in Equation (35). For δ =1, υ = 0.5, N = 7 and different values of M, the modulus of this frequency response is shown in Figure 6.

We see that for ω → 0, the absolute value of the frequency response tends to a constant value. Taking ω = 0 in Equation (35), it follows:

$$\begin{array}{l}H(0){\left[\frac{6}{\mathrm{\Gamma }(N+3)}\frac{1}{{\delta }^{\nu }}\right]}^{-1}=\frac{1}{\mathrm{\Gamma }(2-\nu )}{\displaystyle \sum _{m=1}^{M+N+1}\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}}+\frac{1}{\mathrm{\Gamma }(2-\nu )}{\displaystyle \sum _{m=1}^M\frac{\mathrm{\Gamma }(m-\nu +1)}{\mathrm{\Gamma }(m)}}\\ -\frac{2}{N}\frac{1}{\mathrm{\Gamma }(3-\nu )}{\displaystyle \sum _{m=1}^{M+N+1}\frac{\mathrm{\Gamma }(m-\nu +2)}{\mathrm{\Gamma }(m)}}+\frac{2}{N}\frac{1}{\mathrm{\Gamma }(3-\nu )}{\displaystyle \sum _{m=1}^M\frac{\mathrm{\Gamma }(m-\nu +2)}{\mathrm{\Gamma }(m)}}\end{array}$$

For the summations, it is proven:

$${\displaystyle \sum _{m=1}^K\frac{\mathrm{\Gamma }(m-\alpha )}{\mathrm{\Gamma }(m)}}=\frac{\mathrm{\Gamma }(K+1-\alpha )}{\mathrm{\Gamma }(K)(1-\alpha )}a\ne 1$$

Then, there remains:

$$\begin{array}{l}H(0){\left[\frac{6}{N\mathrm{\Gamma }(N+3)\mathrm{\Gamma }(4+\nu )}\frac{1}{{\delta }^{\nu }}\right]}^{-1}=\\ =(N-2M-N\nu )\frac{\mathrm{\Gamma }(M+N-\nu +3)}{\mathrm{\Gamma }(M+N+1)}+((3-\nu )N+2(M-\nu +2))\frac{\mathrm{\Gamma }(M-\nu +2)}{\mathrm{\Gamma }(M)}\end{array}$$

(39)

The lower bound of the frequency for which the filter does a fractional differentiation can be defined with the help of the following equation:

$${\omega }_l={\left(\frac{6(N-2M-N\nu )}{N\mathrm{\Gamma }(N+3)\mathrm{\Gamma }(4-\nu )}\frac{\mathrm{\Gamma }(M+N-\nu +3)}{\mathrm{\Gamma }(M+N+3)}+\frac{6((3-\nu )N+2(M-\nu +2))}{N\mathrm{\Gamma }(N+3)\mathrm{\Gamma }(4+\nu )}\frac{\mathrm{\Gamma }(M-\nu +2)}{\mathrm{\Gamma }(M)}\right)}^{1/\nu }$$

(40)

In practice, we should take this frequency 10-times higher than computed with Equation (40).

For the maximum frequency, we should compute the value of the frequency for the first maximum of the absolute value of Equation (35). A tedious computation leads to:

$${\omega }_{\max }\approx 2\sqrt{6}\frac{\sqrt{(1-\nu )(6N+\nu +6N\nu +N^2\nu +N^2+9)}}{\delta (6N+\nu +6N\nu +N^2\nu +N^2+9)}$$

Therefore, we can define a bandwidth B of the filter, which is equal to B = ω_max – ω_l.

For M →∞, we had to take the limits for the ratios of the gamma functions. The work in [9] (5.11.13) gives:

$$\begin{array}{l}\frac{1}{z^{a-b}}\frac{\mathrm{\Gamma }(z+a)}{\mathrm{\Gamma }(z+a)}=1+\frac{1}{2z}(a-b)(a+b-1)+\\ +\frac{1}{24z^2}(a-b-1)(a-b)(3{(a+b-1)}^2-(a-b+1))+O(z^{-3})\end{array}$$

Substitution in Equation (39) gives with δ =1:

$$H(0)=\frac{1}{M^{\nu }}\frac{(12N-10\nu -18N\nu +11{\nu }^2-3{\nu }^3+6N{\nu }^2-6N^2\nu +6N^2)}{2N(N+2)\mathrm{\Gamma }(N+1)(3-\nu )\mathrm{\Gamma }(1-\nu )}\left(1+O\left(\frac{1}{M}\right)\right)$$

Because ν > 0, this constant goes to zero, and the frequency response approximates the ideal case for low frequencies.

With this frequency response, it is seen that the filter works only well when the frequencies of the input signal lie inside the band-pass of the filter. The maximum frequency of this band-pass is mainly dependent on the number N and the minimum frequency on the number M.

Many authors (e.g., [15–18]) tried to describe the properties of a discrete fractional filter in the time domain with causal functions. Then, there are always transient effects at the initial time t = 0. There, the filter has other properties than elsewhere, because of the discontinuity of the input signal (and the resulting high frequencies).

When using the frequency response, the realms of lower and of higher frequencies need special attention. For the lower frequencies, the frequency response does not go to zero when using a finite number of points (in practice, this is always the case). This we call the minimum frequency effect. When the input signal is causal, the frequency response has no minimum frequency effect (this is important when, for example, treating differential equations). For the higher frequencies, there is always a maximum depending on the sample frequency (Shannon frequency). However, for the orthogonal derivative, there is an extra maximum that is lower than the Shannon frequency. For the GL filter, there is no such maximum.

7. Conclusions

The first conclusion of this paper is that when using a discrete fractional filter, the properties of the filter should be examined using the frequency response. The second conclusion is that a discrete fractional differentiating orthogonal filter has always a certain bandwidth.

In general, a fractional differentiating filter can be built from two serial filters. The first one has a frequency response with modulus |ω|^υ. The second filter is a so-called low-pass filter. This low-pass filter determines the highest frequency for differentiation of the filter. In the case of the Jacobi filter (continuous case), one has to keep in mind that the frequency response has side lobes. In the discrete case, these side lobes are not important, because there is a maximum sample frequency. The Jacobi filter can be used if the maximum signal frequency is far beyond the side lobes. These sides lobes go to zero for high frequencies. Hence, this is a very good fractional differentiating filter.

As another example for the low-pass filter, we can choose the Butterworth filter. For this filter, the frequency response is:

$$H_n(\omega )={(-i\omega )}^{\nu }\left(\frac{1}{1+{(\omega /{\omega }_0)}^{2n}}\right)$$

which gives a very good frequency response for fractional differentiation. There are no side lobes. In Figure 7, this is demonstrated. There is no simple analogous discrete filter.

The work in [1] makes some remarks about the practical applications of the filters obtained from orthogonal polynomials and the Butterworth filter. In the analog integer case, the Butterworth filter can be much more easily constructed physically In the analog fractional case, a filter should be built with the frequency response (−iω)^υ. This gives a difficult problem. Solutions are always approximations.

In the discrete case, the Butterworth frequency response can be transformed into a filter equation. See [19]. For more details about the Butterworth filter, see [20] (Section 5.2.1) (both analog and discrete), [21] (both analog and discrete), [22] (Section 3.2) (analog) and [19] (Section 12.6) (digital).