Fast Cascaded Recursive Filtering via a Block-Matrix Reformulation

Recursive (IIR) filters realized as cascaded second-order sections (biquads) offer both design generality and robustness against coefficient quantization. However, their inherent sample-to-sample feedback dependency poses a fundamental obstacle to parallel computation. This paper reformulates the biquad difference equation as a banded block-Toeplitz linear system and introduces a stride-$N$ permutation that maps a group of $NL$ samples into a block-tridiagonal structure whose entries are scalar multiples of identity and shift matrices. Within this framework, two parallel algorithms are developed for the recursive solution: a partial LU (PH) factorization that preserves the sparse block structure and a cyclic reduction that is applied to recursive filtering, to the best of our knowledge, for the first time. It reduces the sequential dependency depth from $\mathcal{O}(N)$ to $\mathcal{O}(\log_2 N)$. For a cascade of $K$ biquads, the intermediate permutations between successive sections cancel exactly, so that only a single permutation/de-permutation pair is required for the entire cascade, eliminating $2(K{-}1)$ redundant stages. Exact block-level operation counts are derived for every algorithmic stage and validated against cycle-accurate measurements on three Intel micro-architectures supporting AVX2 SIMD instructions. Experimental results for a 16th-order system show that the proposed multi-block algorithms reduce clock cycles per sample by up to $10\times$ compared to scalar filtering, with both algorithms scaling favorably on newer architectures. On a single Meteor Lake core, cyclic reduction achieves approximately 618 MS/s -- an $8\times$ throughput improvement over scipy.signal.sosfilt.