Discretization Techniques

In Digital Filter-Implementation we have shown how digital filters are actually implemented using difference equations. In the blogs, Filters and Controllers in control system, we have shown a couple of basic filters and the PID controller in the continuous time domain. Finally, in Digital Control Systems 1.01 and 1.02 we shown what comes along with digital systems.

In this article we show how to transfer, transfer functions defined in the continuous-time domain to discrete-time domain; discretization. There are several techniques to accomplish this; i) modeling the sampled system or ii) by using Z-transform mappings.

  1. Modeling the sampled system
    • Zero-order hold
    • First-order hold
    • Other higher order hold circuits
  2. Z-transform mappings
    • Forward euler
    • Backward euler
    • Trapeziod / Tustin
    • Tustin with pre-warp
    • Matched pole zero
    • Modified matched pole zero

In this article we only go in detail in the Z-transform mapping techniques.

Finite difference methods

The first three methods actually use the finite difference methods to approximation the derivative and obtain a mapping from the continuous (s = j\omega_c) to the discrete (z = e^{j\omega_dh}) time domain. Herein, the forward and backward Euler methods use the following approximations

    \[\dot{x}(kh) \approx \frac{x(kh +h) + x(kh)}{h} \rightarrow s = \frac{z - 1}{h}\]

    \[\dot{x}(kh) \approx \frac{x(kh) - x(kh - h)}{h} \rightarrow s = \frac{1 - z^{-1}}{h}\]

Whereas the trapeziod method, also referred as the Tustin’s method, uses the approximation

    \[\frac{\dot{x}(kh) + \dot{x}(kh + h)}{2} \approx \frac{x(kh + 1) - x(kh)}{h} \rightarrow s = \frac{2}{h}\frac{z-1}{z+1}\]

Note that the forward Euler, stable poles in s may be mapped to unstable poles in z. For both backward Euler as Tustin’s method, stable poles in s implicates stable poles in z.

Tustin with pre-warping

When determining discrete-time frequency response of the Tustin rule, H_d(z) with z = e^{j\omega_dh}, to a continuous time transfer function H_c(s), it reveals a relation between the continuous and discrete time domain; Every frequency point in continuous time domain is mapped to the discrete time domain via \omega _{c}={\frac{2}{h}}\text{tan}\left(\omega_dh/2\right) \Leftrightarrow \omega_d = \frac{2}{h}\text{arctan}\left(\omega_ch/2\right).

    \begin{gather*}H_d(z) = H_c(s)|_{s = \frac{h}{2}\frac{z-1}{z+1}} \Leftrightarrow H_d(z) = H_c\left(\frac{h}{2}\frac{z-1}{z+1}\right) \\ \Leftrightarrow H_d(e^{j\omega_dh}) = H_c\left(\frac{h}{2}\frac{e^{j\omega_dh}-1}{e^{j\omega_dh}+1}\right) \Leftrightarrow H_d(e^{j\omega_dh}) = H_c\left(j\frac{2}{h}\text{tan}\left(\frac{\omega_dh}{2}\right)\right)\end{gather*}

The tangent introduces a non-linear mapping where the distortion is small when \omega_d h is small. We can compensate for this by setting \omega_c =\frac{2}{h}\text{tan}\left(\omega_dh / 2\right) for every frequency we have control over. That is to say, first replace the critical frequencies in the continuous-time transfer functions via

    \begin{gather*}(s + \omega) \rightarrow (s + \hat{\omega}) \text{ with } \hat{\omega} = \frac{2}{h}\text{tan}\left(\omega h / 2\right) \\(s^2 + 2\beta\omega s + \omega^2) \rightarrow (s^2 + 2\beta\hat{\omega} s + \hat{\omega}^2) \text{ with } \hat{\omega} = \frac{2}{h}\text{tan}\left(\omega h / 2\right)\end{gather*}

Secondly apply the normal Tustin transformation s = \frac{2}{h}\frac{z-1}{z+1}. Finally, scale the static gain of the discrete transfer function to match the gain of the continuous-time transfer function, e.g., H_c(j\omega_0) = k \cdot H_d\left(e^{j\omega_0 h}\right) with \omega_0 = 0 and solving for k.

Remark that most programs, such as Matlab only offer to replace one critical frequency via their discretization (c2d) function. They apply the transformation

    \[s = \frac{\omega_0}{\text{tan}\left(\omega_0 h / 2\right)}⁡ \frac{z - 1}{z + 1}\]

In which \omega_0 is the frequency chosen by the user to be matched in both continuous as discrete time domain.

Matched Pole Zero

This method maps the poles and zeroes according the the relation z = e^{sh}. That is to say,

    \begin{gather*}(s + \omega) \rightarrow z - e^{-\omega h} \\(s + a)^2 + b^2 \rightarrow z^2 - 2\left(e^{-ah}\text{cos}\left( b h \right)\right)z + e^{-2ah}\end{gather*}

If the numerator (zeros) is of lower order than the denominator (poles), add powers of (z + 1) to the numerator, until numerator and denominator are of equal order. This causes an averaging of the current and past input values, as in Tustin’s method.

However, it might be possible that it is desirable to have the output only depend on the past values. For instance, due to costly computing power. In that case, only add as much powers of (z + 1) such that the numerator is of lower order than the denominator by 1. In which case, this method is referred to as the modified matched pole zero.

