AlgDiff: an open source toolbox for the design, analysis and discretisation of algebraic differentiators

3.1 Derivative estimation

While the early works [15, 16] use differential algebraic manipulations and operational calculus to derive the differentiators, later ones such as [20], [21], [22, 25, 28] use time-domain derivations. This shall be briefly recalled.

Let I t = [ t − T , t ] ⊂ R , T > 0, be an arbitrary closed interval and consider a Lebesgue integrable signal y : I t → R . The derivatives up to a finite order n < min{α, β} + 1, with α , β ∈ R and α, β > − 1, of y can be approximated using truncated generalised Fourier expansions of order N based on orthogonal Jacobi polynomials as

with g N , T , ϑ ( α , β ) defined in (A.7). In the following, the kernel g N , T , ϑ ( α , β ) and the parameter T are called algebraic differentiator and window length, respectively. The scalars α and β parametrise the weight function of the Jacobi polynomials. Section 3.2 discusses their effects on the filter properties.

The estimate y ̂ ( n ) corresponds to a delayed approximation of y ⁽ⁿ⁾ at time t and the delay, denoted by δ _t , can be parametrised by ϑ. It can be expressed as

If ϑ corresponds to a zero of the polynomial P N + 1 ( α , β ) , which yields a small but known delay, the order of the approximation is increased, as first pointed out in [16] (see also [22, 26]). Alternatively, a delay-free approximation can be achieved for ϑ = 1 or a prediction for δ _t < 0 (see also [22, 26, 29]).

3.2 Interpretations

As discussed in [26], different interpretations can be attached to algebraic differentiators. Two important ones are recalled in Table 1. The approximation-theoretic interpretation stems from the time-domain derivation of the differentiators. It is helpful for the analysis of the approximation error and delays. The systems-theoretic interpretation is a direct consequence of the definition of the approximated derivative in (1): The estimate is the output of a linear time invariant filter where the kernel is the nth order derivative of g N , T , ϑ ( α , β ) . Applying an integration by parts in (1) shows that the estimate can be interpreted as the output of a filter with kernel g N , T , ϑ ( α , β ) . The filter output, then, is the sought derivative. This interpretation will be used in the next section to design differentiators using the toolbox.

Table 1:

Extract from [26]: interpretations of g = g N , T , ϑ ( α , β ) and the practical usage of each.

Context	Interpretation	Practical usage
Approximation-theoretic	Polynomial approximation of derivative using a Hilbert reproducing kernel	Analysis of estimation properties and relation to established differentiation methods
Systems-theoretic	Linear time-invariant filtering of the sought derivative:	Analysis of filter properties and parametrisation
	Linear time-invariant filtering of the measured signal:	Discretisation and implementation

4.1 Impulse and step responses

The impulse and step responses of an algebraic differentiator g N , T , ϑ ( α , β ) have been derived in [22] and are summarised in Appendix C. Example 1 shows how the toolbox can be used to easily compute impulse responses, their derivatives and step responses of an algebraic differentiator.

Example 1

(Impulse and step responses). The impulse and step responses of an algebraic differentiator with parameter α = 4, β = 4, N = 0, and T = 0.1 can be computed using the toolbox as shown in Code snippet 1.

Code snippet 1:

Evaluation of the impulse and step responses of an algebraic differentiator.

Figure 2 shows the impulse and step responses of algebraic differentiators for N = 0 and different values of α and β. It can be observed that increasing β and keeping α constant causes a stronger weighting of the impulse response towards 0. This increase implies a stronger weighting of more recent values of the sought derivative and corresponds to the reduction of the estimation delay. In contrast, increasing α and keeping β constant yields a stronger weighting of the impulse response in the direction of T, which leads to a stronger weighting of older values of the sought derivative and, thus, to an increase of the estimation delay. The effects of the parameters on the delay has been analysed analytically in [16, 20, 25, 26], for example.

Figure 2:

Impulse and step responses of algebraic differentiators with parameters α, β ∈ {1, 3, 6}, and N = 0.

The analytical properties of the impulse response and its derivatives can be used to derive tuning guidelines for fault detection algorithms, for example. Approaches have been developed and experimentally validated in [21, 22, 28, 30], for example.

