Infinite impulse response filters are staples of real-time biosignal processing, but direct-form implementations suffer catastrophic numerical instability on mobile hardware. A fourth-order Butterworth bandpass filter represented as a single transfer function will accumulate floating-point rounding errors that destroy passband flatness and stopband attenuation. The solution: decompose high-order filters into cascaded second-order sections—biquads—each with bounded coefficient magnitudes and predictable pole placement.

Photoplethysmography signals from smartphone cameras or wearable PPG sensors demand aggressive filtering. Heart rate sits between 0.5 and 4 Hz, but ambient light flicker at 50 or 60 Hz, motion artifacts, and DC drift dominate the raw signal. A typical processing chain applies a DC-blocking highpass at 0.4 Hz, a lowpass at 5 Hz to reject noise, and a notch filter at the mains frequency. Implementing this as a single eighth-order IIR filter is mathematically elegant but numerically disastrous on ARM Cortex processors with 32-bit floats.

Why Direct Form Fails

Consider a fourth-order Butterworth lowpass with cutoff at 5 Hz and sampling rate 100 Hz. The transfer function has four poles near the unit circle in the z-plane. When you expand the denominator polynomial, coefficients range from 1.0 down to 10-6. Multiplying a signal sample by 10-6 in single precision loses 20 bits of information. Feedback from previous outputs magnifies these errors exponentially. Within seconds, the filter state diverges or oscillates at spurious frequencies.

The root cause is coefficient sensitivity. High-order polynomials are ill-conditioned: tiny changes in coefficients produce large changes in pole locations. A pole that should sit at radius 0.95 and angle 12° might drift to 0.98 and 15° after quantization, shifting the frequency response by hundreds of millihertz—unacceptable for heart rate variability analysis.

Second-Order Section Decomposition

Every high-order IIR filter can be factored into a product of second-order transfer functions. A fourth-order filter becomes two biquads in series; an eighth-order filter becomes four. Each biquad has the canonical form:

H(z) = (b0 + b1·z⁻¹ + b2·z⁻²) / (1 + a1·z⁻¹ + a2·z⁻²)

This structure implements a complex conjugate pole pair and a complex conjugate zero pair. Coefficients remain order-unity: |a1|, |a2|, |b0|, |b1|, |b2| are typically between 0.1 and 2.0. Quantization errors stay localized to each section and do not compound across the cascade.

Pole pairing matters. When decomposing an eighth-order filter into four biquads, you have 24 possible orderings. The optimal strategy pairs poles with similar radii and places sections with the highest Q (sharpest resonance) last in the chain. This minimizes intermediate signal magnitude, preventing overflow in fixed-point implementations and reducing quantization noise in floating-point.

Implementation: Transposed Direct Form II

The transposed direct form II is the standard biquad topology for mobile DSP. It uses two state variables per section and has superior numerical properties compared to direct form I. Here's a production-ready Swift implementation for 32-bit float processing:

struct Biquad {
  var b0, b1, b2, a1, a2: Float
  var s1: Float = 0  // state
  var s2: Float = 0

  mutating func process(_ x: Float) -> Float {
    let y = b0 * x + s1
    s1 = b1 * x - a1 * y + s2
    s2 = b2 * x - a2 * y
    return y
  }
}

For a PPG processing pipeline at 100 Hz sampling rate, a typical cascade includes:

  • One biquad highpass at 0.4 Hz (second-order Butterworth, removes DC drift)
  • Two biquads forming a fourth-order lowpass at 5 Hz (flattens passband, -80 dB at 20 Hz)
  • One biquad notch at 50 Hz with Q=30 (rejects mains interference, -40 dB depth)

The coefficient calculation uses bilinear transform with frequency prewarping. For a second-order Butterworth lowpass at normalized frequency Ω = 2πf/fs:

let K = tan(Ω / 2)
let norm = 1 / (1 + sqrt(2) * K + K * K)
b0 = K * K * norm
b1 = 2 * b0
b2 = b0
a1 = 2 * (K * K - 1) * norm
a2 = (1 - sqrt(2) * K + K * K) * norm

Notch filters require separate design equations. For a notch at frequency f_notch with bandwidth BW:

let Ω = 2 * .pi * f_notch / fs
let BW_rad = 2 * .pi * BW / fs
let β = tan(BW_rad / 2)
let γ = cos(Ω)
b0 = 1 / (1 + β)
b1 = -2 * γ * b0
b2 = b0
a1 = b1
a2 = (1 - β) * b0

Cascade Ordering and Gain Scaling

Ordering biquads affects both numerical stability and computational efficiency. The general rule: place lowpass and highpass sections first to attenuate out-of-band energy early, then apply notch filters, then any peaking or resonant sections. For the PPG chain described above, the optimal order is highpass → lowpass → notch.

Gain scaling prevents intermediate overflow. Each biquad can amplify its input by a factor equal to its peak magnitude response. For a fourth-order lowpass formed from two biquads, the first section might have 6 dB gain at DC. If the input signal has a strong DC component, the intermediate output could exceed ±1.0 even if the final output is well-behaved. Distribute the overall cascade gain across sections by scaling each biquad's b-coefficients:

let cascadeGain: Float = 0.25  // -12 dB
let perSectionGain = pow(cascadeGain, 1.0 / Float(numSections))
biquad.b0 *= perSectionGain
biquad.b1 *= perSectionGain
biquad.b2 *= perSectionGain

Frequency Response Validation

After coefficient calculation, sweep the frequency response to verify passband ripple, stopband attenuation, and group delay. Generate 512 samples of a chirp signal sweeping 0 to 50 Hz, process through the cascade, and compute the FFT magnitude. For a PPG bandpass, expect:

  • Passband 0.5–4 Hz: ±0.5 dB ripple
  • Stopband 6 Hz: