IIR and FIR filters

Deep dive

How it works

In: sample stream x[n], coefficients h (FIR) or b, a (IIR) · Out: filtered stream y[n] · State: last N-1 inputs (FIR) or 2 values per biquad section (IIR) · Cost: N multiply-adds per sample (FIR), about 5 per biquad (IIR)

1. Pick a target: type (low/high/band-pass/stop), cut-off fc/fs, ripple, stopband attenuation, allowed delay.
2. FIR, window method: take the ideal brick-wall impulse response 2 fc sinc(2 fc (n-M)), M=(N-1)/2, multiply by a window (Hamming, Blackman, Kaiser) to cut the infinite tails, normalise DC gain to 1. Equiripple (Parks-McClellan) gives the fewest taps for a given spec.
3. IIR, bilinear method: take an analog prototype (Butterworth, Chebyshev, elliptic), map s -> 2 (z-1)/(z+1) after pre-warping the cut-off with tan(pi fc/fs). Result: biquad sections (see biquad filters).
4. Run it sample by sample; keep N-1 inputs (FIR) or 2 state values per section (IIR, transposed direct form II).

FIR:  x -> [z^-1]-[z^-1]-[z^-1]...   y = h0 x0 + h1 x1 + h2 x2 + ...
IIR:  x -> (+) -> b-taps ----> y      poles from feedback: y also re-enters via a-taps

Code

/// Linear-phase lowpass FIR by windowed sinc. `fc` = cutoff / fs (0..0.5), `n` odd.
pub fn fir_lowpass(fc: f64, n: usize) -> Vec<f64> {
    assert!(n % 2 == 1);
    let m = (n - 1) as f64 / 2.0;
    let mut h: Vec<f64> = (0..n)
        .map(|i| {
            let x = i as f64 - m;
            let sinc = if x == 0.0 { 2.0 * fc } else { (2.0 * PI * fc * x).sin() / (PI * x) };
            let w = 0.54 - 0.46 * (2.0 * PI * i as f64 / (n - 1) as f64).cos(); // Hamming
            sinc * w
        })
        .collect();
    let sum: f64 = h.iter().sum();
    h.iter_mut().for_each(|v| *v /= sum); // exact unity gain at DC
    h
}

/// Streaming FIR: doubled history buffer so the dot product is one contiguous slice.
pub struct Fir {
    h_rev: Vec<f32>, // taps reversed
    hist: Vec<f32>,  // 2n samples
    pos: usize,
}
impl Fir {
    pub fn new(h: &[f64]) -> Self {
        let n = h.len();
        Fir { h_rev: h.iter().rev().map(|&v| v as f32).collect(), hist: vec![0.0; 2 * n], pos: 0 }
    }
    pub fn process(&mut self, x: f32) -> f32 {
        let n = self.h_rev.len();
        self.hist[self.pos] = x;
        self.hist[self.pos + n] = x; // mirror
        self.pos = if self.pos + 1 == n { 0 } else { self.pos + 1 };
        // newest sample is at pos-1+n; window of the last n samples starts at pos
        let w = &self.hist[self.pos..self.pos + n];
        w.iter().zip(&self.h_rev).map(|(a, b)| a * b).sum()
    }
}

/// One-pole lowpass (exponential smoother): y[n] = y[n-1] + a (x[n] - y[n-1]).
pub struct OnePole { a: f32, y: f32 }
impl OnePole {
    pub fn new(fc: f64) -> Self {
        // exact pole placement: pole p = exp(-2 pi fc), a = 1 - p
        OnePole { a: (1.0 - (-2.0 * PI * fc).exp()) as f32, y: 0.0 }
    }
    pub fn process(&mut self, x: f32) -> f32 { self.y += self.a * (x - self.y); self.y }
}

/// 2nd-order Butterworth lowpass (bilinear transform), transposed direct form II.
pub struct Biquad { b: [f64; 3], a: [f64; 2], z: [f64; 2] }
impl Biquad {
    pub fn butter_lp(fc: f64) -> Self {
        let k = (PI * fc).tan();
        let q = std::f64::consts::FRAC_1_SQRT_2;
        let n = 1.0 / (1.0 + k / q + k * k);
        let b0 = k * k * n;
        Biquad { b: [b0, 2.0 * b0, b0], a: [2.0 * (k * k - 1.0) * n, (1.0 - k / q + k * k) * n], z: [0.0; 2] }
    }
    pub fn process(&mut self, x: f64) -> f64 {
        let y = self.b[0] * x + self.z[0];
        self.z[0] = self.b[1] * x - self.a[0] * y + self.z[1];
        self.z[1] = self.b[2] * x - self.a[1] * y;
        y
    }
}

