← Back to Autonomy
Interactive Reference · Embedded Control Systems
Vol. 01   /   Digital Domain

Digital Control,
from s to z, in practice.

A working reference for engineers who ship control code to a microcontroller. Every plot below is generated live from the math — adjust the knobs and watch the difference equations behave. Terminology, practical implementation notes, and the rules of thumb you keep reaching for are collected alongside.

By Majid Mazouchi

§ 01

Sampling & aliasing — the first decision

Everything downstream inherits the consequences of your sampling rate. Choose it wrong and no amount of controller cleverness recovers the information you discarded. Choose it right and the continuous-time intuition you trained on in graduate school mostly still works.

A continuous signal \(x(t)\) sampled at period \(T_s\) becomes \(x[k] = x(kT_s)\). The Nyquist–Shannon theorem says you must sample at strictly more than twice the highest frequency present in \(x(t)\) to reconstruct it. Every sinusoid above \(f_s/2\) folds down to an alias frequency \(|f - n f_s|\) and is indistinguishable from a real signal at that alias. The analog anti-aliasing filter is therefore not optional — it is how you stop reality from lying to your DSP.

Terminology

Sample period \(T_s\)
Time between samples, in seconds. Its reciprocal \(f_s = 1/T_s\) is the sample rate in hertz.
Nyquist frequency
\(f_N = f_s/2\). The hard ceiling on representable signal content.
Aliasing
A signal component at frequency \(f > f_N\) appearing in the sampled data at \(|f - n f_s|\), masquerading as a lower-frequency signal.
Anti-aliasing filter
Analog low-pass filter placed before the ADC that attenuates content above \(f_N\). No digital filter can undo aliasing.
Zero-order hold (ZOH)
DAC behavior: the output stays at the last computed value until the next update. This introduces an average delay of \(T_s/2\).
Fig. 1  ·  Aliasing — watch a signal fold
interactive
  • Closed-loop control: \(f_s \geq 10\)–\(20 \times\) closed-loop bandwidth. For a 1 kHz current loop, 10–20 kHz sampling is the working range.
  • Open-loop monitoring / logging: \(f_s \geq 2.5 \times\) the highest component of interest, with a matching analog filter.
  • Anti-alias filter corner: place at \(\sim\)0.3–0.4 \(f_s\) for a simple first-order RC; closer to \(f_N\) if you use a steep (4th+ order) active filter.
  • PMSM current sensing: synchronous sampling at the PWM carrier (or its midpoint) naturally rejects switching ripple without an aggressive analog filter.
  • Suspicious low-frequency wobble that wasn't there on the scope is almost always an alias. Start by doubling \(f_s\) and see if it moves.
§ 02

The z-transform, as a working tool

For discrete signals the z-transform plays the role Laplace plays for continuous signals. Defined as \(X(z) = \sum_{k=0}^{\infty} x[k] z^{-k}\), it turns difference equations into algebra and delays into multiplications by \(z^{-1}\). You only need to use a handful of facts — the derivations are in every textbook.

The working facts

\(z^{-1}\)
A one-sample delay. The expression \(y[k] = 0.9\,y[k-1] + 0.1\,u[k]\) is \(Y(z) = \dfrac{0.1}{1 - 0.9 z^{-1}} U(z)\).
s-plane to z-plane
Under exact mapping \(z = e^{sT_s}\), the LHP (stable continuous) maps to the inside of the unit circle, the imaginary axis maps to the unit circle, and the RHP maps outside.
Pole magnitude
\(|p| < 1 \Rightarrow\) stable mode, \(|p| = 1 \Rightarrow\) marginally stable, \(|p| > 1 \Rightarrow\) unstable.
Pole angle
The argument of a complex-conjugate pole pair sets the oscillation frequency: \(\omega_d = \arg(p)/T_s\) in rad/s.
Pole radius
Controls decay. A pole at radius \(r\) has time constant \(\tau \approx -T_s/\ln(r)\), so smaller radius = faster decay.

That last pair — pole radius sets decay, pole angle sets oscillation — is the intuition you need when you see a z-domain pole plot. The closer a pole sits to \(z = 1\) (i.e. radius near 1, angle near 0), the slower and more integrator-like the mode. Poles near \(z = -1\) oscillate at the Nyquist rate and are usually a warning sign.

Continuous intuition says “poles in the left half-plane are good.” Discrete intuition says “poles inside the unit circle are good.” Same physics, different picture — and the unit-circle picture is the one you debug against. — Working rule
§ 03

Discretization — four ways to turn \(G(s)\) into \(G(z)\)

You have a continuous design — a notch filter, a lead compensator, a plant model for simulation — and you need a difference equation that runs at \(T_s\). There are four methods in common industrial use. Each corresponds to a specific substitution \(s \mapsto f(z)\) or a specific integration rule.