4.2 Frequency-domain properties

By recalling that the estimate of the nth order derivative of y is a convolution product of the signal and the nth order derivative of the differentiator g N , T , ϑ ( α , β ) , it follows that

with F y ̂ ( n ) , G N , T , ϑ ( α , β ) , and Y the Fourier transforms of y ̂ ( n ) , g N , T , ϑ ( α , β ) , and y, respectively. For the sake of completeness, the closed form of G N , T , ϑ ( α , β ) is recalled in (A.9). Thus, this numerical differentiation scheme falls within the framework of classical approaches for differentiation in the frequency domain, where a smoothing filter, in this case g N , T , ϑ ( α , β ) , is followed by an ideal differentiation operator (ιω) ⁿ . This property has first been observed in [23].

Analysing G N , T , ϑ ( α , β ) is thus helpful to derive tuning guidelines by specifying desired frequency-domain properties. The approximation

of the amplitude spectrum of an algebraic differentiator has been proposed in [21], where ω _c is defined in (A.11) and μ = 1 + min{α, β}. It follows that G N , T , ϑ ( α , β ) can be interpreted as a low pass filter: The frequency ω _c and the set ω < ω c are called cutoff frequency and passband, respectively. The set ω > ω c is called stopband.

Tuning guidelines proposed in [21, 22] can then be inferred from (2): A desired cutoff frequency ω _c and a desired filter order, i.e., the stopband slope given by 20μ dB, can be specified. Then, the filter window length can be computed from (A.11).

It has been shown in [24–26] that increasing N increases the sensitivity to noise. It has also been shown in [24] that the choice α = β yields differentiators with maximum robustness to noise. In addition, according to [22, 25], and [26], choosing α ≠ β reduces stopband ripples in the amplitude spectrum.

Example 2 shows how a differentiator can be designed from desired frequency-domain properties using the toolbox. The computation of the amplitude and phase spectra is also demonstrated. Example 3 compares and discusses the amplitude spectrum of differentiators with and without delay.

Example 2

(Specifying frequency-domain properties). An algebraic differentiator with an amplitude spectrum asymptotically decreasing by 40 dB per decade, i.e., min{α, β} = 1, and a cutoff frequency ω _c = 100 rad s⁻¹ can be designed as described in Code snippet 2. The example also shows how the toolbox can be used to compute the amplitude spectrum, its approximation (2), and to return the resulting window length T. To increase the robustness with respect to corrupting noises the choices N = 0 and α = β are made. Additionally, upper and lower bounds for G N , T , ϑ ( α , β ) can be computed. They have been discussed in [22, 25, 26].

Code snippet 2:

Design of an algebraic differentiator from desired frequency-domain properties.

The amplitude spectrum, the approximation, the bounds, and the phase spectrum are given in Figure 3. Note that for this parametrisation the amplitude spectrum has infinitely many zeros. Thus, its lower bound is equal to zero (see also [22]).

Figure 3:

Amplitude and phase spectra of the differentiator from Code snippet 2. The approximation (2) and the upper bound from [22, 25] are also given.

Example 3

(Differentiators with and without delay). As discussed in Section 3, algebraic differentiators can also be parametrised such that the estimated derivative is delay-free. This is done in Code snippet 3 for a differentiator with N = 1 and α = β = 1. Its amplitude spectrum is also computed. For a comparison a differentiator with delay but the same parameters α, β, and N is also considered.

Code snippet 3:

Comparison of amplitude spectra for differentiators with and without delay.

Figure 4 shows the amplitude spectra of both differentiators from the latter example. The delay-free differentiator shows a significant overshoot. Thus, the resulting estimated derivative might not be usable. The decrease in the estimation quality of these differentiators has been remarked in [31]. The overshoots in the amplitude spectra of delay-free differentiators have first been observed in [23].

Figure 4:

Amplitude and phase spectra of the differentiators with and without delay from Code snippet 3.

It has been observed in [21, 22] that the Fourier transform G N , T , ϑ ( α , β ) of an algebraic differentiator with N = 0 and α = β simplifies to