Likewise, as with the Tustin transform with pre-warping, scale the static gain of the discrete transfer function. Note that this method preserves stability and provides less high frequency response error.

Frequency response example

The figure below shows the frequency response of a continuous time transfer function as well as the discretized transfer functions using the various methods discussed. The continuous-time transfer function PID controller with a second order low-pass filter and two notches located at 950 and 1150 Hz. A sampling frequency of 5000 Hz was used.

The forward and backward Euler methods cannot track the controller correctly at around 300 Hz. Furthermore, they are unable to replicate the fast dynamical behavior of the two notches. As such, this is a clear example of why one should not use this methods. Tustin’s method shows the same dynamics as the continuous time transfer function, but due to the distortion in the frequency domain the notches are not located at the correct frequency. Tustin’s method combined with pre-warping is able to track the notches but loses magnitude when reaching the Nyquist frequency. In addition there is a phase offset. The Matched Pole Zero method is able to match the magnitude response very well, but has a large difference in phase which depending on the performance criteria might become problematic.

The Matlab code below was used to generate the plots.

clear all;
close all;


% Init variables and define continuous-time system
s = sym('s');
z = sym('z');
fs = 5e3;
h  = 1/fs;

Kp    = 1;
fd    = 100;
fi    = 10;
flp   = 400;
wlp   = 2*pi*flp;
blp   = 1;
fp{1} = 950;
fz{1} = 950;
bp{1} = 0.1;
bz{1} = 0.01;
wp{1} = 2*pi*fp{1};
wz{1} = 2*pi*fz{1};
fp{2} = 1150;
fz{2} = 1150;
bp{2} = 0.1;
bz{2} = 0.005;
wp{2} = 2*pi*fp{2};
wz{2} = 2*pi*fz{2};
Gc.sym = (wp{1}^2)/(wz{1}^2)*(s^2 + 2*bz{1}*wz{1}*s + wz{1}^2) / ...
                             (s^2 + 2*bp{1}*wp{1}*s + wp{1}^2) * ...
         (wp{2}^2)/(wz{2}^2)*(s^2 + 2*bz{2}*wz{2}*s + wz{2}^2) / ...
                             (s^2 + 2*bp{2}*wp{2}*s + wp{2}^2) * ...                     
         (Kp*(1/(2*pi*fd)*s + 1 + (2*pi*fi)/s)*(2*pi*flp)^2 / ...
         (s^2 + 2*blp*(2*pi*flp)*s + (2*pi*flp)^2));

[Gc.num,Gc.den] = numden(Gc.sym);
Gc.tf = tf(sym2poly(Gc.num),sym2poly(Gc.den));

% Forward Euler
Gd1.sym = subs(Gc.sym,s,(z-1)/h);
[Gd1.num,Gd1.den] = numden(Gd1.sym);
Gd1.tf = tf(sym2poly(Gd1.num),sym2poly(Gd1.den),h);

% Backward Euler
Gd2.sym = subs(Gc.sym,s,(z-1)/(h*z));
[Gd2.num,Gd2.den] = numden(Gd2.sym);
Gd2.tf = tf(sym2poly(Gd2.num),sym2poly(Gd2.den),h);

% Tustin
Gd3.sym = subs(Gc.sym,s,2*(z-1)/(h*(z+1)));
[Gd3.num,Gd3.den] = numden(Gd3.sym);
Gd3.tf = tf(sym2poly(Gd3.num),sym2poly(Gd3.den),h);

% Tustin with pre-warping
wlp   = 2/h*tan(wlp*h/2);
flp   = wlp/(2*pi);
wp{1} = 2/h*tan(wp{1}*h/2);
wz{1} = 2/h*tan(wz{1}*h/2);
wp{2} = 2/h*tan(wp{2}*h/2);
wz{2} = 2/h*tan(wz{2}*h/2);
Gc2.sym = (wp{1}^2)/(wz{1}^2)*(s^2 + 2*bz{1}*wz{1}*s + wz{1}^2) / ...
                              (s^2 + 2*bp{1}*wp{1}*s + wp{1}^2) * ...
          (wp{2}^2)/(wz{2}^2)*(s^2 + 2*bz{2}*wz{2}*s + wz{2}^2) / ...
                              (s^2 + 2*bp{2}*wp{2}*s + wp{2}^2) * ...                     
          (Kp*(1/(2*pi*fd)*s + 1 + (2*pi*fi)/s)*(2*pi*flp)^2 / ...
          (s^2 + 2*1*(2*pi*flp )*s + (2*pi*flp)^2));   

[Gc2.num,Gc2.den] = numden(Gc2.sym);
Gc2.tf = tf(sym2poly(Gc2.num),sym2poly(Gc2.den));

Gd4.sym = subs(Gc2.sym,s,2*(z-1)/(h*(z+1)));
[Gd4.num,Gd4.den] = numden(Gd4.sym);
Gd4.tf = tf(sym2poly(Gd4.num),sym2poly(Gd4.den),h);

% Matched Pole Zero
Gd5.tf = c2d(Gc.tf,h,'matched');

% Frequency spacing
f = logspace(0,log10(5000),3000);
wc = 2*pi*f;
wd = wc(wc < 2*pi*2500);

% Plot
hold all;
xlim([1 3500]);
grid on;
legend({'Continuous-time','Forward Euler','Backward Euler', ...
    'Tustin','Tustin with pre-warping','Matched Pole Zero'});