Algebraic differentiators are linear time-invariant filters with a finite-duration impulse response. These filters can be approximated as lowpass filters with a known cutoff frequency and a known stopband slope. They depend on 5 parameters:
The effects of the parameters can be very briefly summarized as follows (see the second notebook for a detailed discussion):
In the following example, the first derivative of a signal is estimated.
%matplotlib widget
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
from AlgDiff import *
import numpy as np
import sympy as sp
################################################################
# Define signal x and its measurement y=x+eta with eta
# a random disturbance
################################################################
# Define sampling rate
ts = 0.01
# Define signal and its derivatives
a = sp.symbols('a_0:3')
t = sp.symbols('t')
x = a[0]*sp.exp(-a[1]*t)*sp.sin(a[2]*t)
# Derivative to be estimated
dx = sp.diff(a[0]*sp.exp(-a[1]*t)*sp.sin(a[2]*t),t,1)
d2x = sp.diff(a[0]*sp.exp(-a[1]*t)*sp.sin(a[2]*t),t,2)
aeval = {'a_0':1,'a_1':0.1,'a_2':4}
# Evaluate signal and true derivative
teval = np.arange(0,20,ts)
for ai in a:
x = x.subs({ai:aeval[repr(ai)]})
dx = dx.subs({ai:aeval[repr(ai)]})
d2x = d2x.subs({ai:aeval[repr(ai)]})
xeval = sp.lambdify(t, x, "numpy")
xeval = xeval(teval)
dxeval = sp.lambdify(t, dx, "numpy")
dxeval = dxeval(teval)
d2xeval = sp.lambdify(t, d2x, "numpy")
d2xeval = d2xeval(teval)
# Get measurement
eta = np.random.normal(0,0.01,len(xeval))
y = xeval + eta
################################################################
# Initialize a two differentiators with a desired cutoff
# frequency with different alpha and beta. See later the effect
# of these parameters after the discretization
################################################################
wc = 30
algDiff1 = AlgebraicDifferentiator(N=1,alpha=3,beta=3,T=None,wc = wc, ts=ts)
algDiff2 = AlgebraicDifferentiator(N=1,alpha=1.1,beta=1.1,T=None,wc = wc, ts=ts)
# Discretization method
method = "trapezoidal"
# For a delay-free approximation uncomment the followning lines
# Note: A delay-free approximation is less accurate and less
# robust with respect to disturbances
# algDiff1.set_theta(1,False)
# algDiff2.set_theta(1,False)
################################################################
# Frequency-domain analysis
################################################################
omega = np.linspace(1,np.pi/ts,10**3)
# Get phase and amplitude of Fourier transform of the continous
# and discrete time filter for the estimation of the first
# derivative.
# Note: get_ampAndPhaseFilter return the Fourier transform of
# the continuous-time kernel. To compute the filter for the
# n-th derivative multiply with (j\omega)^n, with j the imaginary unit.
# The function get_ampSpectrumDiscreteFilter returns the amplitude
# of the discrete-time filter for the specified derivative order.
ampCont1,phaseCont1 = algDiff1.get_ampAndPhaseFilter(omega)
ampCont1 = ampCont1*(1j*omega)
phaseCont1 = phaseCont1 + np.angle((1j*omega))
ampDiscrete1,phaseDis1 = algDiff1.get_ampSpectrumDiscreteFilter(omega,1,method=method)
ampCont2,phaseCont2 = algDiff2.get_ampAndPhaseFilter(omega)
ampCont2 = ampCont2*(1j*omega)
phaseCont2 = phaseCont2 + np.angle((1j*omega))
ampDiscrete2,phaseDis2 = algDiff2.get_ampSpectrumDiscreteFilter(omega,1,method=method)
# Plot results
fig, ax = plt.subplots(nrows=1, ncols=1,sharex=False, figsize=(10, 5))
l = 3
fig.suptitle("Frequency analysis of the continuous and the discrete time filters")
ax.plot(omega,abs(ampCont1),label='continous filter 1',linewidth=l)
ax.plot(omega,ampDiscrete1,'-.',label='discrete filter 1',linewidth=l)
ax.plot(omega,abs(ampCont2),label='continous filter 2',linewidth=l)
ax.plot(omega,ampDiscrete2,'-.',label='discrete filter 2',linewidth=l)
ax.set_xlabel(r"$\omega$ in rad/s")
ax.set_ylabel(r"amplitude")
ax.legend()
ax.grid()
plt.text(300, 50, "Discretization changes\n filter characteristic in\n the frequency-domain!", size=20, va="top",
bbox=dict(boxstyle="round",
ec=(0.5, 0.5, 0.5),
fc=(0.9, 0.9, 0.9),
)
)
plt.show()
fig, ax = plt.subplots(nrows=1, ncols=1,sharex=False, figsize=(10, 5))
l = 3
fig.suptitle("Frequency analysis of the continuous and the discrete time filters")
ax.plot(omega,phaseCont1,label='continous filter 1',linewidth=l)
ax.plot(omega,phaseDis1,'-.',label='discrete filter 1',linewidth=l)
ax.plot(omega,phaseCont2,label='continous filter 2',linewidth=l)
ax.plot(omega,phaseDis2,'-.',label='discrete filter 2',linewidth=l)
ax.set_xlabel(r"$\omega$ in rad/s")
ax.set_ylabel(r"phase spectra in rad ")
ax.legend()
ax.grid()
plt.show()
################################################################
# Estimate derivatives
################################################################
# Estimate 0-th derivative
xApp = algDiff1.estimateDer(0,y,method=method)
# Estimate first derivative
dxApp = algDiff1.estimateDer(1,y,method=method)
The differentiator has the parameters: Alpha: 3.000000 Beta: 3.000000 Window length in s: 0.250000 Sampling period in s: 0.010000 Polynomial degree: 1 Estimation delay in s: 0.083333 Cutoff Frequency in rad/s: 30.454033 Cutoff Frequency in Hz: 4.846910 Discrete window length: 25 The differentiator has the parameters: Alpha: 1.100000 Beta: 1.100000 Window length in s: 0.150000 Sampling period in s: 0.010000 Polynomial degree: 1 Estimation delay in s: 0.042110 Cutoff Frequency in rad/s: 30.467980 Cutoff Frequency in Hz: 4.849130 Discrete window length: 15
################################################################
# Plot results
################################################################
fig, (fy,fdy) = plt.subplots(nrows=1, ncols=2,sharex=True, figsize=(10, 5))
fig.subplots_adjust( wspace=0.5)
fy.plot(teval,xeval,label='x')
fy.plot(teval,y,label='y')
fy.plot(teval,xApp,label='output filter')
fy.set_xlabel(r'$t$')
fy.set_ylabel(r'signals')
fdy.plot(teval,dxApp,label='output filter')
fdy.plot(teval,dxeval,label='true signal')
fdy.set_xlabel(r'$t$')
fdy.set_ylabel(r'first der. of signals')
plt.legend()
plt.show()
In most embedded applications of algebraic differentiators, the measured signal $y$ is available at discrete sampling instants only. Then, the finite convolution integral of the estimation (see the extended examples) must be approximated by an appropriate quadrature method. This yields discrete FIR filters, whose implementation is favourable in terms of numerical stability.
Let $t_{\mathrm{s}}$ be the sampling period. For the sake of brevity, the abbreviation $f_i = f(i t_{\mathrm{s}})$, $i\in\mathbb{N}$, for a sample of a function $f$ at time $i t_{\mathrm{s}}$ is used in the following. Then, various discrete-time approximations yield at a step $k$ an estimate of the form \begin{equation} \hat{y}^{(n)}_{k}=\sum_{i=0}^{L-1}w_iy_{k-i},\quad\quad k\geq L, \end{equation} where $L$, the number of filter coefficients, and $w_i$, $i=0,...,L$, the filter coefficients, depend on the used numerical integration method. For more details see the documentation of the class and the corresponding literature.
The following example shows how to get the filter coefficients $w_i$ using a given discretization method.
################################################################
# Saving filer coefficients for the estimation of the first
# derivative using the discretization method specified above
################################################################
coeff = algDiff1.discretize(1,method)
w = coeff[1][method]
print(w)
[ 0. 0.33512151 1.02262171 1.70955167 2.1827285 2.33934223 2.16033687 1.68656557 0.9977201 0.1940344 -0.61923765 -1.34457025 -1.90213953 -2.23702326 -2.32362626 -2.16733177 -1.80337856 -1.29296381 -0.7165719 -0.16452897 0.27521673 0.53208862 0.57064979 0.41114518 0.15281794 0. ]