Parallel Cascaded Recursive Filtering on Multi-Core CPUs and GPUs

The companion of this paper reformulated cascaded second-order (biquad) recursive filtering as a block-tridiagonal linear system and developed two parallel solution algorithms, PH factorization and cyclic reduction, reaching over 600 Megasamples per second on a single SIMD core. This paper scales that framework to multi-core CPUs and GPUs, where a new obstacle appears: the terminal outputs of each signal block group are the initial conditions of the next, so naively distributed groups serialize. The dependency is resolved by superposition -- each group's output splits into a zero-state response, computable immediately, and a homogeneous correction applied when the state arrives -- and by a divide-and-conquer form of cyclic reduction that exposes both terminal blocks before back substitution, as asynchronous state propagation requires. Two implementations pair the two dominant deployment scenarios with opposite treatments of the dependency. For real-time streaming, a wavefront pipeline realized with TBB flow graphs parallelizes across cascade sections, preserves first-in-first-out order, and achieves 3.95x scaling on six performance cores, about 2.4 Gigasamples per second for a 16th-order filter. For batched processing, a single-kernel GPU implementation carries each group through the entire cascade in registers and parallelizes across groups with a decoupled look-back protocol; a communication-based cost model, comprising a memory roof, a barrier price, and a latency-hiding floor, reduces tuning to two parameters and predicts the measured behavior across two GPU generations. The best kernels reach 38.2 Gigasamples per second for a single second-order section on an RTX~3060, 85% of the memory-bandwidth roof, exceed the strongest published parallel recurrence baseline at every filter order, and remain numerically valid at order 16, where the direct-form baseline fails.

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.