MethodSubstitutionWhat it does wellWhat it does badly
Forward Euler s = (z − 1)/T Trivial to implement; explicit update. Can map stable continuous poles outside the unit circle for fast plants or large \(T_s\). Unsafe for stiff systems.
Backward Euler s = (z − 1)/(zT) Always preserves stability; implicit but closed-form for linear. Adds artificial damping; shifts poles toward \(z = 0\). Distorts frequency response noticeably.
Tustin (bilinear) s = (2/T)·(z − 1)/(z + 1) Preserves stability. Best frequency-response match at low \(\omega\). The default for filter translation. Frequency warping: \(\omega_d = (2/T)\tan(\omega_a T/2)\). Pre-warp around the critical frequency when it matters.
Zero-order hold Gd(z) = (1 − z−1)·Z{G(s)/s} Physically exact for a plant driven by a ZOH DAC — which is what you actually have. More work to compute symbolically. Not appropriate for discretizing a filter or compensator.

The distinction that matters in practice: discretize the plant with ZOH (that is what the DAC physically does), and discretize the compensator with Tustin (that preserves the frequency response you designed for). Forward and backward Euler are for situations where \(T_s\) is already much smaller than every time constant and you want the cheapest possible code.

Fig. 2  ·  First-order plant, four methods, one step
interactive
small = safe, large = revealing
Watch this: push \(T_s/\tau\) above \(\sim\)1 and Forward Euler goes unstable — its pole leaves the unit circle — while the other three stay well-behaved. This is the canonical failure mode of explicit Euler on a fast plant.
  • Compensators and filters: Tustin. Pre-warp around the corner frequency if that corner sits above \(\sim\)0.1 \(f_s\).
  • Plant models for controller synthesis: ZOH. Matches the physical DAC behavior.
  • Fast, simple digital integrators inside a control loop: Backward Euler. Unconditionally stable and adds a tiny phase lag you can absorb.
  • Never use Forward Euler on a plant whose fastest time constant is comparable to \(T_s\). Your simulation will lie, then diverge.
  • Pre-warping frequency \(\omega_p\): replace \(2/T\) with \(\omega_p/\tan(\omega_p T/2)\) in the Tustin formula to keep gain and phase exact at \(\omega_p\).
§ 04

Digital PID — the form you actually deploy

Every control textbook prints the ideal PID equation and moves on. The version that runs on an actual microcontroller looks rather different. It has a filtered derivative, anti-windup on the integrator, proper handling of the first sample, and it thinks carefully about whether to update in position form or velocity form.

Position form vs. velocity form

The position form computes the full actuator command each tick:

positionu[k] = Kp·e[k] + Ki·∑e[i]·Ts + Kd·(ef[k] − ef[k-1])/Ts

The velocity form updates the command incrementally:

velocityΔu[k] = Kp·(e[k] − e[k-1]) + Ki·Ts·e[k] + Kd·(ef[k] − 2ef[k-1] + ef[k-2])/Ts
u[k] = u[k-1] + Δu[k]

Velocity form has one beautiful property for industrial practice: integrator windup is structural. If the actuator saturates and you simply don't add \(\Delta u\) to \(u\), the integrator effectively stops accumulating. Position form requires explicit anti-windup logic. On the other hand, velocity form's command drifts on any computational upset — a missed tick leaves \(u\) frozen at the wrong value. Position form fails gracefully on restart because it recomputes from scratch.

PID terminology — beyond the textbook

Derivative filter
The raw derivative \((e[k]-e[k-1])/T_s\) is a high-pass with unbounded gain at Nyquist. Real implementations filter the measurement (not the error) through a first-order low-pass with corner \(\omega_d = 1/T_f\), giving the filtered derivative \(D(s) = K_d s/(T_f s + 1)\). Choose \(T_f \approx K_d/(N K_p)\) with \(N = 8\text{–}20\).
Setpoint weighting
Feed the setpoint into the P and D terms at reduced weight (\(b, c < 1\)) to suppress the derivative kick from step setpoint changes without changing disturbance rejection.
Anti-windup, back-calculation
When \(u_{\rm cmd}\) saturates to \(u_{\rm sat}\), add \(K_t (u_{\rm sat} - u_{\rm cmd})\) to the integrator state. \(K_t \approx \sqrt{K_i / K_d}\) is a common tuning.
Anti-windup, conditional integration
Simply freeze the integrator whenever the actuator is saturated. Simpler to reason about; slightly slower to recover.
Bumpless transfer
When switching control modes or toggling auto/manual, initialize the integrator state so \(u[k]\) matches the current actuator value. Avoids the characteristic thump on mode change.
Derivative-on-measurement
Compute D from \(y\), not from \(e\). Eliminates the derivative spike on a setpoint step (a.k.a. "derivative kick") without any setpoint weighting.
Fig. 3  ·  Digital PID on a 2nd-order plant
interactive · closed-loop