where Γ denotes the Gamma Function and J α + 1 2 is the Bessel function of the first kind and order α + 1/2. These special functions are defined in [32], for example. Thus, the amplitude spectrum exhibits an infinite number of positive zeros. Moreover, the differentiator provides phase linearity. By carefully choosing the parameters α and T such that the zeros of G 0 , T , ϑ ( α , α ) coincide with a specific frequency, the exact annihilation of disturbing harmonics can be performed. This parametrisation has been used in experimental studies in [21, 22, 29, 30, 33]. Therein, the measured signals are corrupted by an unmodelled mechanical eigenmode of the system. Its influence on the control and estimation algorithms has been annihilated using the latter property of algebraic differentiators. In [22], this parametrisation has been used to approximate a matched filter. The following example demonstrates this special parametrisation.

Example 4

(Exact annihilation of frequencies). Assume that the signal the derivative of which has to be estimated is corrupted by a harmonic signal of known frequency ω ₀. Choosing the window length T such that ω ₀ T = 2j _k , with j _k the kth zero of the Bessel function of the first kind and order α + 1/2 yields G 0 , T , ϑ ( α , α ) ( ω 0 ) = 0 , i.e., the effects of the corrupting signal are exactly annihilated. This is demonstrated in Code snippet 4 for ω ₀ = 100 rad s⁻¹ and α = 2.

Code snippet 4:

Parametrisation of a differentiator to annihilate a known frequency.

Figure 5 shows the amplitude and phase spectra of this differentiator and confirms the exact annihilation of the disturbing harmonic with frequency ω ₀. However, choosing an appropriate zero and a numerical value for α remain as the degrees of freedom for every specific application. The works [21, 22, 29, 30, 33] have proposed an approach where these parameters are chosen such that a satisfying compromise between disturbance rejection and tolerable window length can be achieved. The size T of the window lengths is a crucial parameter since it influences the estimation delay, the computational burden, and the memory requirements of the implemented differentiators. Available hardware may imply constraints on the latter properties.

Figure 5:

Amplitude and phase spectrum of the differentiator from Code snippet 4 for the annihilation of a harmonic disturbance with angular frequency ω ₀.

4.3 Discrete-time implementation

In most applications, the measured signal is available at discrete time instants only and the integral (1) must be approximated by an appropriate quadrature method. Thus, the estimated derivatives are outputs of digital FIR filters.

From the Nyquist–Shannon sampling theorem, it follows that when discretising a continuous-time filter, the gain at the Nyquist frequency must be sufficiently small to keep aliasing effects low. As proposed in [22], a possibility to avoid these effects is to specify the weakest relative attenuation

of an algebraic differentiator g N , T , ϑ ( α , β ) at the Nyquist frequency ω _N = π/t _s, where t _s and n _max are the sampling period and the highest required derivative order, respectively (see also [26]). Example 5 demonstrates the design of a differentiator with a desired relative attenuation k _N,min.

Example 5

(Specifying the relative attenuation k _N,min). To achieve at least a desired relative attenuation of k _N,min for a given cutoff frequency ω _c, the parameters α and β have to be chosen such that

min { α , β } = ln k N , min ln ω c / ω N + n max − 1 .

A differentiator is designed in Code snippet 5 for the approximation of a first order derivative with the cutoff frequency equal to 90 rad s⁻¹ and a relative attenuation equal to 10⁻³ = −30 dB. The considered sampling period is equal to 10 ms, the parameters α and β are chosen equal, and N = 0.

Code snippet 5:

Design of a differentiator with a specified relative attenuation.

The resulting differentiator has a filter window length equal to 100 ms, i.e., 10 sampling periods. The parameters α and β are equal to 5.53. The estimation delay is equal to 5 sampling periods.

In the following, it is assumed that the window length T is an integral multiple of the sampling period, i.e., T = n _s t _s, and equidistant sampling is considered for simplicity. The abbreviation y _k = y(kt _s), k ∈ N , for a sample of a signal y at time kt _s is used. Then, a discrete-time approximation

