Digital IIR Filter Designer & Biquad Coefficient Synthesizer

1. Fundamental Principles of Digital Infinite Impulse Response (IIR) Filters

Infinite Impulse Response (IIR) filters are recursive discrete-time linear time-invariant (LTI) systems whose impulse response $h[n]$ theoretically extends indefinitely in time. The difference equation governing a causal $N$-th order IIR filter relates the output sequence $y[n]$ to present and past inputs $x[n]$ as well as past feedback outputs:

Taking the bilateral Z-transform yields the rational system transfer function $H(z)$:

Where $z_k$ represent the zeros and $p_k$ represent the poles of the system. For a causal digital filter to exhibit Bounded-Input Bounded-Output (BIBO) stability, all system poles must reside strictly within the interior of the unit circle on the complex Z-plane ($|p_k| < 1.0$ for all $k$).

Compared to Finite Impulse Response (FIR) filters, IIR filters achieve extraordinarily sharp transition band roll-offs with dramatically lower computational complexity—often requiring 5 to 20 times fewer arithmetic operations (multiply-accumulate operations per sample) than an equivalent FIR filter. However, this recursive feedback comes with trade-offs: non-linear phase delay and potential numerical instability under finite-precision coefficient quantization.

2. Classical Approximation Families: Butterworth vs. Chebyshev vs. Elliptic

Analog prototype filter theory provides several classical mathematical polynomials to optimize the magnitude response $|H(jOmega)|$:

• Butterworth Approximation: Maximally flat in both passband and stopband. The magnitude response is monotonically decreasing with zero ripple:
  |H(jOmega)|^2 = frac{1}{1 + (Omega / Omega_c)^{2N}}
  The analog poles lie evenly spaced on a circle of radius $Omega_c$ in the left half of the complex s-plane ($s_k = Omega_c e^{j pi (2k + N - 1)/(2N)}$). Butterworth filters are preferred in high-fidelity audio and biosignal acquisition where passband ripple cannot be tolerated.
• Chebyshev Type I Approximation: Exhibits equal-amplitude ripple in the passband while maintaining monotonic attenuation in the stopband:
  |H(jOmega)|^2 = frac{1}{1 + epsilon^2 T_N^2(Omega / Omega_c)}
  Where $T_N(x)$ is the $N$-th order Chebyshev polynomial and $epsilon = sqrt{10^{R_p/10} - 1}$ is the passband ripple factor. The poles lie on an ellipse in the s-plane, yielding a much steeper cut-off skirt than Butterworth for the same order $N$.
• Chebyshev Type II (Inverse Chebyshev): Maximally flat in the passband with equiripple behavior in the stopband:
  |H(jOmega)|^2 = frac{1}{1 + epsilon^2 [T_N(Omega_r / Omega)]^{-2}}
  Inverse Chebyshev filters introduce transmission zeros directly on the $jOmega$ imaginary axis (mapping onto the Z-plane unit circle), completely nulling specific interference frequencies while keeping the passband flat.

3. The Bilinear Transform and Frequency Pre-Warping

The most robust technique for converting continuous-time analog prototypes $H_a(s)$ into discrete-time digital filters $H(z)$ is the Bilinear Transform (BLT). It is derived from the trapezoidal numerical integration rule, mapping the left half of the s-plane onto the interior of the unit circle in the z-plane via the algebraic substitution:

Because the infinite imaginary axis $s = jOmega$ is mapped onto the finite perimeter of the unit circle $z = e^{jomega}$, the mapping between analog frequency $Omega$ (rad/s) and digital frequency $omega$ (rad/sample) is highly non-linear:

Without compensation, digital filters experience severe frequency warping at frequencies approaching the Nyquist boundary ($f_s / 2$). To guarantee that the digital filter exhibits the desired cutoff frequency $f_c$, the analog cutoff frequency must be pre-warped prior to pole-zero synthesis:

4. Cascaded Second-Order Sections (Biquads) and Numerical Precision

In theoretical textbooks, digital filters are often represented in Direct Form I or Direct Form II as high-order polynomials. In embedded digital signal processors, FPGAs, and audio DSP engines, implementing an IIR filter with order $N ge 4$ as a single direct-form transfer function is universally avoided due to coefficient quantization sensitivity.

