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.
• Course 1 Lesson 48 — 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.

Wiring & bench setup

Same chain as Lab 6.1 minus the timing GPIO: MCP4725 source (tones, and the Part C/D impulse) → ADC A0 = PA0, biquad cascade in firmware, on-chip DAC A2 = PA4 out; scope CH1 = raw input, CH2 = filtered output — CH2 is also where you’ll watch a destabilized impulse response grow and a limit cycle ring.

flowchart LR
  SRC["MCP4725<br/>source tone / impulse<br/>(Lab 3.3 generator)"]
  MCU["NUCLEO-L476RG<br/>ADC A0 = PA0<br/>biquads → DAC PA4"]
  SC1["Siglent CH1<br/>raw input"]
  SC2["Siglent CH2<br/>filtered output"]
  SRC -- "OUT → A0 (PA0)" --> MCU
  MCU -- "A2 (PA4)" --> SC2
  SRC -.-> SC1

flowchart LR
  SRC["MCP4725<br/>source tone / impulse<br/>(Lab 3.3 generator)"]
  MCU["NUCLEO-L476RG<br/>ADC A0 = PA0<br/>biquads → DAC PA4"]
  SC1["Siglent CH1<br/>raw input"]
  SC2["Siglent CH2<br/>filtered output"]
  SRC -- "OUT → A0 (PA0)" --> MCU
  MCU -- "A2 (PA4)" --> SC2
  SRC -.-> SC1

Pin map:

From	To	Pin/jack
Nucleo 3V3 / GND	breadboard + / − rails → MCP4725 VDD / GND	3V3, GND headers
Nucleo D15 / D14	MCP4725 SCL / SDA	PB8 / PB9 (I2C1)
MCP4725 OUT	Nucleo A0 and scope CH1 probe tip	PA0 (ADC1_IN5)
Nucleo A2	scope CH2 probe tip	PA4 (DAC1_OUT1)
CH1 + CH2 ground clips	breadboard − rail	probes at 10× (switch and menu)

Breadboard layout is unchanged from the Lab 6.1 sketch (MCP4725 on the rails) — leave it built. Keep every test signal digitally biased to mid-scale (code 2048 ± amplitude) as in 6.1; an unstable filter will rail PA4 across 0–3.3 V, which is fine into a scope probe and nothing else.

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.

Project & environment setup

Firmware — reuse the Module 6 project firmware/m6-dsp/ (created in Lab 6.1). Nothing new to enable; confirm the .ioc still has:

CubeMX page	Setting
Software Packs	ARM CMSIS-DSP present, ARM_MATH_CM4, hard-float — arm_biquad_cascade_df1_* must link (setup essentials)
ADC1 + DMA + TIM2	IN5 (PA0), circular half-word DMA, TRGO at \(f_s = 16\) kHz — as Lab 6.1
Analog → DAC1	OUT1 (PA4) — the reconstruction output
Connectivity → I2C1 / USART2	400 kHz MCP4725 bus / 115200 8-N-1 VCP

Host — design and pole inspection in the course venv (Toolchain):

source venv/bin/activate   # numpy + scipy (butter/sos2zpk/freqz) + matplotlib + pyserial
mkdir -p labs/lab-6-2/host labs/lab-6-2/captures

One script labs/lab-6-2/host/pole_quantize.py (you write it): the Part A butter → sos2zpk design, the Part C progressive coefficient quantization with the pole-radius plot, and export of the CMSIS coefficient set (signs flipped per Step 3) as a C header.

Keep this lab’s reconciliation in labs/lab-6-2/host/analysis.ipynb — the notebook convention — and export final figures next to it.

Where results go:

Artifact	Path
Bench note	labs/lab-6-2/notes.md
Biquad coefficients (float + Q15, CMSIS sign convention)	firmware/m6-dsp/Core/Inc/biquad_coeffs.h
Pole–zero map + quantized pole-drift plot	labs/lab-6-2/host/pole-drift.png
Scope CSV, stable vs. destabilized impulse response	labs/lab-6-2/captures/impulse-responses.csv
Limit-cycle trace (impulse then zeros, Q15 vs Q31)	labs/lab-6-2/captures/limit-cycle.csv

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 Lesson 48): 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 Lesson 48: the unit circle is the boundary the ROC must contain, and everything you measured is that boundary made physical.

Cross-platform ports & language variants

See the syllabus Implementation tracks for the framing; this is the biquad-specific version. The IIR is firmly sequential: each output feeds back into the next, so it sits in the latency-bound class where the STM32 is home and the GPU is a poor fit.

