Discrete-time systems and the z-transform

Inside a robot's controller there are no differential equations, only loops of code that run once per sample: read the sensors, update some variables, write the outputs. Such code implements a difference equation, the discrete-time counterpart of a differential equation, and its behaviour can be analysed with tools that parallel the continuous ones exactly. The z-transform plays the role of the Laplace transform, the unit delay \( z^{-1} \) plays the role of the integrator 1/s, the unit circle plays the role of the imaginary axis, and the inside of the circle plays the role of the left half-plane. This lesson builds that toolkit, then tackles the practical problem every embedded engineer faces: converting a continuous design, such as a filter or controller drawn in s, into code. Different conversion methods, forward Euler, backward Euler, Tustin and the matched exact method, have different accuracy and stability properties, and choosing badly can make a perfectly good design go unstable on the robot.

Difference equations

A linear time-invariant discrete system relates its current output to current and past inputs and past outputs:

This is directly code: keep the last few inputs and outputs in variables, and each sample compute the weighted sum. If all the a coefficients are zero the output depends only on inputs, the impulse response has finite length, and the system is an FIR filter, like the moving average of lesson 3. If any a is nonzero the output feeds back on itself, the impulse response goes on forever, and the system is recursive, or infinite impulse response (IIR). The first-order low-pass filter \( y[n] = \alpha y[n-1] + (1 - \alpha)x[n] \) used for sensor smoothing in ROB-220 is the simplest IIR system. Its impulse response is \( (1 - \alpha)\alpha^n \), which decays if \( |\alpha| < 1 \) and explodes if \( |\alpha| > 1 \).

The z-transform

The z-transform of a sequence x[n] that is zero for n < 0 is

where z is a complex variable. It is the discrete analogue of the Laplace transform; in fact, for a sampled signal it equals the Laplace transform with the substitution \( z = e^{sT} \). Its key property is the delay rule: delaying a sequence by one sample multiplies its transform by \( z^{-1} \),

so \( z^{-1} \) is the unit delay, the memory element of digital systems, drawn as a block in digital filter diagrams. Applying the transform to the general difference equation with zero initial conditions gives the discrete transfer function

which can be read straight off the code and vice versa. A few transform pairs suffice for most work: the unit impulse \( \delta[n] \) transforms to 1, the unit step to \( z/(z - 1) \), and the geometric sequence \( a^n \) to \( z/(z - a) \). Convolution of sequences becomes multiplication of transforms, as in continuous time.

Poles, the unit circle and stability

The mapping \( z = e^{sT} \) carries the s-plane to the z-plane in a way that explains everything about discrete stability. The imaginary axis \( s = j\omega \) maps to the unit circle \( |z| = 1 \); the left half-plane maps to the inside of the circle; the right half-plane maps to the outside. A continuous pole at \( s = -a \) becomes a discrete pole at \( z = e^{-aT} \), a real number between 0 and 1. A complex pair with damped frequency \( \omega_d \) becomes a pair at angle \( \pm\omega_dT \). Hence the discrete stability rule: an LTI discrete system is stable if and only if all its poles lie strictly inside the unit circle. Poles near the origin give fast decay; poles near the circle give slow decay; a pole at z = 1 is a discrete integrator (accumulator); a pole on the negative real axis gives a sequence that alternates sign every sample, a behaviour with no continuous counterpart, often seen as chatter in poorly tuned digital loops.

Frequency response follows the same pattern: the discrete frequency response is \( H(e^{j\omega T}) \), the transfer function evaluated around the unit circle. Going once round the circle, from angle 0 to π, covers frequencies from 0 to the Nyquist frequency \( f_s/2 \); going further simply repeats, which is the z-domain view of aliasing.

From s to z: discretizing continuous designs

Most filters and controllers are designed in continuous time, then converted to code. The conversion approximates differentiation or integration over one sample. Four methods are common.

• Forward Euler: \( \dot{y} \approx (y[n+1] - y[n])/T \), equivalent to \( s \to (z - 1)/T \). It is the simplest, but it maps part of the stable left half-plane outside the unit circle, so a stable continuous design can become unstable if T is too large relative to its fastest time constant. For a pole at −a it requires \( aT < 2 \).
• Backward Euler: \( \dot{y} \approx (y[n] - y[n-1])/T \), equivalent to \( s \to (z - 1)/(zT) \). It maps the whole left half-plane inside the circle, so it never destabilizes a stable design, but it adds phase lag and damping.
• Tustin (bilinear transform): \( s \to \dfrac{2}{T}\dfrac{z - 1}{z + 1} \), which corresponds to trapezoidal integration. It maps the left half-plane exactly onto the inside of the circle and preserves stability, with better accuracy than either Euler method. It compresses the frequency axis, warping high frequencies, so filters with critical frequencies, such as notches, are designed with frequency prewarping, which makes the discrete response exact at one chosen frequency.
• Zero-order-hold (exact) discretization: computes the discrete system that gives exactly the same samples as the continuous system driven through a zero-order hold, mapping each pole by \( z = e^{sT} \). It is the standard for discretizing plant models for digital control design, and it is what simulation libraries do for state-space models (lesson 12).

import numpy as np
a = 100.0
for T in [0.001, 0.01, 0.025]:
    aT = a*T
    euler = 1 - aT
    tustin = (2 - aT)/(2 + aT)
    exact = np.exp(-aT)
    print(f"T = {T*1000:5.1f} ms: pole euler {euler:7.4f}  tustin {tustin:7.4f}  exact {exact:7.4f}")

T =   1.0 ms: pole euler  0.9000  tustin  0.9048  exact  0.9048
T =  10.0 ms: pole euler  0.0000  tustin  0.3333  exact  0.3679
T =  25.0 ms: pole euler -1.5000  tustin -0.1111  exact  0.0821

At 1 ms all three agree; at 10 ms forward Euler places the pole at zero, a filter with no memory at all; at 25 ms it puts the pole at −1.5, outside the unit circle, so the 'filter' is unstable, while Tustin and the exact method stay inside. The lesson for embedded code: either sample fast relative to the dynamics, or use a method that preserves stability.