On paper a digital filter has infinite precision: its coefficients are real numbers and its arithmetic is exact. On a microcontroller it does not. Coefficients land on a finite grid, every multiply has to be rounded back to the register width, and sums can overflow. These are the finite word-length effects, and they are where a design that looks perfect in double precision can degrade, distort, or even go unstable once it meets the metal.
This page is the theory companion to the practical fixed-point code in the biquad embedded page. It first sets out how fixed-point numbers are represented, then collects the four things that bite, says how big each one is, and shows the structural choices that tame them. The recurring lesson is the same one that motivates cascaded second-order sections: second-order sections are not a stylistic preference, they are what keeps a fixed-point filter alive.
There are really only a handful of distinct effects:
Quantization noise from representing signals and coefficients on a grid.
Coefficient sensitivity: the grid moves the poles and zeros, sometimes catastrophically.
Round-off noise from rounding every intermediate product.
Overflow and the limit cycles that nonlinearity can sustain.
We take them one at a time. The reference throughout is Manolakis and Ingle (Manolakis and Ingle 2011, ch. 15).
Quantization noise and the 6 dB per bit rule
Round a signal to a grid of step \(\Delta\) and you make an error of at most \(\Delta/2\). Model that error as a uniform random variable on \([-\Delta/2, \Delta/2]\) and it has variance
\[\sigma_e^2 = \frac{\Delta^2}{12}.\]
For a \(B\)-bit converter spanning the signal range, \(\Delta\) halves with every extra bit, so the noise power drops by a factor of four. Against a full-scale sinusoid this gives the headline rule of thumb for the signal-to-quantization-noise ratio,
about 6 dB, or one bit, of dynamic range per bit of word length. The demo quantizes a full-scale sine to a range of word lengths and compares the measured SQNR to the formula.
n = np.arange(40000)sine =0.99* np.sin(2* np.pi * n /200) # near-full-scale, no overload from 6 bits upbits = np.arange(6, 17)measured = np.array([sqnr_db(sine, quantize(sine, b)) for b in bits])theory = theoretical_sqnr_db(bits)assert np.max(np.abs(measured - theory)) <1.0# within 1 dB across the sweepfig, ax = plt.subplots(figsize=(8, 4))ax.plot(bits, theory, 'k--', label='6.02B + 1.76 dB')ax.plot(bits, measured, 'o', color='C0', label='measured')ax.set_xlabel('word length B [bits]'); ax.set_ylabel('SQNR [dB]')ax.set_title('Quantization noise: 6 dB per bit')ax.legend(fontsize=9); ax.grid(True, alpha=0.3)fig.tight_layout()plt.show()
Figure 1: Measured signal-to-quantization-noise ratio of a full-scale sinusoid versus word length, against the 6.02B + 1.76 dB rule. The match is within a fraction of a dB.
This is the budget you start from: a 12-bit ADC gives you about 74 dB to work with, a 16-bit one about 98 dB. Everything below eats into it.
Fixed-point representation: Q-formats
When the target has no floating-point unit (a Cortex-M0+, an AVR) or when power and throughput demand the integer path, signals and coefficients are stored as scaled integers in Q-format. A Q\(m\).\(n\) number is a two’s-complement integer read as having \(n\) fractional bits; the two common audio choices both represent a value in \([-1, 1)\):
Format
Bits
Range
Resolution
Use case
Q15
16
\([-1,\ 1 - 2^{-15})\)
\(3.05 \times 10^{-5}\)
Low-power, small filters
Q31
32
\([-1,\ 1 - 2^{-31})\)
\(4.66 \times 10^{-10}\)
High-quality audio
A coefficient stored in Q15 is an integer in \([-32768, 32767]\) standing for \(\text{integer}/32768\), and all arithmetic stays in the integer domain. The programmer’s job is to track where the binary point lands after each operation: a Q15 times Q15 product is Q30 (it has 30 fractional bits in 32), so it must be shifted back before it is stored or summed. Coefficients that exceed \(\pm 1\) (a high-gain peaking EQ, for instance) do not fit Q15 at all and must be pre-scaled, with a compensating gain applied elsewhere. This bookkeeping is the price of dropping the FPU, and it is why the biquad embedded page leans on CMSIS-DSP’s ready-made Q15 and Q31 filters rather than rolling the arithmetic by hand.
Coefficient quantization: the poles move
Quantizing the signal adds noise. Quantizing the coefficients does something more dangerous: it changes the filter. Each coefficient is forced onto the grid,
and those perturbations move the poles and zeros. A sensitivity analysis (Manolakis and Ingle 2011, eq. 15.70) shows that the shift of pole \(p_i\) depends on \(1/\prod_{j\neq i}(p_i - p_j)\): when poles are tightly clustered, as they are in a narrowband or high-order filter, that denominator is tiny and the poles are exquisitely sensitive. Implement such a filter as one high-order direct-form transfer function and quantization can fling a pole clean outside the unit circle.
The cure is structural. Split the filter into second-order sections: each section holds one conjugate pole pair, its two coefficients stay in \([-2, 2]\), and the pole of one section no longer depends on the coefficients of the others. The demo takes a 5th-order elliptic bandpass and quantizes its coefficients two ways.
order, wn = signal.ellipord([0.3, 0.4], [0.2, 0.5], 0.1, 60)b, a = signal.ellip(order, 0.1, 60, wn, btype='band')sos = signal.ellip(order, 0.1, 60, wn, btype='band', output='sos')n_bits =10c_max =2** np.ceil(np.log2(np.max(np.abs(a)) +1.0)) # direct form needs many integer bitsa_q = quantize_coeffs(np.asarray(a), n_bits, c_max=c_max)b_q = quantize_coeffs(np.asarray(b), n_bits, c_max=c_max)sos_q = np.array([quantize_coeffs(s, n_bits, c_max=2.0) for s in sos])# The whole point, checked at render time:assert is_stable(a) andnot is_stable(a_q) and sos_is_stable(sos_q)print(f"direct form needs c_max = {c_max:.0f} (|a| up to {np.max(np.abs(a)):.1f}); "f"quantized direct form stable? {is_stable(a_q)}; quantized SOS stable? {sos_is_stable(sos_q)}")def poles_ba(a):return np.roots(a)def poles_sos(sos):return np.concatenate([np.roots(s[3:6]) for s in np.atleast_2d(sos)])fig, axes = plt.subplots(1, 2, figsize=(10, 5))theta = np.linspace(0, 2* np.pi, 400)for ax, p_inf, p_q, title in [ (axes[0], poles_ba(a), poles_ba(a_q), f'Direct form, {n_bits}-bit'), (axes[1], poles_sos(sos), poles_sos(sos_q), f'Cascade (SOS), {n_bits}-bit')]: ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.7) ax.plot(p_inf.real, p_inf.imag, 'o', mfc='none', mec='C7', label='infinite precision') ax.plot(p_q.real, p_q.imag, 'x', color='C3', label='quantized') ax.set_title(title); ax.set_xlabel('Re'); ax.set_ylabel('Im') ax.set_aspect('equal'); ax.set_xlim(-1.6, 1.6); ax.set_ylim(-1.4, 1.4) ax.axhline(0, color='k', lw=0.3); ax.axvline(0, color='k', lw=0.3) ax.grid(True, alpha=0.3); ax.legend(fontsize=8, loc='upper left')fig.tight_layout()plt.show()
direct form needs c_max = 64 (|a| up to 32.0); quantized direct form stable? False; quantized SOS stable? True
Figure 2: A 5th-order elliptic bandpass under 10-bit coefficient quantization. As one direct-form (b, a) transfer function (left) three of its conjugate pole pairs are pushed outside the unit circle and the filter is unstable. As cascaded second-order sections (right) the same word length leaves every pole in place.
The direct form fails not only because its poles are clustered but because its coefficients have a huge dynamic range (the denominator coefficients of a 10th-degree polynomial run into the tens), so at a fixed word length the grid spacing on them is coarse. The cascade form keeps every coefficient small and every pole independent. This is the numerical argument behind the biquad chapter, made quantitative.
Round-off noise and the accumulator
Even with the poles safely placed, every multiply inside the filter produces a full-width product that has to be rounded back to the register, and each rounding injects a fresh \(\Delta^2/12\) of noise. For a length-\((M{+}1)\) FIR filter that rounds after every tap, the noise powers add and the output noise variance is
\[\sigma_g^2 = (M+1)\,\frac{\Delta^2}{12},\]
so a longer filter is a noisier one. The fix is the double-width accumulator: keep the running sum at full precision and round only once, at the very end. Then a single quantizer contributes and the output noise collapses to
\[\sigma_g^2 = \frac{\Delta^2}{12},\]
independent of filter length. This is exactly why DSP cores carry a wide accumulator: a \(B\)-bit by \(B\)-bit product is about \(2B\) bits, so the running sum is kept at double width plus a few guard bits, and often much wider (CMSIS-DSP’s Q15 biquads accumulate in 64 bits). The same point is made in passing on the biquad embedded page and quantified here.
rng = np.random.default_rng(0)M1 =33# number of tapsh = signal.firwin(M1, 0.25)B =8delta = quantization_step(B)x = (2* rng.random(20000) -1) *0.9# full-ish scale random inputexact = np.convolve(x, h) # full, infinite-precision reference# round after every tap: M+1 quantizers, summed in one aligned frameper_tap = np.zeros(len(exact))for k inrange(M1): contrib = np.zeros(len(exact)) contrib[k:k +len(x)] = h[k] * x per_tap += quantize(contrib, B)single = quantize(exact, B) # wide accumulator: sum exactly, round onceinterior =slice(M1, len(x)) # region where all taps are activevar_per_tap = np.var((per_tap - exact)[interior])var_single = np.var((single - exact)[interior])pred_single = delta**2/12assert0.5< var_single / pred_single <2.0# single quantizer hits the floorassert var_per_tap >8* var_single # per-tap is far noisierprint(f"single-quantizer noise var {var_single:.2e} (Delta^2/12 = {pred_single:.2e})")print(f"per-tap noise var {var_per_tap:.2e}")print(f"wide accumulator is {var_per_tap / var_single:.0f}x quieter for {M1} taps")fig, ax = plt.subplots(figsize=(6, 3.5))ax.bar(['per-tap\n(M+1 quantizers)', 'wide accumulator\n(1 quantizer)'], [var_per_tap, var_single], color=['C3', 'C0'])ax.axhline(pred_single, color='k', linestyle='--', linewidth=0.8, label=r'$\Delta^2/12$ floor')ax.set_yscale('log'); ax.set_ylabel('output noise variance')ax.set_title('Round-off noise: where you round matters')ax.legend(fontsize=8); ax.grid(True, alpha=0.3, axis='y')fig.tight_layout()plt.show()
single-quantizer noise var 5.11e-06 (Delta^2/12 = 5.09e-06)
per-tap noise var 1.03e-04
wide accumulator is 20x quieter for 33 taps
Figure 3: Output round-off noise of a length-33 FIR filter. Rounding once in a wide accumulator matches the Delta^2/12 floor almost exactly; rounding after every tap is roughly an order of magnitude noisier (the (M+1) model is an upper estimate, since the small per-tap products do not fully exercise the quantizer).
Overflow and scaling
Round-off noise sets the floor; overflow sets the ceiling. In two’s-complement fixed point an intermediate sum that exceeds the range wraps around to a large value of the wrong sign, a far more destructive error than rounding because it is unbounded. A sufficient condition to avoid it at a filter output is
You guarantee that by scaling the input down by a factor \(S\). How hard you scale is a real trade-off, because there are three common bounds of increasing conservativeness (left to right):
The rightmost (the \(\ell_1\) norm) never overflows for any bounded input but is so conservative it throws away dynamic range; the leftmost (\(\ell_2\), energy) is the loosest; the peak-gain bound in the middle is the practical choice for narrowband and sinusoidal signals. And here is the tension that has no free lunch: scaling down to dodge overflow lowers the signal against a fixed noise floor, so it costs SNR. Pick the scaling that just barely avoids overflow for the inputs you actually expect, and no more.
When a value does exceed the range despite scaling, the hardware’s overflow behaviour decides how destructive it is. Plain two’s-complement arithmetic wraps: a value one step past the maximum rolls over to the most negative value, tearing the waveform apart. Saturating arithmetic instead clamps to the nearest representable value, which clips the peaks but keeps the signal recognisable. This is why audio fixed-point code saturates rather than wraps (Yiu 2014, ch. 21):
import numpy as npimport matplotlib.pyplot as pltINT16_MIN, INT16_MAX =-32768, 32767n = np.arange(1000)ideal =45000.0* np.sin(2* np.pi *3* n /1000) # exceeds the int16 range# Two's complement wrap: a value past the boundary rolls to the far railwrapped = ((ideal.astype(np.int64) +32768) %65536) -32768# Saturation: clamp to the representable rangesaturated = np.clip(ideal, INT16_MIN, INT16_MAX)# Self-check: wrap injects a near-full-scale jump; saturation never leaves the rangeassert wrapped.max() - wrapped.min() >60000assert saturated.max() <= INT16_MAX and saturated.min() >= INT16_MINfig, axes = plt.subplots(3, 1, figsize=(7, 5), sharex=True)panels = [ (ideal, "Ideal result (exceeds the 16-bit range)"), (wrapped, "Wrapped: two's complement overflow tears the waveform"), (saturated, "Saturated: clamped to the range, mild clipping"),]for ax, (data, title) inzip(axes, panels): ax.plot(n, data, color="C0", lw=1) ax.axhline(INT16_MAX, color="C7", ls="--", lw=0.8) ax.axhline(INT16_MIN, color="C7", ls="--", lw=0.8) ax.set_title(title, fontsize=9) ax.set_ylim(-50000, 50000)axes[-1].set_xlabel("sample n")fig.tight_layout()plt.show()
Figure 4: A signal that exceeds the 16-bit range, handled two ways. Two’s complement wrap (middle) rolls each overflow to the opposite rail, a full-scale discontinuity. Saturation (bottom) clamps to the range, leaving only mild clipping. Dashed lines mark the int16 limits.
No saturating MAC on the Cortex-M4
The Cortex-M4 provides saturating add and subtract (QADD, QSUB) but no saturating multiply-accumulate. A saturating MAC therefore takes two steps, a multiply then a saturating add, costing one extra cycle (Yiu 2014, ch. 21). In tight fixed-point loops the usual practice is the other way around: scale intermediate values down so the accumulator cannot overflow in the first place, and saturate only at the final output.
Limit cycles
The strangest effect appears only because rounding makes the filter nonlinear. A stable IIR filter should decay to zero once its input stops. With a quantizer in the feedback loop it can instead get stuck, circulating a small nonzero value forever: a zero-input limit cycle. It happens when the pole sits close enough to the unit circle that rounding \(a\,y[n-1]\) lands back on \(y[n-1]\), so the decay that should shrink the state never quite does. The dead band is small, on the order of a few quantization steps, but for audio it is an audible idle tone and for a control loop it is a steady-state error that will not go away.
near = first_order_zero_input(a=0.9, y0=0.5, n_bits=6, n_samples=60)far = first_order_zero_input(a=0.5, y0=0.5, n_bits=6, n_samples=60)assertabs(near[-1]) > quantization_step(6) and far[-1] ==0.0fig, ax = plt.subplots(figsize=(9, 3.5))ax.step(np.arange(len(near)), near, where='mid', color='C3', label='a = 0.9 (limit cycle)')ax.step(np.arange(len(far)), far, where='mid', color='C0', label='a = 0.5 (decays to 0)')ax.axhline(0, color='k', lw=0.4)ax.set_xlabel('n'); ax.set_ylabel('y[n]'); ax.set_title('Zero-input limit cycle (6-bit rounding in the loop)')ax.legend(fontsize=9); ax.grid(True, alpha=0.3)fig.tight_layout()plt.show()
Figure 5: Zero-input response of the rounding recursion y[n] = Q{a y[n-1]} from the same nonzero start. With the pole near the unit circle (a = 0.9) the state locks into a nonzero limit cycle; with the pole well inside (a = 0.5) it decays cleanly to zero.
Limit cycles are a structure problem too. Some realizations (error-feedback forms, certain wave and lattice structures) are provably free of them; ordinary direct forms are not. As with coefficient sensitivity, the choice of structure decides whether the effect bites.
Open questions
Predicting limit cycles is genuinely hard. Because the quantizer makes the system nonlinear, there is no general closed-form test for whether a given filter will sustain a limit cycle, or how large its dead band will be. Bounds exist for first- and second-order sections, but for an arbitrary structure the honest answer is still to simulate and listen.
The noise model is a convenient fiction. Treating quantization error as white, uniform, and uncorrelated with the signal is excellent for busy signals and poor for slowly varying or near-constant ones, where the error becomes signal-dependent and tonal. Dither (deliberately adding a small amount of noise before quantizing) decorrelates the error from the signal, trading a modest increase in the noise floor (about 3–5 dB for TPDF) for a clean, distortion-free spectrum. See the dither topic for the full treatment: dither types, the SQNR cost, and when the trade is worth it.
Best structure is filter-specific. Cascade beats direct form for sensitivity, but the ordering and pairing of sections, and which structure minimizes round-off noise, depend on the particular pole-zero constellation. There is no single ordering that is optimal for every filter.
Gold, Bernard, and Charles M. Rader. 1969. Digital Processing of Signals. McGraw-Hill.
Manolakis, Dimitris G., and Vinay K. Ingle. 2011. Applied Digital Signal Processing: Theory and Practice. Cambridge University Press.
Oppenheim, Alan V., and Ronald W. Schafer. 2010. Discrete-Time Signal Processing. 3rd ed. Pearson.
Proakis, John G., and Dimitris G. Manolakis. 2007. Digital Signal Processing: Principles, Algorithms, and Applications. 4th ed. Pearson.
Yiu, Joseph. 2014. The Definitive Guide to ARM Cortex-M3 and Cortex-M4 Processors. 3rd ed. Oxford: Newnes.
Source Code
---title: "Finite Word-Length Effects"subtitle: "What happens to a filter when the bits run out"---On paper a digital filter has infinite precision: its coefficients are real numbers and its arithmetic is exact. On a microcontroller it does not. Coefficients land on a finite grid, every multiply has to be rounded back to the register width, and sums can overflow. These are the **finite word-length effects**, and they are where a design that looks perfect in double precision can degrade, distort, or even go unstable once it meets the metal.This page is the theory companion to the practical fixed-point code in the [biquad embedded page](../../basics/09-biquad/embedded.qmd). It first sets out how fixed-point numbers are represented, then collects the four things that bite, says how big each one is, and shows the structural choices that tame them. The recurring lesson is the same one that motivates [cascaded second-order sections](../../basics/06-filter-design.qmd#second-order-sections): **second-order sections are not a stylistic preference, they are what keeps a fixed-point filter alive.**There are really only a handful of distinct effects:- **Quantization noise** from representing signals and coefficients on a grid.- **Coefficient sensitivity**: the grid moves the poles and zeros, sometimes catastrophically.- **Round-off noise** from rounding every intermediate product.- **Overflow** and the **limit cycles** that nonlinearity can sustain.We take them one at a time. The reference throughout is Manolakis and Ingle [@manolakis2011applied, ch. 15].```{python}#| echo: falseimport numpy as npimport matplotlib.pyplot as pltfrom scipy import signalfrom finite_wordlength import (quantize, quantization_step, sqnr_db, theoretical_sqnr_db, quantize_coeffs, is_stable, sos_is_stable, first_order_zero_input)```<hr>## Quantization noise and the 6 dB per bit ruleRound a signal to a grid of step $\Delta$ and you make an error of at most $\Delta/2$. Model that error as a uniform random variable on $[-\Delta/2, \Delta/2]$ and it has variance$$\sigma_e^2 = \frac{\Delta^2}{12}.$$For a $B$-bit converter spanning the signal range, $\Delta$ halves with every extra bit, so the noise power drops by a factor of four. Against a full-scale sinusoid this gives the headline rule of thumb for the signal-to-quantization-noise ratio,$$\text{SQNR} \approx 6.02\,B + 1.76\ \text{dB},$$about **6 dB, or one bit, of dynamic range per bit of word length**. The demo quantizes a full-scale sine to a range of word lengths and compares the measured SQNR to the formula.```{python}#| label: fig-sqnr#| fig-cap: "Measured signal-to-quantization-noise ratio of a full-scale sinusoid versus word length, against the 6.02B + 1.76 dB rule. The match is within a fraction of a dB."n = np.arange(40000)sine =0.99* np.sin(2* np.pi * n /200) # near-full-scale, no overload from 6 bits upbits = np.arange(6, 17)measured = np.array([sqnr_db(sine, quantize(sine, b)) for b in bits])theory = theoretical_sqnr_db(bits)assert np.max(np.abs(measured - theory)) <1.0# within 1 dB across the sweepfig, ax = plt.subplots(figsize=(8, 4))ax.plot(bits, theory, 'k--', label='6.02B + 1.76 dB')ax.plot(bits, measured, 'o', color='C0', label='measured')ax.set_xlabel('word length B [bits]'); ax.set_ylabel('SQNR [dB]')ax.set_title('Quantization noise: 6 dB per bit')ax.legend(fontsize=9); ax.grid(True, alpha=0.3)fig.tight_layout()plt.show()```This is the budget you start from: a 12-bit ADC gives you about 74 dB to work with, a 16-bit one about 98 dB. Everything below eats into it.<hr>## Fixed-point representation: Q-formatsWhen the target has no floating-point unit (a Cortex-M0+, an AVR) or when power and throughput demand the integer path, signals and coefficients are stored as scaled integers in **Q-format**. A Q$m$.$n$ number is a two's-complement integer read as having $n$ fractional bits; the two common audio choices both represent a value in $[-1, 1)$:| Format | Bits | Range | Resolution | Use case ||--------|------|-------|------------|----------|| Q15 | 16 | $[-1,\ 1 - 2^{-15})$ | $3.05 \times 10^{-5}$ | Low-power, small filters || Q31 | 32 | $[-1,\ 1 - 2^{-31})$ | $4.66 \times 10^{-10}$ | High-quality audio |A coefficient stored in Q15 is an integer in $[-32768, 32767]$ standing for $\text{integer}/32768$, and all arithmetic stays in the integer domain. The programmer's job is to track where the binary point lands after each operation: a Q15 times Q15 product is Q30 (it has 30 fractional bits in 32), so it must be shifted back before it is stored or summed. Coefficients that exceed $\pm 1$ (a high-gain peaking EQ, for instance) do not fit Q15 at all and must be pre-scaled, with a compensating gain applied elsewhere. This bookkeeping is the price of dropping the FPU, and it is why the [biquad embedded page](../../basics/09-biquad/embedded.qmd) leans on CMSIS-DSP's ready-made Q15 and Q31 filters rather than rolling the arithmetic by hand.<hr>## Coefficient quantization: the poles moveQuantizing the *signal* adds noise. Quantizing the *coefficients* does something more dangerous: it changes the filter. Each coefficient is forced onto the grid,$$\hat{a}_k = a_k + \Delta a_k, \qquad \hat{b}_k = b_k + \Delta b_k,$$and those perturbations move the poles and zeros. A sensitivity analysis [@manolakis2011applied, eq. 15.70] shows that the shift of pole $p_i$ depends on $1/\prod_{j\neq i}(p_i - p_j)$: when poles are **tightly clustered**, as they are in a narrowband or high-order filter, that denominator is tiny and the poles are exquisitely sensitive. Implement such a filter as one high-order direct-form transfer function and quantization can fling a pole clean outside the unit circle.The cure is structural. Split the filter into [second-order sections](../../basics/06-filter-design.qmd#second-order-sections): each section holds one conjugate pole pair, its two coefficients stay in $[-2, 2]$, and the pole of one section no longer depends on the coefficients of the others. The demo takes a 5th-order elliptic bandpass and quantizes its coefficients two ways.```{python}#| label: fig-coeff-quant#| fig-cap: "A 5th-order elliptic bandpass under 10-bit coefficient quantization. As one direct-form (b, a) transfer function (left) three of its conjugate pole pairs are pushed outside the unit circle and the filter is unstable. As cascaded second-order sections (right) the same word length leaves every pole in place."order, wn = signal.ellipord([0.3, 0.4], [0.2, 0.5], 0.1, 60)b, a = signal.ellip(order, 0.1, 60, wn, btype='band')sos = signal.ellip(order, 0.1, 60, wn, btype='band', output='sos')n_bits =10c_max =2** np.ceil(np.log2(np.max(np.abs(a)) +1.0)) # direct form needs many integer bitsa_q = quantize_coeffs(np.asarray(a), n_bits, c_max=c_max)b_q = quantize_coeffs(np.asarray(b), n_bits, c_max=c_max)sos_q = np.array([quantize_coeffs(s, n_bits, c_max=2.0) for s in sos])# The whole point, checked at render time:assert is_stable(a) andnot is_stable(a_q) and sos_is_stable(sos_q)print(f"direct form needs c_max = {c_max:.0f} (|a| up to {np.max(np.abs(a)):.1f}); "f"quantized direct form stable? {is_stable(a_q)}; quantized SOS stable? {sos_is_stable(sos_q)}")def poles_ba(a):return np.roots(a)def poles_sos(sos):return np.concatenate([np.roots(s[3:6]) for s in np.atleast_2d(sos)])fig, axes = plt.subplots(1, 2, figsize=(10, 5))theta = np.linspace(0, 2* np.pi, 400)for ax, p_inf, p_q, title in [ (axes[0], poles_ba(a), poles_ba(a_q), f'Direct form, {n_bits}-bit'), (axes[1], poles_sos(sos), poles_sos(sos_q), f'Cascade (SOS), {n_bits}-bit')]: ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.7) ax.plot(p_inf.real, p_inf.imag, 'o', mfc='none', mec='C7', label='infinite precision') ax.plot(p_q.real, p_q.imag, 'x', color='C3', label='quantized') ax.set_title(title); ax.set_xlabel('Re'); ax.set_ylabel('Im') ax.set_aspect('equal'); ax.set_xlim(-1.6, 1.6); ax.set_ylim(-1.4, 1.4) ax.axhline(0, color='k', lw=0.3); ax.axvline(0, color='k', lw=0.3) ax.grid(True, alpha=0.3); ax.legend(fontsize=8, loc='upper left')fig.tight_layout()plt.show()```The direct form fails not only because its poles are clustered but because its coefficients have a huge dynamic range (the denominator coefficients of a 10th-degree polynomial run into the tens), so at a fixed word length the grid spacing on them is coarse. The cascade form keeps every coefficient small and every pole independent. This is the numerical argument behind the [biquad](../../basics/09-biquad/index.qmd) chapter, made quantitative.<hr>## Round-off noise and the accumulatorEven with the poles safely placed, every multiply inside the filter produces a full-width product that has to be rounded back to the register, and each rounding injects a fresh $\Delta^2/12$ of noise. For a length-$(M{+}1)$ FIR filter that rounds after every tap, the noise powers add and the output noise variance is$$\sigma_g^2 = (M+1)\,\frac{\Delta^2}{12},$$so a longer filter is a noisier one. The fix is the **double-width accumulator**: keep the running sum at full precision and round only once, at the very end. Then a single quantizer contributes and the output noise collapses to$$\sigma_g^2 = \frac{\Delta^2}{12},$$independent of filter length. This is exactly why DSP cores carry a **wide accumulator**: a $B$-bit by $B$-bit product is about $2B$ bits, so the running sum is kept at double width plus a few guard bits, and often much wider (CMSIS-DSP's Q15 biquads accumulate in 64 bits). The same point is made in passing on the [biquad embedded page](../../basics/09-biquad/embedded.qmd) and quantified here.```{python}#| label: fig-roundoff#| fig-cap: "Output round-off noise of a length-33 FIR filter. Rounding once in a wide accumulator matches the Delta^2/12 floor almost exactly; rounding after every tap is roughly an order of magnitude noisier (the (M+1) model is an upper estimate, since the small per-tap products do not fully exercise the quantizer)."rng = np.random.default_rng(0)M1 =33# number of tapsh = signal.firwin(M1, 0.25)B =8delta = quantization_step(B)x = (2* rng.random(20000) -1) *0.9# full-ish scale random inputexact = np.convolve(x, h) # full, infinite-precision reference# round after every tap: M+1 quantizers, summed in one aligned frameper_tap = np.zeros(len(exact))for k inrange(M1): contrib = np.zeros(len(exact)) contrib[k:k +len(x)] = h[k] * x per_tap += quantize(contrib, B)single = quantize(exact, B) # wide accumulator: sum exactly, round onceinterior =slice(M1, len(x)) # region where all taps are activevar_per_tap = np.var((per_tap - exact)[interior])var_single = np.var((single - exact)[interior])pred_single = delta**2/12assert0.5< var_single / pred_single <2.0# single quantizer hits the floorassert var_per_tap >8* var_single # per-tap is far noisierprint(f"single-quantizer noise var {var_single:.2e} (Delta^2/12 = {pred_single:.2e})")print(f"per-tap noise var {var_per_tap:.2e}")print(f"wide accumulator is {var_per_tap / var_single:.0f}x quieter for {M1} taps")fig, ax = plt.subplots(figsize=(6, 3.5))ax.bar(['per-tap\n(M+1 quantizers)', 'wide accumulator\n(1 quantizer)'], [var_per_tap, var_single], color=['C3', 'C0'])ax.axhline(pred_single, color='k', linestyle='--', linewidth=0.8, label=r'$\Delta^2/12$ floor')ax.set_yscale('log'); ax.set_ylabel('output noise variance')ax.set_title('Round-off noise: where you round matters')ax.legend(fontsize=8); ax.grid(True, alpha=0.3, axis='y')fig.tight_layout()plt.show()```<hr>## Overflow and scalingRound-off noise sets the floor; **overflow** sets the ceiling. In two's-complement fixed point an intermediate sum that exceeds the range wraps around to a large value of the wrong sign, a far more destructive error than rounding because it is unbounded. A sufficient condition to avoid it at a filter output is$$|y[n]| = \left|\sum_k h[k]\,x[n-k]\right| \le \sum_k |h[k]|\,|x[n-k]| < 1.$$You guarantee that by scaling the input down by a factor $S$. How hard you scale is a real trade-off, because there are three common bounds of increasing conservativeness (left to right):$$\Bigl(\textstyle\sum_k |h[k]|^2\Bigr)^{1/2} \;\le\; \max_\omega |H(e^{j\omega})| \;\le\; \sum_k |h[k]|.$$The rightmost (the $\ell_1$ norm) never overflows for *any* bounded input but is so conservative it throws away dynamic range; the leftmost ($\ell_2$, energy) is the loosest; the peak-gain bound in the middle is the practical choice for narrowband and sinusoidal signals. And here is the tension that has no free lunch: **scaling down to dodge overflow lowers the signal against a fixed noise floor, so it costs SNR.** Pick the scaling that just barely avoids overflow for the inputs you actually expect, and no more.When a value does exceed the range despite scaling, the hardware's overflow behaviour decides how destructive it is. Plain two's-complement arithmetic **wraps**: a value one step past the maximum rolls over to the most negative value, tearing the waveform apart. **Saturating** arithmetic instead clamps to the nearest representable value, which clips the peaks but keeps the signal recognisable. This is why audio fixed-point code saturates rather than wraps [@yiu2014, ch. 21]:```{python}#| label: fig-wrap-vs-saturate#| fig-cap: "A signal that exceeds the 16-bit range, handled two ways. Two's complement wrap (middle) rolls each overflow to the opposite rail, a full-scale discontinuity. Saturation (bottom) clamps to the range, leaving only mild clipping. Dashed lines mark the int16 limits."import numpy as npimport matplotlib.pyplot as pltINT16_MIN, INT16_MAX =-32768, 32767n = np.arange(1000)ideal =45000.0* np.sin(2* np.pi *3* n /1000) # exceeds the int16 range# Two's complement wrap: a value past the boundary rolls to the far railwrapped = ((ideal.astype(np.int64) +32768) %65536) -32768# Saturation: clamp to the representable rangesaturated = np.clip(ideal, INT16_MIN, INT16_MAX)# Self-check: wrap injects a near-full-scale jump; saturation never leaves the rangeassert wrapped.max() - wrapped.min() >60000assert saturated.max() <= INT16_MAX and saturated.min() >= INT16_MINfig, axes = plt.subplots(3, 1, figsize=(7, 5), sharex=True)panels = [ (ideal, "Ideal result (exceeds the 16-bit range)"), (wrapped, "Wrapped: two's complement overflow tears the waveform"), (saturated, "Saturated: clamped to the range, mild clipping"),]for ax, (data, title) inzip(axes, panels): ax.plot(n, data, color="C0", lw=1) ax.axhline(INT16_MAX, color="C7", ls="--", lw=0.8) ax.axhline(INT16_MIN, color="C7", ls="--", lw=0.8) ax.set_title(title, fontsize=9) ax.set_ylim(-50000, 50000)axes[-1].set_xlabel("sample n")fig.tight_layout()plt.show()```::: {.callout-note title="No saturating MAC on the Cortex-M4"}The Cortex-M4 provides saturating add and subtract (`QADD`, `QSUB`) but no saturating multiply-accumulate. A saturating MAC therefore takes two steps, a multiply then a saturating add, costing one extra cycle [@yiu2014, ch. 21]. In tight fixed-point loops the usual practice is the other way around: scale intermediate values down so the accumulator cannot overflow in the first place, and saturate only at the final output.:::<hr>## Limit cyclesThe strangest effect appears only because rounding makes the filter **nonlinear**. A stable IIR filter should decay to zero once its input stops. With a quantizer in the feedback loop it can instead get stuck, circulating a small nonzero value forever: a **zero-input limit cycle**. It happens when the pole sits close enough to the unit circle that rounding $a\,y[n-1]$ lands back on $y[n-1]$, so the decay that should shrink the state never quite does. The dead band is small, on the order of a few quantization steps, but for audio it is an audible idle tone and for a control loop it is a steady-state error that will not go away.```{python}#| label: fig-limit-cycle#| fig-cap: "Zero-input response of the rounding recursion y[n] = Q{a y[n-1]} from the same nonzero start. With the pole near the unit circle (a = 0.9) the state locks into a nonzero limit cycle; with the pole well inside (a = 0.5) it decays cleanly to zero."near = first_order_zero_input(a=0.9, y0=0.5, n_bits=6, n_samples=60)far = first_order_zero_input(a=0.5, y0=0.5, n_bits=6, n_samples=60)assertabs(near[-1]) > quantization_step(6) and far[-1] ==0.0fig, ax = plt.subplots(figsize=(9, 3.5))ax.step(np.arange(len(near)), near, where='mid', color='C3', label='a = 0.9 (limit cycle)')ax.step(np.arange(len(far)), far, where='mid', color='C0', label='a = 0.5 (decays to 0)')ax.axhline(0, color='k', lw=0.4)ax.set_xlabel('n'); ax.set_ylabel('y[n]'); ax.set_title('Zero-input limit cycle (6-bit rounding in the loop)')ax.legend(fontsize=9); ax.grid(True, alpha=0.3)fig.tight_layout()plt.show()```Limit cycles are a structure problem too. Some realizations (error-feedback forms, certain wave and lattice structures) are provably free of them; ordinary direct forms are not. As with coefficient sensitivity, the choice of structure decides whether the effect bites.<hr>## Open questions- **Predicting limit cycles is genuinely hard.** Because the quantizer makes the system nonlinear, there is no general closed-form test for whether a given filter will sustain a limit cycle, or how large its dead band will be. Bounds exist for first- and second-order sections, but for an arbitrary structure the honest answer is still to simulate and listen.- **The noise model is a convenient fiction.** Treating quantization error as white, uniform, and uncorrelated with the signal is excellent for busy signals and poor for slowly varying or near-constant ones, where the error becomes signal-dependent and tonal. Dither (deliberately adding a small amount of noise before quantizing) decorrelates the error from the signal, trading a modest increase in the noise floor (about 3–5 dB for TPDF) for a clean, distortion-free spectrum. See the [dither topic](../dither/index.qmd) for the full treatment: dither types, the SQNR cost, and when the trade is worth it.- **Best structure is filter-specific.** Cascade beats direct form for sensitivity, but the *ordering* and *pairing* of sections, and which structure minimizes round-off noise, depend on the particular pole-zero constellation. There is no single ordering that is optimal for every filter.<hr>## ReferencesThe treatment follows Manolakis and Ingle [@manolakis2011applied, ch. 15]; see also Oppenheim and Schafer [@oppenheim2010discrete] and Proakis and Manolakis [@proakis2007digital]. The coupled-form structure that reduces pole sensitivity near the real axis is due to Gold and Rader [@goldrader1969].::: {#refs}:::