The roots of high-order polynomials are hypersensitive to tiny perturbations in polynomial coefficients. In a 6th or 8th order filter, quantizing coefficients to 16-bit or 32-bit floats can displace complex pole pairs across the unit circle, turning a stable design into an explosive oscillator.

To guarantee numerical stability, the overall transfer function $H(z)$ is decomposed into a series of cascaded Second-Order Sections (SOS), commonly called biquad filters:

Each section pairs a complex conjugate pole pair with its nearest zero pair, ensuring that no polynomial degree exceeds 2, confining pole sensitivity strictly to local coefficients.

5. Step-by-Step Worked Numerical Example: 2nd Order Butterworth Lowpass

Consider designing a 2nd-order ($N=2$) Butterworth lowpass filter for a digital audio stream at sampling rate $f_s = 48{,}000text{ Hz}$ with cutoff frequency $f_c = 4{,}800text{ Hz}$:

1. Calculate Normalized Digital Cutoff Frequency: $$omega_d = 2pi frac{4800}{48000} = 0.2pi approx 0.6283185text{ rad/sample}$$
2. Pre-warp the Analog Frequency: $$K = tanleft(frac{omega_d}{2}right) = tan(0.1pi) = tan(18^circ) approx 0.3249197$$
3. Determine Analog Normalized Poles: For $N=2$ Butterworth, normalized poles in the s-plane lie at angles $pm 135^circ$: $$s_{1,2} = -frac{1}{sqrt{2}} pm j frac{1}{sqrt{2}} approx -0.707107 pm j 0.707107$$ The analog prototype denominator is $s^2 + sqrt{2} s + 1$.
4. Apply Bilinear Substitution: Replacing $s leftarrow frac{1}{K} frac{1 - z^{-1}}{1 + z^{-1}}$ into the normalized prototype yields the digital biquad coefficients: $$text{norm} = 1 + sqrt{2} K + K^2 = 1 + 1.4142136(0.32492) + (0.32492)^2 approx 1.565103$$ $$b_0 = frac{K^2}{text{norm}} = frac{0.105573}{1.565103} approx 0.067455$$ $$b_1 = 2 cdot b_0 approx 0.134910, quad b_2 = b_0 approx 0.067455$$ $$a_1 = frac{2(K^2 - 1)}{text{norm}} = frac{2(0.105573 - 1)}{1.565103} approx -1.142981$$ $$a_2 = frac{1 - sqrt{2} K + K^2}{text{norm}} = frac{1 - 0.459505 + 0.105573}{1.565103} approx 0.412802$$
5. Verify Pole Locations: Roots of $1 - 1.142981 z^{-1} + 0.412802 z^{-2} = 0$ give $z = 0.57149 pm j 0.29363$. The pole radius is $r = sqrt{0.57149^2 + 0.29363^2} = sqrt{0.3266 + 0.0862} = sqrt{0.4128} approx 0.6425 < 1.0$. Because $r < 1$, the filter is completely stable.

6. Common DSP Engineering Pitfalls & Best Practices

• Neglecting Frequency Pre-warping: At low ratios ($f_c ll f_s$), warping is negligible ($tan(x) approx x$). But when $f_c > 0.15 f_s$, unwarped designs miss the designated cutoff point by 15% to 50%, altering band-edge transitions dramatically.
• Accumulator Overflow in Direct Form I vs Direct Form II: Direct Form II shares state delay elements, saving RAM, but can suffer from internal node overflow in fixed-point DSPs. Direct Form I is preferred in fixed-point implementations because the summation occurs at the output node where headroom can be managed with extended guard bits.
• Transposed Direct Form II for Audio Processing: Transposed Direct Form II structures exhibit significantly lower round-off noise and prevent zero-input limit cycles when poles are clustered near $z = 1$ (low cutoff frequencies).
• Phase Distortion in Phase-Critical Applications: Unlike symmetric FIR filters, IIR filters have non-linear phase response and non-constant group delay. In audio mastering and ECG acquisition, non-linear phase can smear transients. When linear phase is mandatory with IIR filters, offline zero-phase forward-backward filtering (such as MATLAB/SciPy filtfilt) is used.