The tests (not shown) check unity DC gain, symmetric taps, 40+ dB attenuation well inside the stopband, impulse response equal to the taps, and the biquad's -3 dB point at fc and decaying impulse response. FIR group delay here is 25 samples (N=51).

Complexity and constants

• FIR: N multiply-adds and N loads per output; state 4N bytes (f32), doubled by the mirror buffer. Latency (N-1)/2 samples.
• IIR biquad: 5 multiplies, 4 adds, 2 state values; latency not constant.
• Measured here (Apple M4 Pro, release build, scalar Rust, order of): 51-tap f32 FIR 8 ns/sample, 255-tap 84 ns/sample, one f64 biquad 2.6 ns/sample. The biquad is a serial dependency chain, the FIR a dot product.
• Kaiser's estimate: N is about (A - 8) / (14.36 df) taps for A dB attenuation and normalised transition width df (order of).

Efficiency

• Bound: FIR is compute and load bandwidth (N MACs per sample); IIR is latency (the feedback loop cannot be pipelined across samples), not bandwidth.
• Biggest lever is choosing IIR when phase does not matter: about 5 MACs instead of tens to hundreds.
• Process blocks, not samples: a FIR over a block is a SIMD-friendly sliding dot product (CMSIS-DSP, NEON/AVX use this); for N above roughly 64 and long blocks, FFT convolution costs O(log N) per sample.
• Use symmetry: a linear-phase FIR needs only N/2 multiplies (pre-add x[n-k] + x[n-(N-1-k)]).
• Decimating filters (see sample rate conversion): compute only the outputs you keep, so cost drops by the decimation factor; half-band filters have about half of their taps zero.
• On small MCUs, f32 or Q15/Q31 integer arithmetic instead of f64 saves cycles and energy; I did not measure it.
• Multiple independent accumulators break the f32 sum dependency; I did not benchmark it.

Implementation notes

• FIR history as a doubled buffer (written at pos and pos+N) gives a contiguous slice for the dot product with no modulo and no branch in the hot loop.
• Reverse the taps once at construction so that the dot product runs forward over memory.
• IIR: transposed direct form II needs 2 states per biquad and has good numerical behaviour in float. Direct form I/II of order above 2 loses precision as poles cluster near 1: cascade second-order sections (SOS), pole pairs sorted by Q.
• Subnormals: decaying IIR state in float slows some CPUs 10-100x (order of). Flush denormals (FTZ/DAZ) or add a tiny DC offset.
• Filter start-up: initialise state to the first sample (steady state) to avoid a step transient; filtfilt pads and initialises states for the same reason (Gustafsson).
• Cut-off near 0 or near Nyquist needs f64 coefficients even if the filter runs in f32.
• Parameter changes at run time: smooth the coefficient change, or crossfade two filters, to avoid clicks.

Pitfalls

• fc is in cycles/sample (0-0.5) in the code above, not Hz or rad/sample: mixing the three is the usual bug.
• Even tap count gives half-sample delay and a forced zero at Nyquist, so it cannot be a high-pass or band-stop.
• Quantised IIR coefficients move poles; a pole just inside the unit circle may land outside and oscillate forever.
• Forgetting that filtering changes alignment: FIR delay (N-1)/2 must be compensated before comparing with the raw signal or detecting peaks.
• A single high-order IIR polynomial (above roughly order 6, order of) is numerically fragile in float; design as SOS.
• Zero-phase filtfilt is not causal and cannot run in real time; it also doubles the filter order in the magnitude response.
• Assuming "IIR is cheaper" at equal group-delay flatness: equalising IIR phase usually costs more than a FIR.

Variants and alternatives

In production

SciPy scipy.signal (firwin, butter, sosfilt, filtfilt), ARM CMSIS-DSP on microcontrollers (arm_fir_f32, arm_biquad_cascade_df2T_f32: FIR docs), audio EQ and crossovers (RBJ cookbook biquads), ECG/PPG preprocessing (PPG heart rate), software radios (digital down-conversion, sample rate conversion).

Implementations

Sources

• Alan V. Oppenheim, Ronald W. Schafer, Discrete-Time Signal Processing, Pearson, 3rd ed. 2010 (ch. 7).
• Julius O. Smith III, Introduction to Digital Filters with Audio Applications, W3K 2007. online
• T. W. Parks, J. H. McClellan, Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase, IEEE Trans. Circuit Theory 19(2), 1972. doi
• Fredrik Gustafsson, Determining the Initial States in Forward-Backward Filtering, IEEE Trans. Signal Processing 44(4), 1996. doi
• Robert Bristow-Johnson, Audio EQ Cookbook. W3C copy