Lab 6.2 — IIR Biquad Stability

← Course 2 syllabus · Module 6 · Prev: « Lab 6.1 · Next: Lab 6.3 »

Goal

Build an infinite impulse response (IIR) filter as a cascade of second-order sections (biquads) on the STM32, and confront the one thing FIR never made you worry about: stability. You’ll place poles and zeros in the z-plane, watch a pole move as you quantize the coefficients, and see a filter that was stable in double precision go unstable — or ring in a limit cycle — once the coefficients live in 16-bit fixed point. The payoff is judgment: knowing why real DSP firmware ships IIR filters as cascaded biquads in a numerically robust form, and how to check that a design will survive quantization before you flash it. This is the difference between “the math says it’s a filter” and “it’s still a filter on the actual silicon.”

Recommended reading

• Lyons Ch. 6 — IIR filters: the difference equation, poles and zeros, the biquad, and the practical stability and coefficient-quantization pitfalls. Primary reading.
• O-DSP Ch. 6 — structures for discrete-time systems: direct forms, cascade/parallel realizations, and the effects of finite word length (coefficient quantization, limit cycles) — the rigorous treatment.
• Course 1 Week 9 — complex analysis: the z-transform’s region of convergence, poles, and the stability criterion. The theory you are spending here (a causal LTI system is stable iff its ROC includes the unit circle iff all poles are inside it).
• Kuo — fixed-point IIR implementation: Q-format scaling of biquad coefficients, accumulator width, and limit-cycle avoidance.

Equipment & parts

• STM32 Nucleo-64 (NUCLEO-L476RG) with the Module 5 ADC + DMA acquisition project as the front end.
• MCP4725 DAC (or on-chip DAC) to output the filtered signal.
• Signal source into the ADC: DAC-synthesized tones / tone sweep (no built-in generator on the base Siglent).
• Siglent SDS1104X-E scope for magnitude/phase and for catching oscillation or limit-cycle ringing.
• Host with numpy/scipy to design the filter, compute pole locations, and quantize coefficients.

Safety & don’t-break-it

• Same 0–3.3 V analog-pin rule as Lab 6.1 — the input signal must stay inside V_DDA. Nothing new electrically here, but an unstable IIR can rail the DAC output; keep an RC or a scope-only view and don’t drive anything downstream that dislikes a full-scale oscillation.
• An unstable filter is a firmware hazard, not a smoke hazard. When a pole leaves the unit circle the output grows without bound until it saturates — the danger is that you ship it. Always test with worst-case input (an impulse and a full-scale step) and watch for growth or sustained oscillation before trusting a design.
• Overflow oscillation in fixed point can appear even with poles inside the circle: a wrap-around in the accumulator feeds back and self-sustains. Use saturating arithmetic (CMSIS-DSP Q15/Q31 biquads saturate) and design with headroom.

Background

A biquad is a second-order IIR section with transfer function

\[ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}, \]

realized by the difference equation (Direct Form I shown; feedback in \(a\), feedforward in \(b\)):

\[ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]. \]

The zeros are the roots of the numerator, the poles are the roots of the denominator:

\[ p_{1,2} = \frac{-a_1 \pm \sqrt{a_1^2 - 4 a_2}}{2}. \]

For complex-conjugate poles \(p = r e^{\pm j\theta}\) the denominator is \(1 - 2r\cos\theta\,z^{-1} + r^2 z^{-2}\), so the pole radius and angle read directly off the coefficients:

\[ r = \sqrt{a_2}, \qquad \cos\theta = \frac{-a_1}{2\sqrt{a_2}}. \]

Stability criterion. A causal LTI system is stable iff its impulse response is absolutely summable, which for a rational \(H(z)\) means every pole lies strictly inside the unit circle, \(|p_i| < 1\). For a biquad that is exactly \(r = \sqrt{a_2} < 1\), i.e. \(|a_2| < 1\), together with \(|a_1| < 1 + a_2\) (the stability triangle). Equivalently (Course 1 Week 9): the region of convergence of a causal sequence is the exterior of a circle through the outermost pole, and it must include \(|z| = 1\) for the DTFT — the frequency response — to exist. A pole on the circle (\(r = 1\)) is a marginally stable oscillator; outside (\(r > 1\)) the output diverges.

Coefficient quantization moves the poles. The coefficients \(a_1, a_2\) are stored with finite precision. Rounding them to \(B\) bits perturbs the denominator and therefore moves the pole locations. A high-Q pole sits very close to the unit circle (\(r \to 1\)), so a tiny coefficient error can push \(r\) past 1 and destabilize an otherwise-good design. The pole sensitivity is worst for narrowband, high-order filters realized as a single high-order section — which is exactly why you never implement a high-order IIR as one Direct Form section. Cascading second-order biquads keeps each section’s poles individually well-conditioned. Quantitatively, quantizing \(a_2\) by \(\Delta a_2\) moves the radius by roughly \(\Delta r \approx \Delta a_2 / (2r)\), so near \(r=1\) the design margin you need is set by your word length.

Limit cycles. In fixed point, the rounding/truncation at each feedback multiply is a nonlinearity. Even with a stable linear filter and zero input, this can trap the output in a small, self-sustaining oscillation — a limit cycle — that never decays to zero. Its amplitude is bounded by a few LSBs (dead-band), but for an audio or control output it’s an audible/observable defect. Wider accumulators, magnitude truncation, or dithering suppress them; you’ll provoke one deliberately.

Cascade of biquads. A general order-\(2K\) IIR is factored into \(K\) biquads in series:

\[ H(z) = \prod_{k=1}^{K} \frac{b_{0k} + b_{1k} z^{-1} + b_{2k} z^{-2}}{1 + a_{1k} z^{-1} + a_{2k} z^{-2}}. \]