STM32 bare-metal (C, and Rust). In C the cascade is arm_biquad_cascade_df1_f32 (or _q15, which needs the postShift coefficient scaling). This is where the entire payload of the lab lives: coefficient quantization nudging a high-Q pole across the unit circle, limit cycles, and fixed-point overflow oscillation are all phenomena of the constrained, fixed-point target — you provoke and measure them here, with the DWT cycle counter for the per-sample cost (setup essentials). In Rust (#![no_std]) the recursion is a hand-rolled sosfilt-style loop; Rust’s checked arithmetic surfaces the Q15 overflow-oscillation decision that C leaves silent — you must choose wrapping vs. saturating, which is exactly the right thing to think about here.

Raspberry Pi 5 (Linux userspace, C or NumPy). scipy.signal.sosfilt, or the same sequential C, runs trivially. But float64 makes the whole quantization-and-limit-cycle curriculum vanish — poles don’t move enough to matter, and there is no dead-band. The lab literally cannot be taught on the Pi; that is the point.

Jetson Orin Nano. A deliberate poor fit: the feedback dependency does not parallelize, and parallel-scan reformulations of a recursion rarely pay off for a mere two biquads against kernel-launch overhead. There is nothing for the GPU to do.

Jetson Orin Nano — detailed procedure (embedded Linux)

Short by design — the value of this port is measuring the mismatch, then writing down why. Reuse the Lab 6.1 Jetson harness (edge/ CMake project, .npy test vectors, SciPy reference, p50/p99 timing under jetson_clocks + chrt); board config in the Jetson setup essentials.

1. mkdir -p labs/lab-6-2/edge; export this lab’s test signal and the scipy.signal.sosfilt reference output for your cascade.
2. CPU float32: compile your portable-C biquad cascade (the same shared/ kernel the firmware uses) with -O3 and run it through the harness. It verifies against SciPy and the per-sample cost is tiny — note that NEON barely helps here, because the recursion serializes; compare with how much it helped the FIR.
3. The vanishing curriculum, demonstrated: repeat the lab’s high-Q experiment (the pole radius that destabilized the Q15 build on the STM32) in float32 and float64 on the Jetson. Verify it doesn’t destabilize, and quantify why from the host-side pole analysis you already did in Part A. One paragraph in notes.md: which failure modes of this lab exist on which arithmetic, with the measured evidence.
4. GPU non-run: don’t skip it silently — time a cupyx.scipy.signal.sosfilt call on the same data once, record that the launch+copy overhead exceeds the entire CPU filtering time for any realistic block, and write the number down. “We measured that it loses” teaches more than “everyone says it loses.”

Raspberry Pi 5 differences: identical CPU story (performance governor); no GPU non-run to record.

Measure and compare (fill Measured on each platform):

Platform / build	Per-sample latency	Numerical-robustness content present?	Predicted	Measured
STM32 bare-metal, C (Q15, DWT cycles)	bounded, deterministic	yes — pole drift, limit cycles, overflow	home target	…
STM32 bare-metal, Rust	≈ C	yes — overflow decision made explicit	≈ C	…
Pi 5, sosfilt / float64 C	fast, jittery tail	no — float64 erases it	trivial	…
Jetson GPU	launch-overhead-bound	n/a (won’t parallelize)	net loss	…

The lesson to teach: numerical robustness is a bare-metal / fixed-point discipline, not something float64 Linux userspace ever forces on you.

Same STM32: bare-metal vs RTOS

The runtime axis has a middle rung worth measuring on the MCU itself: run the biquad cascade under FreeRTOS and compare against the bare-metal build. On a filter this cheap the comparison is almost purely about structure, which is the clean way to see it, and it sets up Lab 7.2.

• Bare-metal (above): the DMA half/full-transfer callback runs arm_biquad_cascade_df1_f32 on the ready half inline, in ISR context — a handful of MACs per sample, done before the callback returns.
• FreeRTOS (C): turn the DMA callback into a fast signal and move the cascade into a prioritized DSP task. The HAL_ADC_Conv*CpltCallback only does osSemaphoreRelease(sem) (or queues the ready-buffer index with osMessageQueuePut) and returns; a high-priority task blocks on osSemaphoreAcquire(sem, osWaitForever), runs the biquad cascade, and passes the output to a separate DAC task. Setup: enable FREERTOS → CMSIS_V2 and move the HAL timebase to a spare timer (TIM17) per the setup essentials, then osSemaphoreNew / osThreadNew.
• Rust (RTIC / Embassy): in RTIC, the DMA-complete IRQ is a hardware task that spawns a lower-priority software task carrying the buffer index; the software task runs the recursion. In Embassy, an async task awaits the DMA future and filters on wake. Same ISR→task deferral, statically scheduled.
• What you’ll see: the biquad’s per-block cost is tiny, so the ~few-µs context switch is negligible against the block budget either way — the RTOS overhead effectively vanishes. That means the choice here is purely structural, not a performance question: you move to tasks only to separate the filter cleanly from acquisition/output (the Lab 7.2 pattern), never because the scheduler buys or costs you real time. Measure it with the DWT counter and confirm the deadline margin is essentially unchanged.

Build (same STM32)	Per-block latency/jitter added	Deadline margin	Structural benefit	Measured
Bare-metal, inline in DMA callback	none (filters in ISR)	full	monolithic ISR	…
FreeRTOS, DMA→semaphore→DSP task	+ one context switch (negligible vs budget)	essentially unchanged	filter isolated from acquisition/output	…
Rust RTIC, DMA hw task→sw task	≈ FreeRTOS	essentially unchanged	compile-time-checked task priorities	…

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.