A few things worth trying in the demo: push \(K_i\) hard with a tight actuator limit and watch the overshoot; untick anti-windup and watch it get dramatically worse; swap derivative-on-measurement off and watch the control signal spike on every setpoint change. These are the classroom pathologies that show up in real systems.

  • Never ship a PID without anti-windup. The cost is a few lines of code; the cost of omitting it is a customer incident.
  • Always filter the derivative. Unfiltered D amplifies measurement noise by \(1/T_s\) — on a 10 kHz loop that is \(10{,}000\times\) gain on the noise floor.
  • Derivative on measurement unless you have a specific reason not to. Setpoint steps should not ring the actuator.
  • Integrator reset on mode changes, enable/disable, and whenever the loop is effectively open (e.g., actuator fault). Stale integrator state is the most common cause of "why did it kick so hard on startup?"
  • Parallel form is easier to tune intuitively; standard form (\(K_p(1 + 1/(T_i s) + T_d s)\)) is easier to compare to textbook methods. Pick one and stick with it in your codebase.
  • Saturation awareness: always know whether \(u[k]\) hit the limit. Log it. A loop that saturates 20% of the time is telling you something — usually that \(K_i\) is too high or the actuator is undersized.
§ 05

Stability, read from the unit circle

In continuous time you check stability by looking at whether poles sit in the left half-plane. In discrete time you check whether poles sit inside the unit circle. The mapping is not just topological — the character of the response is legible directly from the pole location. Radius controls decay, angle controls oscillation frequency, and how close you get to \(z = 1\) (DC) versus \(z = -1\) (Nyquist) tells you what kind of mode you're looking at.

Fig. 4  ·  Pole placer — drag the pole, read the response
interactive · drag the red dot

Reading the unit circle

Pole near \(z = 1\)
Slow, integrator-like, positive step response. A pole at exactly \(z = 1\) is a pure integrator (marginally stable).
Pole near \(z = -1\)
Alternating-sign mode at the Nyquist rate. If unintended, this usually means your sample rate is too low for the dynamics.
Complex pair at radius \(r\), angle \(\theta\)
Damped oscillation: decay per sample = \(r\), oscillation frequency = \(\theta/T_s\) rad/s. Equivalent continuous \(\zeta\) and \(\omega_n\) can be read via \(z = e^{sT_s}\).
Pole on the unit circle
Sustained oscillation — not stable, not unstable. Comes up in oscillators and resonant filters; every other use is probably a bug.
Pole outside the unit circle
Unstable. Each sample the mode multiplies by \(r > 1\). If your code produces such a pole, something is wrong — a sign error in a filter coefficient, an Euler discretization with too large \(T_s\), or a positive feedback loop.
§ 06

Quantization — when the number line has gaps

Finite word length is where the ideal math and the silicon disagree. Every ADC reading, every multiplication, every state variable lives on a grid with spacing LSB. If you do not design around it, you get limit cycles, biased estimates, and silent loss of dynamic range.

Terminology

LSB
Least significant bit — the resolution of your representation. For an \(N\)-bit signed integer on range \([-V, V]\), LSB = \(2V/2^N\).
Q-format
Fixed-point convention: Qm.n means \(m\) integer bits and \(n\) fraction bits. Total word length is \(m + n + 1\) for signed. Range \(\approx [-2^m, 2^m)\), resolution \(2^{-n}\).
Quantization SNR
For a full-scale sinusoid, \(\mathrm{SNR} \approx 6.02 N + 1.76\) dB. Each additional bit is 6 dB.
Dither
Small pseudo-random signal added before quantization to decorrelate quantization error from the signal. Turns a structured limit cycle into white noise.
Overflow vs. saturation
On overflow, uncontrolled wraparound; on saturation, the value is clipped to the representable range. Always use saturating arithmetic for signals; never let a motor controller wrap.
Limit cycle
A small-amplitude oscillation that persists because the quantizer turns a decaying mode into a non-decaying one. Common in narrow-band digital filters with insufficient word length.
Fig. 5  ·  Quantization noise vs. word length
interactive
  • State variables need more bits than signals. A 12-bit ADC feeding an IIR filter typically needs 24+ bit internal state to avoid limit cycles.
  • Multiplications grow word length. Two Qm.n values multiply to Q(2m).(2n). Decide immediately whether you truncate, round, or widen — and be consistent.
  • Scale for headroom. A signal nominally ±1 should be stored at around 50–70% of full scale, not 99%, so transients don't saturate. Budget 3–6 dB of headroom.
  • Round, don't truncate, unless you are certain the bias is acceptable. Truncation biases toward −∞; rounding is unbiased.
  • Integrators are the worst offenders. Small quantization errors accumulate. Either use extra bits on integrator state, or use a leaky integrator, or reset periodically.
  • Float-to-fixed migration: instrument the float simulation to log the min/max of every signal over a representative workload, then pick Q-formats with 2× headroom. Never guess.