The section ordering and the pole–zero pairing affect the internal dynamic range and thus overflow/noise — a real design choice CMSIS-DSP leaves to you.

Procedure

Part A — Design and inspect poles on the host.

1. Design a low-pass IIR, e.g. a 4th-order Butterworth at \(f_c = 2\) kHz, \(f_s = 16\) kHz, as second-order sections:

   import numpy as np, scipy.signal as sig
   sos = sig.butter(4, 2000, btype="low", fs=16000, output="sos")
   z, p, k = sig.sos2zpk(sos)
   print(np.abs(p))    # pole radii — all must be < 1

2. Plot the pole–zero map (scipy.signal.tf2zpk / a unit-circle plot) and confirm every pole radius \(|p| < 1\). Note the highest-Q section (the pole closest to the circle) — that’s the one quantization will threaten.

3. Note the CMSIS-DSP biquad coefficient sign convention: arm_biquad_cascade_df1_f32 expects \(\{b_0, b_1, b_2, a_1, a_2\}\) per stage with the feedback terms already negated relative to scipy’s sos (which stores \(+a_1, +a_2\)). Convert carefully — a sign error here is the classic “why is my filter unstable” bug.

Part B — Float biquad cascade on the STM32.

4. From the Module 5 DMA project, in the block callback convert ADC codes to centered float and run the cascade:

   #define NSTAGE 2
   arm_biquad_casd_df1_inst_f32 S;
   static float32_t coeffs[5*NSTAGE];  /* {b0,b1,b2,a1,a2} per stage */
   static float32_t state[4*NSTAGE];   /* DF1: 2 x + 2 y per stage */
   arm_biquad_cascade_df1_init_f32(&S, NSTAGE, coeffs, state);
   /* per block: */
   arm_biquad_cascade_df1_f32(&S, xblock, yblock, BLOCK);

5. Sweep tones and verify the magnitude response matches the host freqz(sos). This is your known-good baseline.

Part C — Watch quantization move a pole.

6. On the host, quantize the highest-Q section’s \(a_1, a_2\) to progressively fewer bits (e.g., simulate Q15: np.round(a*32768)/32768, then fewer). Recompute \(r = \sqrt{a_2}\) after each rounding and plot the pole drifting toward (and, for a deliberately high-Q design, past) the unit circle.
7. Flash a version using the quantized coefficients. For a mild Butterworth it stays stable but the response shifts slightly; then design a high-Q peaking/narrowband biquad (pole very near the circle) and repeat — flash it, feed an impulse (a single full-scale sample), and observe on the scope whether the response decays (stable) or grows/sustains (destabilized by rounding).

Part D — Provoke a limit cycle in fixed point.

8. Switch to arm_biquad_cascade_df1_q15 (or _q31) with Q15 coefficients and a high-Q, lightly-damped biquad. Drive it with an impulse, then remove the input (feed zeros) and watch the output on the scope / in a serial log. A limit cycle shows as a small, non-decaying oscillation or a stuck nonzero value where the linear filter predicts decay to 0.
9. Compare Q15 vs. Q31: the wider word pushes the limit-cycle dead-band down and usually kills the visible oscillation.

Deliverable & expected results

Capture: the host pole–zero plot (float vs. quantized), a scope shot of a stable decay vs. a destabilized growing/oscillating impulse response, and the fixed-point limit-cycle trace (impulse then silence). Log the measured pole radius at each quantization level.

Quantity	Predicted	Measured
All pole radii, 4th-order Butterworth \(f_c=2\) kHz	\(|p| < 1\) (largest ≈ 0.7–0.8)	…
Highest-Q pole radius \(r=\sqrt{a_2}\), float	(from design)	…
Same pole radius after Q15 quantization	shifts by \(\approx \Delta a_2/(2r)\)	…
Stability of high-Q design after quantization	may cross \(r>1\) → unstable	…
Limit-cycle amplitude, Q15, zero input	a few LSBs (non-zero)	…
Limit-cycle amplitude, Q31	≈ 0 (below LSB)	…

Analysis & reconciliation

Reconcile the measured pole radius against the coefficient you actually flashed: compute \(r = \sqrt{a_2}\) from the quantized \(a_2\) and confirm it matches where the pole appears to be (stable decay rate, or growth). Tie the observed destabilization back to the sensitivity estimate \(\Delta r \approx \Delta a_2/(2r)\) — a high-Q pole (\(r \approx 0.99\)) with a \(2^{-15}\) coefficient error moves by a tiny amount, but if it was already within that of the circle, it crosses. This is the concrete reason the field uses cascaded biquads: each section’s poles are individually far from the circle relative to the whole high-order polynomial’s roots, so the same word length buys far more margin. For the limit cycle, confirm the residual oscillation amplitude scales with the LSB (halves roughly per extra bit) and vanishes with the wider accumulator — evidence it’s a finite-word-length effect, not a design instability. Cross-check the whole thing against Course 1 Week 9: the unit circle is the boundary the ROC must contain, and everything you measured is that boundary made physical.

Going further

• Implement the same filter as a single 4th-order Direct Form II section and quantize it — watch it destabilize at a word length where the biquad cascade is still fine. This is the whole argument for biquads, shown by experiment.
• Try Direct Form II transposed (arm_biquad_cascade_df2T_f32) and compare its internal-noise/overflow behavior to DF1 — the transposed form has different round-off noise properties.
• Reorder the cascade sections and pair poles with zeros differently; measure the change in output noise floor (ties to Lab 6.4).
• Design a notch (a pair of zeros on the unit circle plus nearby poles) to kill a 60 Hz interferer and confirm the notch depth vs. how close you dare place the poles.