of (1) can be achieved. Therein, θ, L, and w _i depend on the numerical integration method used as discussed in [22, 26, 28]. For example, the mid-point rule yields θ = 1 2 , w i = t s g i + θ ( n ) , and L = n _s. The normalisation factor Φ has been first introduced in [22] to guarantee that the DC component of the sought derivative is preserved. The effects of normalisation are demonstrated in the following example using the functions of the toolbox.

Example 6

(Effect of the normalisation factor). Per default, the toolbox implements the normalisation presented in (3). The code in Code snippet 6 demonstrates the effects of the normalisation by comparing the amplitude spectra of two differentiators: The first one is discretised with the normalisation and the second is without.

Code snippet 6:

Discretisation of differentiators with and without normalisation.

Figure 6 shows the spectra of two continuous-time differentiators discretised with and without the normalisation factor. It can be seen that the DC component is only preserved when the normalisation is used.

Figure 6:

Amplitude of two continuous-time differentiators compared to those of discretised ones with and without normalisation.

Algebraic differentiators and their tuning guidelines have been discussed in the continuous-time domain. Thus, the discretised ones have to preserve the specified properties. Since most tuning guidelines discussed here rely on the frequency-domain interpretation of the differentiators, the analyses of the discretisation effects can be performed by comparing the Fourier transform of the continuous and discrete differentiators. This can be done using the cost function

introduced in [22] (see also [26, 28]), where Ω has to be chosen according to the frequency interval of interest. In the latter cost function, ω ↦ F g ( n ) ( ω ) is the Fourier transform of the continuous-time differentiator g = g N , T , ϑ ( α , β ) . The discrete-time filter from (3) is denoted g ̂ ( n ) and its Fourier transform is

The cost function J compares the discretisation error with the noise amplification of the continuous-time differentiator (the reader is referred to [22, 26] for further details). Thus, it can be used to decide if the former can be neglected compared to the latter. The cost function can be easily computed using the functions in the toolbox as demonstrated in Example 7.

Example 7

(Discretisation error). The cost function (4) is computed in Code snippet 7 for a differentiator with the parameters α = β = 7, N = 0, and T = 20t _s.

Code snippet 7:

Evaluation of the discretisation error.

The resulting cost J is approximately equal to 1.6 ⋅ 10⁻¹⁰. Thus, the discretisation effects can be neglected. The variation of J with respect to the parameters of the differentiator has been analysed in [22, 28], for example.

4.4 Estimation of derivatives using AlgDiff

Example 8 demonstrates the numerical differentiation of a given signal using the implemented functions and a differentiator with desired frequency-domain properties.

Example 8

(Numerical differentiation). An algebraic differentiator with the cutoff frequency ω _c = 20 rad s⁻¹ and a desired relative attenuation k _N, min = 10⁻³ shall be designed for the estimation of the first order derivative of a signal sampled with a sampling period equal to 20 ms. As in Example 2, the choice N = 0 and α = β is made to increase the robustness with respect to disturbances. The resulting differentiator has the parameters α = β = 3.35 and a window length equal to 14t _s. The estimation delay corresponds to 7t _s. The cost function J is equal to 7.9 ⋅ 10⁻⁶. Thus, discretisation effects can be neglected. The parametrisation of the differentiator and the estimation of the derivative of a noisy signal t ↦ y(t) = sin(t) + η(t) is shown in Code snippet 8. The additive disturbance η is drawn from a Gaussian distribution with standard deviation equal to 0.02 such that the signal to noise ratio is SNR = log 10 ∑ i y 2 ( t i ) ∑ i η 2 ( t i ) ≈ 31 dB .

Code snippet 8:

Approximation of a derivative of a given signal using a differentiator designed with desired frequency-domain properties.

The estimation results are shown in Figure 7 and compared to the results of a simple forward difference method. For a comparison with more advanced methods the reader is referred to [25, 34, 35]. This example shows the advantages of algebraic differentiators and their systematic tuning guidelines. Furthermore, it highlights the user friendliness of the toolbox AlgDiff.

Figure 7:

Estimation of the first order derivative x ̇ of a sinusoidal signal x using a measurement y using the forward-difference method (FD) and the algebraic differentiator (Alg.) designed in Code snippet 8.