Digital Control
1. At a glance
Digital control is control law implemented on a clocked digital processor — an MCU, DSP, FPGA, or PLC — sampling sensor signals at discrete times t = kT_s, computing an updated command, and writing the result to actuators before the next sample. It is the bridge between the textbook continuous-time control engineer learns first (transfer functions G(s) in the Laplace domain, root-locus, Bode plots, Nyquist) and the firmware that actually runs on real silicon in real products.
Every modern control loop is digital. Analog op-amp PID and lead-lag networks survive only in niches: very-high-bandwidth RF AGC, some audio analog effects, sub-µs laser-diode current servos, and legacy industrial process equipment built before the 1990s. Everything else — motor drives, robot joint loops, aircraft fly-by-wire, missile autopilots, vehicle ABS/ESC/cruise, chemical-plant PID, drone attitude, drive-by-wire, switching-power-supply voltage loops, hard-disk head positioning — is sampled-data control running on a microcontroller, DSP, FPGA, or PLC. Analog continues to appear inside loops (the gate driver and current sense in a motor drive are analog), but the outer regulating brain is digital.
The discipline lives at a precise boundary. Above it sits continuous control theory (the s-plane, Bode plots, classical and modern state-space methods — see [[Engineering/classical-control]] and [[Engineering/state-space-methods]]). Below it sits embedded software (RTOS, ISR latency, fixed-point math, deterministic timing — see [[Engineering/microcontrollers]] and [[Engineering/realtime-embedded]]). Digital control is the discretisation step in between, and the single largest source of field-deployed control bugs. Aliasing because somebody forgot the anti-alias filter, integrator wind-up because the actuator saturated, limit cycles from coefficient quantisation, missed deadlines because an unrelated ISR went long, derivative kick on a setpoint step that broke a coupling — all are digital-control failure modes, not theory failures.
Three issues dominate practical work. Sample rate selection (too slow loses phase margin; too fast wastes MCU cycles and amplifies quantisation noise). Discretisation method (Tustin vs ZOH vs backward Euler — which preserves which property of the continuous design). Computational delay (the half-sample of ZOH lag plus the actual compute time eats phase margin and is often the reason a loop that simulated fine oscillates on hardware).
2. First principles
Sampling: from continuous to discrete time
A continuous signal x(t) is sampled at uniform intervals T_s (sample period, units: seconds). The sample rate is f_s = 1/T_s (Hz) or ω_s = 2π/T_s (rad/s). The resulting discrete sequence is x[k] = x(kT_s).
Between samples the digital-to-analog converter holds the most recent command constant — this is the zero-order hold (ZOH). The composite plant seen by the discrete controller is therefore “ZOH + continuous plant + sampler”, not the bare continuous plant. This trio is what discretisation captures.
The sampling theorem (Nyquist–Shannon, 1948)
A bandlimited continuous signal with maximum frequency component ω_max can be reconstructed perfectly from its samples iff ω_s > 2·ω_max. The bound ω_N = ω_s/2 is the Nyquist frequency. Any energy above ω_N folds back into the [0, ω_N] band as aliases: a sinusoid at frequency ω > ω_N appears, after sampling, as a sinusoid at |ω − n·ω_s| for the integer n that brings the result into [0, ω_N]. Aliasing is irreversible — once folded, an alias is indistinguishable from a real signal at the lower frequency.
The fix is an anti-alias filter (AAF): an analog low-pass filter ahead of the ADC with corner frequency comfortably below ω_N (typically ω_c ≈ 0.4·ω_s for a 2nd–4th-order filter, lower for steeper roll-off). The AAF removes the high-frequency content that would otherwise alias; once a signal is in the digital domain the AAF cannot be added retroactively.
Z-transform
The Z-transform plays the role for sampled-data systems that the Laplace transform plays for continuous systems. For a discrete sequence x[k]:
X(z) = Σ_{k=0}^{∞} x[k] · z^{-k}
The complex variable z lives on the complex plane; z = e^{jωT_s} traces the unit circle as ω ranges over [0, ω_s). The mapping from s-plane to z-plane induced by sampling is z = e^{sT_s}: the LHP of s (Re(s) < 0, stable continuous poles) maps to the interior of the unit circle (|z| < 1, stable discrete poles). The s-plane jω axis maps to the unit circle |z| = 1. The origin s = 0 maps to z = 1.
Key properties:
- Linearity: Z{a·x[k] + b·y[k]} = a·X(z) + b·Y(z).
- Time-shift: Z{x[k − n]} = z^{-n} · X(z). The unit delay z^{-1} is the elementary building block of every digital filter.
- Convolution: Z{x[k] * h[k]} = X(z) · H(z). A linear time-invariant discrete system has transfer function H(z) = Y(z)/X(z).
- Initial value theorem: x[0] = lim_{z→∞} X(z).
- Final value theorem: x[∞] = lim_{z→1} (1 − z^{-1}) · X(z), valid when poles of (1 − z^{-1})X(z) are strictly inside |z| < 1.
Stability rule: a discrete LTI system is stable iff all poles of H(z) lie inside the open unit disc |z| < 1. Continuous-time stability (poles in the open LHP) and discrete-time stability look different but mean the same physical thing: bounded input produces bounded output.
Common Z-transform pairs (unit-pulse δ[k] at z = 1):
| x[k] (k ≥ 0) | X(z) | ROC |
|---|---|---|
| δ[k] | 1 | all z |
| 1 (unit step) | z/(z − 1) | |z| > 1 |
| k | T_s · z / (z − 1)² | |z| > 1 |
| a^k | z/(z − a) | |z| > |a| |
| k · a^k | a·z/(z − a)² | |z| > |a| |
| sin(ωkT_s) | z·sin(ωT_s) / (z² − 2z·cos(ωT_s) + 1) | |z| > 1 |
| cos(ωkT_s) | z·(z − cos(ωT_s)) / (z² − 2z·cos(ωT_s) + 1) | |z| > 1 |
Difference equations and the discrete transfer function
The discrete equivalent of a continuous differential equation is a difference equation:
a₀·y[k] + a₁·y[k−1] + … + a_n·y[k−n] = b₀·u[k] + b₁·u[k−1] + … + b_m·u[k−m]
Taking Z-transforms and rearranging gives the rational transfer function:
H(z) = Y(z)/U(z) = (b₀ + b₁·z^{-1} + … + b_m·z^{-m}) / (a₀ + a₁·z^{-1} + … + a_n·z^{-n})
This is precisely the form an MCU computes: tap-delayed inputs and outputs combined with multiply-accumulates. The discrete frequency response is obtained by evaluating H(z) on the unit circle:
H(e^{jωT_s}) for ω ∈ [0, ω_s/2]
unlike continuous Bode plots which extend to infinite frequency, the discrete Bode is periodic with period ω_s and bandlimited to ω_N — there is no information above Nyquist.
3. Practical math / design equations
Discretisation: turning G(s) into G(z)
Six standard methods. They differ in how they approximate the integral s^{-1} or the differential operator s.
| Method | Mapping (substitute into G(s)) | Pros | Cons | When to use |
|---|---|---|---|---|
| Zero-Order Hold (ZOH) | (1 − z^{-1}) · Z{G(s)/s}, evaluated as a step-invariant transform | Exactly matches the step response of plant + ZOH; physical — what a DAC really does | Frequency response distorted near Nyquist | Discretising a plant model for digital-control design; default for plant identification |
| First-Order Hold (FOH) | Ramp-invariant transform | Better for ramp-like inputs | More algebra; rarely used in practice | Specialty trajectory-tracking work |
| Tustin / bilinear | s ← (2/T_s) · (z − 1)/(z + 1) | Maps LHP→unit disc one-to-one; stable continuous → stable discrete; preserves frequency response qualitatively | Frequency warping: ω_c → (2/T_s)·tan(ω_c·T_s/2) | Discretising a controller designed in continuous time; most common choice |
| Tustin with pre-warping | Same as Tustin, scaled so one chosen frequency ω_p matches | Bode magnitude and phase exact at ω_p | Only one frequency is matched exactly | Notch filters, narrow-band controllers where ω_p is the critical frequency |
| Forward (rectangular) Euler | s ← (z − 1)/T_s | Trivial coefficient computation; cheap | Stable continuous can become unstable discrete if T_s too large | Quick prototyping; do not ship |
| Backward Euler | s ← (z − 1)/(z · T_s) | Always preserves stability (stable continuous → stable discrete) | Adds half-sample phase lag at all frequencies | Conservative; common for integrators |
| Matched pole-zero | each s-plane pole/zero p_i → z-plane pole/zero z_i = e^{p_i · T_s}; gain matched at DC or another frequency | Preserves pole/zero structure of the continuous design | Hand work for each pole; not automated in some tools | High-Q resonant controllers, observers |
In MATLAB, the function c2d(sys_continuous, Ts, method) accepts 'zoh', 'foh', 'tustin', 'matched', 'impulse'. In Python, scipy.signal.cont2discrete((num, den), dt, method) accepts 'zoh', 'bilinear' (alias for Tustin), 'euler', 'backward_diff'. The python-control package mirrors MATLAB’s API: control.c2d(sys, Ts, method).
Tustin pre-warping corrects the frequency warping at one chosen design frequency ω_p:
s ← (ω_p / tan(ω_p · T_s / 2)) · (z − 1)/(z + 1)
This is the right discretisation for a continuous notch designed at exactly ω_p: the digital notch sits at exactly ω_p, not at a slightly mapped neighbour. Outside ω_p the warping is still present, but the dominant feature is correct.
Discrete PID — positional vs velocity form
A continuous PID is u(t) = K_p · [e(t) + (1/T_i)·∫e(τ)dτ + T_d·de/dt], where e = r − y. Three discretisation choices, each yielding a different algorithm.
Positional (absolute) form, backward Euler:
- P term: u_p[k] = K_p · e[k]
- I term: u_i[k] = u_i[k−1] + K_p · (T_s/T_i) · e[k]
- D term (with first-order filter at N·1/T_d, typical N = 8…20 to limit noise): u_d[k] = (T_d / (T_d + N·T_s)) · u_d[k−1] − (K_p · T_d · N / (T_d + N·T_s)) · (y[k] − y[k−1])
Total: u[k] = u_p[k] + u_i[k] + u_d[k].
Note the derivative is taken on −y[k], not on the error e[k] — this is derivative-on-PV and avoids “derivative kick” when the setpoint is stepped (a step in r produces an infinite continuous derivative; differentiating only y avoids the spike).
Velocity (incremental) form:
Δu[k] = u[k] − u[k−1] = K_p · {(e[k] − e[k−1]) + (T_s/T_i) · e[k] + (T_d/T_s) · (−y[k] + 2·y[k−1] − y[k−2])}
Then u[k] = u[k−1] + Δu[k] outside the algorithm (with saturation logic on u).
Trade-offs:
| Aspect | Positional | Velocity |
|---|---|---|
| Bumpless transfer (auto ↔ manual) | Requires tracking | Trivial — u[k−1] already current |
| Anti-windup | Explicit clamp + back-calc on u_i | Implicit — saturating u stops Δu from accumulating |
| Setpoint changes | Need derivative-on-PV; need setpoint weighting | Same |
| Plant offsets, integrator drift | Cleaner — integrator state is bounded explicitly | Δu state drifts with integer overflow if not careful |
| Two-degree-of-freedom (β, γ weighting) | Standard in process-control textbooks (ISA form) | Possible but uglier |
Setpoint weighting (two-degree-of-freedom PID):
u[k] = K_p · [(β·r[k] − y[k]) + (T_s/T_i) · Σ(r − y) + (T_d/T_s) · (γ·r[k] − y[k]) − …]
β and γ ∈ [0, 1] reduce the influence of the setpoint on the P and D terms. β ≈ 0.3 and γ = 0 are common — they soften step response without changing disturbance rejection (which is dominated by the y feedback). The ISA-PID block in most DCS systems (Emerson DeltaV, Honeywell Experion, Yokogawa Centum) implements this form.
Anti-windup
When the actuator saturates (PWM hits 100 %, valve hits fully open, motor drive hits current limit), the integrator keeps accumulating error — and when the saturation eventually clears, the integral state has wound up beyond reason and the loop overshoots wildly. Three standard mitigations:
-
Clamping (conditional integration): stop integrating when u is saturated and the error sign would drive it further into saturation.
if (!(u_sat && (sign(e) == sign(u_sat)))) u_i += Kp * Ts/Ti * e; -
Back-calculation: drive the integrator state toward what would make u = u_sat. Add a term K_aw · (u_sat − u_unsat) to du_i/dt, with K_aw typically 1/T_t where T_t = √(T_i · T_d) (Åström’s rule of thumb).
-
Conditional integration with deadband: simpler clamp variant; ignore integration when |e| < threshold.
Skeleton (positional, derivative-on-PV, back-calculation anti-windup):
static float u_i = 0.0f, y_prev = 0.0f, u_d = 0.0f;
const float Kp = 0.5f, Ti = 0.2f, Td = 0.0f, N = 10.0f;
const float Ts = 0.001f; // 1 kHz
const float U_MAX = +1.0f, U_MIN = -1.0f;
const float Kaw = 1.0f / sqrtf(Ti * (Td + 1e-9f));
float pid_step(float r, float y) {
float e = r - y;
// P
float u_p = Kp * e;
// D on PV, with filter
float a = Td / (Td + N * Ts);
float b = Kp * Td * N / (Td + N * Ts);
u_d = a * u_d - b * (y - y_prev);
y_prev = y;
// Pre-saturation
float u_unsat = u_p + u_i + u_d;
// Saturate
float u = u_unsat;
if (u > U_MAX) u = U_MAX;
else if (u < U_MIN) u = U_MIN;
// Integral with back-calculation anti-windup
u_i += Kp * (Ts / Ti) * e + Kaw * (u - u_unsat) * Ts;
return u;
}Sample-rate selection
Rules of thumb collected from Åström/Wittenmark, Franklin/Powell/Workman, and field practice:
| Application | f_s / closed-loop bandwidth ratio | Comment |
|---|---|---|
| Process control (slow chemical loops) | 5–10× | Slow plants, samples of seconds; PLCs at 100 ms—1 s |
| Servo positioning, motor speed loops | 20–30× | Typical: 1–10 kHz current loop, 1 kHz speed loop |
| Motor current (FOC) | 50–100× | 10–20 kHz inner current loop on TI C2000 / STM32G4 |
| Audio DSP | ≫ 100× | 48 kHz / 96 kHz / 192 kHz — bandwidth set by source, not loop dynamics |
| Power converters (DC-DC, inverters) | 5–10× of switching frequency | Sample synced to PWM; 50–500 kHz typical |
| Safety-critical (flight control, brake-by-wire) | 50–200× | DO-178C / ISO 26262 require ample margin; redundant sampling |
Lower bound: at fewer than 5× CL bandwidth the half-sample ZOH lag dominates and the digital loop has much less phase margin than the continuous design implied. Upper bound: at more than 200× CL bandwidth quantisation noise and finite-precision coefficient sensitivity dominate; pole locations migrate to or outside the unit circle.
Quantisation and ADC resolution
For an N-bit ADC with full-scale range V_FS, the quantisation step is q = V_FS / 2^N. Assuming uniform quantisation noise:
- Signal-to-quantisation-noise ratio (SQNR): SQNR = 6.02·N + 1.76 dB (for a full-scale sinusoid).
- Effective Number of Bits (ENOB): ENOB = (SINAD − 1.76 dB) / 6.02. Always ENOB < N due to thermal noise, integral nonlinearity (INL), differential nonlinearity (DNL), and clock jitter.
| Part | Type | N | ENOB | Notes |
|---|---|---|---|---|
| STM32G4 on-chip | SAR | 12-bit | ~10.5 @ 4 MSPS | Free with the MCU |
| STM32H7 on-chip | SAR | 16-bit (with oversampling) | ~13.5 | 3.6 MSPS native, 16-bit via hardware oversample |
| TI ADS131M08 | Σ-Δ, 8-ch simultaneous | 24-bit | ~21 @ 8 kSPS | Industrial energy / motor sense |
| TI ADS127L11 | Σ-Δ | 24-bit | ~23 @ 4 kSPS | Precision instrumentation |
| Microchip MCP3461R | Σ-Δ | 24-bit | ~21 @ 122 SPS | Energy metering, weighing |
| Analog Devices LTC2387 | SAR pipelined | 18-bit | ~17.5 @ 15 MSPS | Industrial imaging, ultrasound |
| TI ADS5263 | Pipelined | 16-bit | ~13.5 @ 100 MSPS | Software-defined radio |
For control, an ENOB ≥ 10 over the operating range of the controlled variable is the practical lower bound — below that, limit cycles appear (the loop oscillates by ±1 LSB of the ADC because the controller cannot resolve small errors). The fix is either more ADC bits, deliberate dither, or accepting the LSB-scale oscillation.
Computational and ZOH delay
The ZOH itself introduces a half-sample of phase lag: φ_ZOH(ω) = −ω·T_s/2 radians. At the closed-loop crossover ω_c this subtracts from the phase margin budget computed in continuous time:
ΔPM = (180°/π) · ω_c · T_s / 2
Example: f_s = 1 kHz (T_s = 1 ms) and ω_c = 100 rad/s ⇒ ΔPM = 2.86°. At f_s = 10 kHz ⇒ ΔPM = 0.29°. At f_s = 100 Hz with the same plant ⇒ ΔPM = 28.6° — eats most of a typical 45° design margin.
Computational delay is additional: if the ADC is read at the start of the cycle, the control is computed, and the result is written at the next cycle boundary, there is one full T_s of dead time on top of the half-sample ZOH lag. Total: −1.5·ω·T_s/2 = −0.75·ω·T_s. If the result is written at the end of the same cycle as soon as the compute finishes, the total is intermediate. The discipline of writing PWM at the next PWM update event (deterministic) trades a bit of average lag for predictability — and predictability is what STA-equivalent timing analysis of the loop relies on.
4. Reference data
Implementation architectures
| Architecture | Determinism | Typical f_s | Use case |
|---|---|---|---|
| Bare-metal MCU, timer ISR | Hard | up to ~50 kHz on M4, 200 kHz+ on M7 with TCM | Motor drives, regulator loops, PMSM FOC |
| MCU + DMA double-buffer | Hard | up to ADC native rate (MHz) | High-rate sensor sampling, batch processing |
| RTOS periodic task | Soft-hard | up to ~1–10 kHz depending on RTOS jitter | General embedded control where priority structure matters |
| Cooperative state machine | Soft | up to ~1 kHz | Simple firmware without RTOS |
| DSP (TI C2000, ADI SHARC) | Hard | 100–200 kHz control, MHz signal proc | High-end motor drives, power electronics |
| FPGA control (Xilinx Kintex / Lattice ECP5) | Hard, sub-µs jitter | 1 MHz+ | Multi-axis motion, ultra-low-latency servos, hardware-in-the-loop |
| CompactRIO / cRIO + FPGA personality | Hard | 1 kHz–100 kHz | NI lab / aerospace test rigs; FPGA personality + RT host |
| dSPACE MicroAutoBox / SCALEXIO | Hard | 1–10 kHz typical | Automotive rapid prototyping, autocode from Simulink |
| PLC (Allen-Bradley ControlLogix, Siemens S7-1500) | Soft (ms scan) | 10–100 ms cycle | Process industries, factory automation; IEC 61131-3 ST PID blocks |
| OPAL-RT HYPERSIM, RT-LAB | Hard | 25–100 kHz | Real-time power-system simulation, HIL |
MCU peripheral primitives for control loops
| Peripheral | Function | Examples |
|---|---|---|
| Advanced timer / motor-control timer | Center-aligned PWM, ADC trigger, dead-time, fault input | STM32 TIM1/TIM8, TI C2000 EPWM, NXP S32K eTimer / FlexPWM |
| Synchronised ADC | Sample at PWM mid-point (zero-current ripple) | STM32 ADCs in dual/triple mode + injected channels, C2000 ADC SOC |
| Quadrature encoder interface | Hardware decode + index capture | STM32 TIM in encoder mode, C2000 eQEP, AVR mega328 (no native; uses ext ints) |
| DMA | Move ADC results to RAM without CPU intervention | All Cortex-M, C2000, ESP32 |
| Hardware divide / FPU | Single-cycle FP for biquads, transforms | M4F/M7F single-precision, C2000 FPU32 + TMU |
| CAN-FD / FlexCAN | Networked control loops (vehicle, robotics) | NXP S32, STM32 FDCAN, Renesas RH850 |
| Lockstep cores | Safety-critical fault detection | Infineon AURIX TC3xx, NXP S32K3, Renesas RH850/U2A |
Specific control SoCs in production
| Part | Use case | Notable features |
|---|---|---|
| TI TMS320F2837xD (C2000) | Motor / inverter control | Dual C28x @ 200 MHz + dual CLA real-time co-processors, TMU (trig accel), 16-bit Σ-Δ ADCs, 200 ps EPWM |
| TI TMS320F28004x | Mid-tier motor control | Single C28x @ 100 MHz, CLA, 16-bit ADCs, USD 4–8 |
| STM32G474 | Digital power and motors | Cortex-M4F @ 170 MHz, 12-bit ADCs @ 4 MSPS, HRTIM with 184 ps PWM resolution, 4 op-amps + 4 comparators on-die |
| STM32H723 | High-end motor / servo | Cortex-M7 @ 550 MHz, 16-bit ADCs (oversampled), L1 cache, MDMA |
| NXP S32K344 | Automotive ASIL-D | Lockstep Cortex-M7 ×3, ECC RAM/Flash, FlexCAN-FD, HSE security; AEC-Q100 |
| ADI ADSP-CM41x | Digital power / motors | Cortex-M4 + Σ-Δ ADCs on-die; popular in solar inverters |
| Renesas RX72T | Industrial servo | RXv3 core @ 200 MHz, 12-bit ADC at 1 MSPS, hardware divider for sin/cos |
| Microchip dsPIC33CK | Motor and digital power | dsPIC core with DSP MAC, fast 12-bit ADCs, high-resolution PWM |
| Xilinx Zynq UltraScale+ MPSoC | Multi-axis motion + Linux UI | Quad-A53 + dual-R5 + FPGA fabric — control runs on R5 or FPGA |
Tools and toolchains
| Tool | Function |
|---|---|
| MATLAB Control System Toolbox | c2d, bode, nichols, nyquist, margin, rlocus for design and discretisation |
| MATLAB/Simulink + Embedded Coder | Model-based design, autocode generation to C for STM32, C2000, AURIX, S32K, MicroAutoBox |
Python scipy.signal, python-control | Open-source equivalent of MATLAB CST |
| PLECS Coder (Plexim) | Power-electronics circuit + control codegen for C2000, STM32, FPGA |
| PSCAD / EMTDC | Electromagnetic-transient power-system simulation with controller models |
| Tracealyzer (Percepio), SystemView (Segger) | Real-time visualisation of task and ISR scheduling jitter |
| Cortex-M DWT cycle counter | Sub-cycle profiling on-target |
| LDRA Testbed, Polyspace, Coverity | MISRA C / CERT C / ISO 26262 static analysis |
| TI Code Composer Studio, IAR EWARM, Keil MDK, STM32CubeIDE | Target IDEs |
| OPAL-RT HYPERSIM, NI VeriStand, dSPACE ControlDesk | Hardware-in-the-loop platforms |
| Speedgoat real-time targets | Simulink-native HIL hardware |
5p. Theory — discretisation derivations
ZOH equivalent
Given continuous plant G(s), preceded by a DAC that holds the input constant for T_s, the discrete equivalent is:
G_ZOH(z) = (1 − z^{-1}) · Z{ G(s)/s }
The inner expression Z{G(s)/s} is the Z-transform of the step response of G(s) sampled every T_s; multiplying by (1 − z^{-1}) implements the unit-step-input/unit-pulse-input difference. This makes ZOH step-invariant: the discrete system has exactly the same response to a step input as the continuous system + DAC, at every sample instant. Between samples the digital model says nothing — the real DAC just holds.
For a single pole G(s) = K/(s + a): G_ZOH(z) = K·(1 − e^{−aT_s}) / (a·(z − e^{−aT_s})).
Bilinear (Tustin) transform
The trapezoidal-rule approximation to ∫: y(t) ≈ y(t − T_s) + (T_s/2)·(u(t) + u(t − T_s)). Z-transforming, the integrator 1/s maps to:
1/s ≈ (T_s/2) · (z + 1)/(z − 1)
equivalently s ≈ (2/T_s)·(z − 1)/(z + 1). This is one-to-one between the LHP and the unit disc — every stable G(s) maps to a stable G(z) and vice-versa. The cost is frequency warping: a continuous feature at ω_c maps to a discrete feature at ω_d = (2/T_s) · arctan(ω_c · T_s / 2). For ω_c·T_s ≪ 1, ω_d ≈ ω_c; near Nyquist the mapping pulls features toward ω_s/2.
Forward and backward Euler
Forward Euler: dy/dt ≈ (y[k+1] − y[k])/T_s, so s ≈ (z − 1)/T_s. The LHP maps to the half-plane Re(z) < 1, which does not fit inside the unit disc — a stable G(s) with a pole at s = −1/T_s maps to z = 0 (still stable) but a pole at s = −3/T_s maps to z = −2 (unstable). The method fails when sample rate is not fast enough relative to plant dynamics.
Backward Euler: dy/dt ≈ (y[k] − y[k−1])/T_s, so s ≈ (z − 1)/(z · T_s). The LHP maps inside |z − 1/2| < 1/2 — a circle strictly inside the unit disc — so every stable G(s) is stable when discretised. The cost is half-sample phase lag at all frequencies.
6p. Worked examples
Example A — Speed loop on an STM32G474
Plant: DC motor, transfer function from voltage to speed G(s) = K_m / (τ·s + 1), with K_m = 100 (rad/s)/V and τ = 0.05 s. Continuous PI design target: bandwidth ω_c = 30 rad/s, PM 60°. Designed continuous controller: K_p = 0.5, T_i = 0.1 s.
Sample at f_s = 1 kHz (T_s = 1 ms). Use Tustin discretisation.
For the I term, Tustin gives 1/(s·T_i) → (T_s/(2·T_i)) · (1 + z^{-1})/(1 − z^{-1}). With T_s/(2·T_i) = 0.005:
u_i[k] = u_i[k−1] + 0.005 · K_p · (e[k] + e[k−1])
Phase lag from ZOH at ω = 30 rad/s: φ_ZOH = −30 · 0.001 / 2 · (180°/π) = −0.86°. Negligible vs the 60° margin budget. The discretised loop will look indistinguishable from the continuous design in simulation.
If sampling instead at f_s = 100 Hz (T_s = 10 ms), φ_ZOH = −8.6° — the margin drops from 60° to ~51°, still acceptable. At f_s = 30 Hz (T_s = 33 ms), φ_ZOH ≈ −28° — half the margin gone, and a closed-loop step response will overshoot considerably more than designed.
Example B — Servo lead-lag + notch at 80 Hz, discretised at 10 kHz
Continuous notch at ω_n = 2π·80 = 503 rad/s, depth 30 dB, Q = 10:
N(s) = (s² + (ω_n/(Q·10^{30/20}))·s + ω_n²) / (s² + (ω_n/Q)·s + ω_n²)
Discretise with Tustin with pre-warping at ω_p = ω_n. Pre-warping substitutes:
α = ω_n / tan(ω_n · T_s / 2) = 503 / tan(503 · 1e−4 / 2) = 503 / tan(0.02515) = 503 / 0.02515 ≈ 20000
(very close to 2/T_s = 20000 since ω_n·T_s/2 = 0.025 is small). Substituting s ← α·(z − 1)/(z + 1) into N(s) and clearing denominators gives a discrete biquad b0 + b1·z^{-1} + b2·z^{-2} / (1 + a1·z^{-1} + a2·z^{-2}) whose magnitude response has a notch at exactly 80 Hz (the pre-warping target) — without pre-warping it would sit at ~79.7 Hz, a 0.4 % shift inaudible for most servo work but visible in lab measurements.
Implementation: Direct Form II Transposed (DF2T) biquad. Of the four canonical biquad forms (DF1, DF2, DF1T, DF2T), DF2T is preferred for fixed-point implementation because the two state variables are stored in a high-precision accumulator (32-bit or 64-bit) between samples — coefficient quantisation effects are minimised. CMSIS-DSP arm_biquad_cascade_df2T_f32() implements this form natively.
Quantising the coefficients to Q1.15 (16-bit fixed-point with 1 sign bit and 15 fractional bits, range [−1, +1)) introduces ±2^{-16} ≈ ±15 ppm per coefficient. For a high-Q notch this can shift the notch frequency by 0.1 Hz and reduce notch depth by 1–2 dB. Float coefficients on an M4F eliminate this entirely; the choice between fixed and float is BOM cost on bare M4 vs M4F (negligible today).
Example C — Computational delay handling, 4 kHz control loop
T_s = 250 µs. Profile (DWT cycle counter, STM32G474 @ 170 MHz):
- Wait for ADC EOC: ~50 µs (8500 cycles)
- Compute Park/Clarke transform, two current PIs, inverse Park, SVPWM: 80 µs (13600 cycles)
- Write to TIM1 CCR registers: 1 µs
Total: 131 µs after the sample instant. The PWM update happens at the next TIM1 update event (synchronous with the sampling), so effective delay is 250 µs (one full sample). Discrete model includes one extra unit delay z^{-1} on the actuator path.
Phase lag at ω_c = 2π · 200 Hz (current loop crossover) = 1257 rad/s: φ_delay = −1257 · 250e−6 · (180°/π) = −18.0°. Combined with the ZOH half-sample (φ_ZOH = −9.0°), total is −27° eaten by digital implementation alone — must be subtracted from the continuous design’s phase margin.
Mitigation that is often not worth it but is sometimes essential: apply control synchronously to the current sample instead of waiting one full cycle (“skid”). Sample at t = kT_s, compute during kT_s + ε, write to PWM CCR at kT_s + (T_s − Δ) for Δ ≈ 1 µs before the PWM reload. This recovers approximately one T_s of delay but requires absolutely deterministic compute time — any overrun and the previous frame’s PWM is held one extra cycle, which is much worse than the deterministic-delay scheme. Use only with rigorously bounded compute and a watchdog.
7p. Edge cases and gotchas
Missing anti-alias filter. A 50/60 Hz mains hum aliasing into a 2 kHz current loop sampled at 4 kHz appears at 50 Hz (in-band) and corrupts the loop. The bench symptom is mysterious low-frequency oscillation that disappears when mains power is on a UPS. The cure is an analog RC filter or differential active filter ahead of every ADC channel that doesn’t already have a Σ-Δ ADC (Σ-Δ have built-in oversampling decimation that acts as a sharp digital AAF, but the analog modulator still benefits from a simple RC roll-off).
Sample jitter from ISR contention. If the ADC-EOC ISR can be delayed by another ISR’s execution, the effective T_s is not constant — and a constant-T_s discretisation assumes it is. Symptom: spurs at the contention rate (e.g. SysTick at 1 kHz interferes with a 10 kHz loop and produces a 1 kHz spur in motor current). Fix: priority the control ISR above all else, or use a synchronous DMA path that does not depend on the CPU.
Coefficient quantisation in fixed-point. A pole at z = 0.99 represented in Q1.15 quantises to either 0x7EB8 = 0.98987… or 0x7EB9 = 0.98993… — the closest representation differs from intended by 60 ppm. For a single pole this is harmless; for a cascade of biquads or a 6th-order IIR the cumulative effect can move closed-loop poles outside the unit circle. Mitigation: factor the controller into 2nd-order sections (biquads), use DF2T form, place poles away from the unit circle when possible, or move to floating-point coefficients.
Integer overflow in the integrator. A 32-bit signed integer accumulator with the integrator term Kp·(Ts/Ti)·e[k] = (a small Q1.31 number) added every cycle can roll over after ~2^31 / increment samples. Always saturate the integrator state explicitly (clamp to bounds before storing), and prefer signed saturating arithmetic (Cortex-M4 has SSAT / USAT single-cycle instructions).
Time-varying sample period. When the “sample” is an encoder edge time rather than a fixed clock — e.g., variable-step velocity estimation in a low-speed motor controller — the discretisation coefficients must change every sample. Naïve approach: recompute coefficients on the fly. Better: keep the controller in fixed T_s and use an observer to fill in samples between encoder edges.
Sampling synchronisation between cascaded loops. A current loop at 10 kHz feeding a speed loop at 1 kHz: the speed loop’s command becomes the current loop’s reference. If the speed loop runs asynchronously from the current loop, there is a sampling-rate-conversion jitter at the boundary — typically harmless if the outer loop is ≥10× slower (Åström’s rule), but problematic at ratios under 5×. Always trigger the inner loop deterministically from the same timer as the outer loop, so the rate ratio is exact and the sample-instant relationship is fixed.
Operator display rates. A 1 kHz current waveform displayed at 10 Hz screen update is aliased: the human sees a slowly-moving pattern that is not the actual waveform. Either decimate-with-filter (proper anti-alias before downsampling) or display peak/RMS/min/max statistics rather than instantaneous samples.
Loop coupling through shared CPU. Two loops at different rates on one MCU share compute time. If the high-rate loop occasionally takes longer (e.g., a state-machine transition), the low-rate loop can miss its deadline — and a missed deadline in feedback control is far worse than added latency, because the discrete model no longer matches reality. Always run worst-case compute analysis; instrument with a GPIO toggle at ISR entry/exit and watch on a scope.
Watchdog in safety-critical loops. Any control loop running on safety-critical hardware (vehicle, medical, aerospace) needs an independent watchdog timer (WDT) that the control task pets each cycle. A WDT timeout fires a non-maskable interrupt or a hardware reset that drops the actuator to a known safe state (motor freewheels, valve closes, fuel cuts). ISO 26262 ASIL-D requires a windowed WDT — must be petted neither too early nor too late, so a stuck loop running at the wrong rate is also detected.
Discrete state-space form
For multi-input multi-output (MIMO) problems, the discrete state-space representation is the workhorse:
x[k+1] = A_d · x[k] + B_d · u[k] y[k] = C_d · x[k] + D_d · u[k]
The discrete matrices come from the continuous (A, B, C, D) via the ZOH equivalent:
A_d = e^{A·T_s}, B_d = (∫₀^{T_s} e^{A·τ} dτ) · B, C_d = C, D_d = D
In MATLAB: sysd = c2d(ss(A,B,C,D), Ts, 'zoh'). In Python: scipy.signal.cont2discrete((A,B,C,D), Ts). The matrix exponential is computed by Padé approximation with scaling-and-squaring; in real-time code the discrete matrices are computed once offline and stored as constants.
The discrete LQR control law u[k] = −K · x[k] with K = (R + B_d^T · P · B_d)^{-1} · B_d^T · P · A_d, where P solves the discrete algebraic Riccati equation, gives optimal regulation in the LQ sense. MATLAB: [K, S, e] = dlqr(Ad, Bd, Q, R). Discrete Kalman filtering uses the same matrices and a dual Riccati. State-space methods are covered in depth in the companion note [[Engineering/state-space-methods]].
Stability checks in z
A discrete polynomial p(z) = z^n + a_{n−1}·z^{n−1} + … + a₀ has all roots inside the unit disc iff it passes the Jury test (the discrete analogue of Routh-Hurwitz). For order 2 the conditions simplify nicely:
| Order 2 polynomial z² + a₁·z + a₀ stable iff |
|---|
| |a₀| < 1 |
| 1 + a₁ + a₀ > 0 |
| 1 − a₁ + a₀ > 0 |
These three inequalities trace the stability triangle in the (a₀, a₁) plane: a triangle with vertices (−1, 0), (1, −2), (1, 2). All discrete 2nd-order systems must fall inside it. Visualising controller tuning in this plane is a useful intuition tool.
Closed-loop pole placement
Given plant G(z) and controller K(z), the closed-loop transfer function is T(z) = G·K / (1 + G·K). The closed-loop poles are roots of 1 + G(z)·K(z) = 0. Discrete root-locus traces these as the controller gain varies. Unlike continuous root-locus which extends to infinity, discrete root-locus is bounded: closed-loop poles must stay inside |z| < 1 for stability, and ideally inside |z| < ρ < 1 for fast settling (typical design target ρ ≈ 0.8…0.95, corresponding to settling time ≈ 4·T_s / (1 − ρ)).
8p. Tools and software
The most-used flows today:
-
Model-based-design teams — Simulink for plant + controller modelling, Control System Toolbox for analysis, Embedded Coder for autocode generation to C for the target MCU. Verification by Software-in-the-Loop (SIL) on the host, then Processor-in-the-Loop (PIL) on the actual silicon, then Hardware-in-the-Loop (HIL) on dSPACE / Speedgoat / NI VeriStand. Standard in automotive, aerospace, defence; common in industrial automation.
-
Embedded-first teams — controller designed in MATLAB or Python offline, coefficients pasted as
const floatarrays into hand-written C using CMSIS-DSP for biquads and matrix math. Faster iteration for small teams; less infrastructure overhead; lower-cost toolchain. Common in startups, consumer products, hobbyist motor drives, smaller industrial-controls vendors. -
FPGA control teams — controller designed in MATLAB/Python, exported to fixed-point Q-format coefficients, hand-coded in SystemVerilog or VHDL with HLS for inner DSP kernels (Vitis HLS). Targets Xilinx Kintex / UltraScale+ / Versal, Intel Stratix / Agilex, Lattice ECP5 / CertusPro-NX. Used for sub-µs latency loops (HFT, high-end multi-axis motion, large-machine drives).
-
PLC control teams — IEC 61131-3 Structured Text (or Function Block Diagram) authored in Studio 5000 / TIA Portal / Codesys. PID function blocks are vendor-supplied with auto-tuning; engineer’s job is tuning parameters and structuring the program for safety and maintainability. Industrial process control, factory automation, water/wastewater, building HVAC.
For mixed-signal / power-electronics specifically: PLECS (Plexim) and PSIM (Powersim) provide circuit-level switching-converter simulation plus controller models, and code-generate to TI C2000 / STM32 / Microchip dsPIC. PLECS Coder is widely used in solar inverter, EV charger, and motor-drive industries.
For HIL: dSPACE SCALEXIO / MicroAutoBox, NI VeriStand + PXI, OPAL-RT HYPERSIM, Speedgoat baseline / performance. Plant model runs in real time at 10–100 kHz on the HIL box; the target ECU runs against simulated I/O over CAN, FlexRay, ARINC 429, analog, digital. This is the only practical way to test corner cases (sensor failures, plant degradation, network outages) without endangering hardware or operators.
For verification and safety: LDRA Testbed, Polyspace (MathWorks), Coverity (Synopsys) cover MISRA C / CERT C / ISO 26262 / DO-178C static-analysis requirements. VectorCAST generates unit and integration tests with MC/DC coverage required by DO-178C DAL-A and ISO 26262 ASIL-D.
11. Cross-references
[[Engineering/classical-control]]— the continuous-time design that becomes the input to discretisation; root-locus, Bode, Nyquist, PID tuning.[[Engineering/state-space-methods]]— discrete LQR, observers, and Kalman filters live naturally in the z-domain; the companion methodology for MIMO digital control.[[Engineering/microcontrollers]]— substrate for control implementation: ADC, PWM, timers, NVIC, DMA, FPU.[[Engineering/digital-logic]]— what FPGA control loops are built from; sample timing comes from synchronous digital design.[[Engineering/op-amps]]— anti-alias filters, current sense, signal conditioning ahead of every ADC.[[Engineering/electric-motors]]— the most common controlled plant; FOC and direct-torque control depend on accurate digital implementation.[[Engineering/power-electronics]]— switching-converter voltage and current loops; synchronous PWM sampling.[[Engineering/mpc-control]](planned) — model predictive control; QP solved every sample.[[Engineering/realtime-embedded]](planned) — RTOS, ISR priority, deadline management.[[Engineering/fpga-design]](planned) — FPGA implementation of high-rate control loops.[[Robotics/motors-electric]](planned) — FOC, SVPWM, sensorless estimation.[[Robotics/bayesian-estimation]](planned) — Kalman / particle filters running at sample rate.[[Languages/Tier3/scientific]](planned) — toolchain for control design and autocode generation.
12. Citations
- Åström, K. J. & Wittenmark, B. (1997). Computer-Controlled Systems: Theory and Design (3rd ed.). Prentice Hall. The canonical reference; sample-rate selection, anti-windup back-calculation derivation, two-degree-of-freedom PID forms.
- Franklin, G. F., Powell, J. D. & Workman, M. L. (1997). Digital Control of Dynamic Systems (3rd ed.). Addison-Wesley. Comprehensive coverage of discretisation methods, ZOH design, sample-rate effects.
- Phillips, C. L., Nagle, H. T. & Chakrabortty, A. (2014). Digital Control System Analysis and Design (4th ed.). Pearson. Practical engineering examples, fixed-point implementation, state-space.
- Ogata, K. (1995). Discrete-Time Control Systems (2nd ed.). Prentice Hall. Classroom standard for Z-transform, state-space discrete control.
- Houpis, C. H. & Lamont, G. B. (1992). Digital Control Systems: Theory, Hardware, Software (2nd ed.). McGraw-Hill. Bridge from theory to hardware implementation.
- Forsythe, W. & Goodall, R. M. (1991). Digital Control: Fundamentals, Theory and Practice. McGraw-Hill. UK-traditional treatment; strong on practical PID and observer implementation.
- Skogestad, S. (2003). “Simple analytic rules for model reduction and PID controller tuning.” Journal of Process Control, 13(4), 291–309. SIMC PID tuning rules used widely in process industries; informs sample-rate-vs-bandwidth selection.
- Cervin, A., Henriksson, D., Lincoln, B., Eker, J. & Årzén, K.-E. (2003). “How does control timing affect performance?” IEEE Control Systems Magazine, 23(3), 16–30. Jitter and missed-deadline effects on closed-loop performance.
- IEEE Std 754-2019. IEEE Standard for Floating-Point Arithmetic. Defines single and double-precision IEEE 754 formats used in float-based controllers.
- IEEE Std 1057-2017. IEEE Standard for Digitizing Waveform Recorders. ADC characterisation methodology, ENOB, SINAD, INL, DNL definitions.
- IEC 61131-3:2013. Programmable Controllers — Part 3: Programming Languages. Defines Structured Text, Function Block Diagram, Ladder Logic for PLCs; standard PID function block.
- MISRA C:2012, with Amendment 1:2016 and Amendment 2:2020. Guidelines for the Use of the C Language in Critical Systems. Rules for safety-critical embedded C; mandatory for ISO 26262 software development.
- ISO 26262:2018 (Parts 1–12). Road Vehicles — Functional Safety. Part 6 covers software development for automotive systems including ASIL-classified control software; mandates windowed WDT, fault detection, lockstep cores, etc. for ASIL-D loops.
- DO-178C / RTCA (2011). Software Considerations in Airborne Systems and Equipment Certification. The aerospace equivalent; DAL-A control software requirements.
- IEC 61508:2010 (Parts 1–7). Functional Safety of Electrical / Electronic / Programmable Electronic Safety-related Systems. Industrial SIL framework underpinning ISO 26262 and DO-178C.
- ARM (2014). Cortex-M4 Technical Reference Manual and CMSIS-DSP documentation. Implementation reference for biquad DF2T, fixed-point Q-format conventions, SSAT/USAT instructions.
- Texas Instruments. TMS320C28x Optimizing C/C++ Compiler User’s Guide and InstaSPIN-FOC application notes. C2000 toolchain and motor-control library reference.
- STMicroelectronics. AN4013 STM32 cross-series timer overview, AN5208 STM32G4 HRTIM, X-CUBE-MCSDK Motor Control SDK. STM32-specific reference for digital-power and motor implementation.