§ 07

Computational delay — the sample you always lose

On a real MCU the ADC reading arrives at time \(kT_s\), the control code runs for some fraction of a sample period, and the DAC update is applied at or before \((k+1)T_s\). In the idealized textbook picture the update is instantaneous. In practice you incur a delay between \(0\) and \(T_s\), and if the update is deferred to the next PWM boundary it is almost always very close to exactly one sample. That delay is a frequency-dependent phase lag that eats your phase margin.

lag at loop bandwidth∠(e−jωTs) = −ωTs rad = −360·f·Ts deg

So at a loop bandwidth equal to \(f_s/10\), one sample of computational delay costs \(36°\) of phase margin. That is enormous. Three options to manage it:

  1. Absorb it: design the controller on the model that already includes the one-sample delay. Your continuous-time phase-margin target must be met by the discrete model with \(z^{-1}\) in the loop, not the naive continuous plant.
  2. Predict around it: use a Smith predictor or a state observer that predicts \(y[k+1]\) from \(y[k]\) and \(u[k]\), and compute the control based on the prediction. Works well when the plant model is accurate.
  3. Remove it: use double-update PWM (update at midpoint and edge) or split the computation across two timers so the update latency is \(\ll T_s\). This is how high-end motor control ASICs achieve 100 kHz+ current loops with low apparent delay.
Fig. 6  ·  The phase margin eaten by one sample of delay
interactive
  • Assume one sample of delay by default in your control synthesis. Be pleasantly surprised if it turns out to be less.
  • Keep \(f_{\rm bw} \leq f_s/10\) and a one-sample delay costs you only ≈36° of phase margin — manageable.
  • Below \(f_s/5\) bandwidth, delay compensation becomes essential. Above \(f_s/3\), you are in prediction-observer territory or you are not really closing the loop you think you are.
  • Measure, don't assume: toggle a GPIO at the top and bottom of your ISR. The scope will tell you the actual delay in nanoseconds, which is what goes into the model.
  • Jitter is often worse than constant delay. A 50% jitter on a 100 μs ISR is a 50 μs uncertainty — treat it as additional phase noise in your stability budget.
§ 08

The cheat sheet

Everything above, condensed. This is the page you screenshot and tape above your desk.

Sample rate

fs ≥ 10–20× closed-loop bandwidth. Anti-alias filter corner at 0.3–0.4 fs.

Stability region

Inside the unit circle: |z| < 1. Radius sets decay, angle sets oscillation frequency ω = arg(z)/Ts.

Discretize the plant

ZOH. That is what the DAC physically does. Use for controller synthesis and simulation.

Discretize the compensator

Tustin. Pre-warp if the critical frequency is above 0.1 fs.

Never use Forward Euler

…on a plant with time constants near Ts. It can map stable continuous poles outside the unit circle.

Filtered derivative

Low-pass corner N·Kp/Kd with N = 8–20. Compute D from measurement, not error.

Anti-windup

Back-calculation gain Kt ≈ √(Ki/Kd). Always. No exceptions.

Integrator reset

On mode change, enable, disable, fault recovery, bumpless transfer. Stale state = startup kick.

Quantization SNR

6.02·N + 1.76 dB for a full-scale sinusoid. Each bit = 6 dB.

Word length budget

Signals: ADC width + 2 bits. State: 2× signal width. Accumulators: whatever prevents saturation over worst-case.

Computational delay

Assume 1 sample by default. Each sample of delay = 360·f·Ts° phase at frequency f.

Bandwidth ceiling

Practical fbw ≤ fs/10. Above fs/5 you need delay compensation or a Smith predictor.

Tustin formula

s → (2/Ts)·(z−1)/(z+1). Pre-warp: replace 2/Ts with ω/tan(ωTs/2).

ZOH of 1st-order

For G(s) = 1/(τs+1): G(z) = (1−a)/(z−a) with a = exp(−Ts/τ).

Always instrument

Log: saturation flags, integrator state, control effort, loop execution time. The bug you don't measure is the bug you don't fix.

Float-to-fixed

Instrument in float. Log per-signal min/max. Pick Q-format with 2× headroom. Test on worst-case stimulus.