arXiv CSopen access

Parallel Cascaded Recursive Filtering on Multi-Core CPUs and GPUs

1

Parallel Cascaded Recursive Filtering on Multi-Core CPUs and GPUs

arXiv:2607.23763v1 [eess.SP] 26 Jul 2026

Haotian Zhai, Student Member, IEEE, and Bernd-Peter Paris, Senior Member, IEEE

Abstract—The companion paper of this two-part series 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.95× scaling on six performance cores, about 2.4 Gigasamples per second for a 16th-order filter. For batched processing, a singlekernel GPU implementation carries each group through the entire cascade in registers and parallelizes across groups with a decoupled lookback 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. Index Terms—IIR, recursive filter, cascaded second-order sections, multi-block filtering, parallel algorithms, cyclic reduction, decoupled lookback, multi-core CPU, GPU, latency hiding.

I. I NTRODUCTION

R

ECURSIVE filters are among the most fundamental building blocks in digital signal processing. Compared to non-recursive (FIR) filters, they typically require far fewer coefficients to meet a given magnitude-response specification, translating directly into fewer computations per unit time. This efficiency makes them especially attractive in highthroughput settings such as real-time digital filtering in modern communication receivers [1], speech and audio processing [2], high-resolution image restoration [3], and video processing [4]. However, the feedback dependency in a recursive filter, where each output sample depends on previous outputs, creates a sequential bottleneck that limits computational speed. Approaches to this bottleneck span half a century, from FFTbased approximations of the impulse response [5], [6] through block filtering [7], [8] to SIMD implementations on vector This work was supported in part by the National Science Foundation under grant 2029836. H. Zhai and B.-P. Paris are with the Department of Electrical and Computer Engineering, George Mason University, Fairfax, VA 22030 USA (e-mail: [email protected]; [email protected]).

DSP architectures and commodity CPUs [9], [10], [11], [12], [13]. In the first part of this two-part series [14], we addressed it for the most practically relevant realization—the cascade of second-order sections—through a block-matrix reformulation, multi-block filtering, in which a stride-N permutation maps a group of N L samples into a block-tridiagonal system solved by two parallel algorithms, PH factorization and cyclic reduction. On a single SIMD core, the resulting implementation reaches 616 MS/s for a 16th-order filter, an 8× improvement over the standard scipy.signal.sosfilt routine [15]. A single core, however, is no longer where high-throughput computing resides: commodity CPUs integrate a dozen or more heterogeneous cores, and GPUs keep tens of thousands of threads resident to hide latency [16]. Scaling multi-block filtering across such units meets an obstacle that does not exist within one core: the inter-group state dependency. The last two outputs of each signal block group are the initial conditions of the next group in every section [14], so a naive distribution of groups across processing units serializes them. Moreover, two different deployment scenarios impose different disciplines: real-time streaming, as in a software-defined radio receive chain, requires low latency and first-in-first-out output order, while batched processing of recorded data requires maximum throughput and imposes no order. The two scenarios favor different architectures, and, as this paper shows, different algorithmic treatments of the same dependency. Prior work offers partial answers. Building on the classical prefix-scan constructions [17], [18], the modern single-pass GPU scan resolves the analogous dependency of first-order recurrences through decoupled lookback, in which each block publishes a local result and inspects its predecessors’ published states [19]. Maleki et al. extended scan-based evaluation to general mth-order linear recurrences [20], but in the directform realization, whose per-thread state grows with the order and whose coefficients are numerically fragile for high-order filters [21]. GPU implementations of block IIR filtering launch separate kernels per computation stage, with every intermediate result passing through global memory [22], and GPU recursive filtering for image processing has relied on overlapped blocks tailored to the two-dimensional setting [3]. On CPUs, taskbased runtimes such as Threading Building Blocks provide pipelined parallelism across cores [23], but published parallel IIR work targets the direct form rather than the cascade. No existing work, to our knowledge, integrates inter-group state propagation with group-level parallel filtering algorithms for the cascaded second-order realization on either architecture. This paper fills that gap. We extend the multi-block filtering framework of the first part [14] to multi-core CPUs and GPUs, covering both deployment scenarios. The main contributions are as follows. 1) Inter-group state propagation: We derive the state-

2

bearing forms of PH factorization and cyclic reduction: of section k+1. The cascaded realization confines the effect by superposition, each group’s output splits into a zero- of coefficient quantization to a single pole pair per section and state response, computable before the group’s state is is numerically more robust than the direct form. known, and a homogeneous correction applied when the state arrives, allowing groups to execute concurrently. B. Multi-Block Filtering 2) Divide-and-conquer cyclic reduction: We reformulate The input sequence is partitioned into blocks of L samcyclic reduction so that both terminal blocks of a group ples, and N consecutive blocks form a signal block group are available before back substitution—the property that of N L samples. Within a group, the samples are perasynchronous state propagation requires. muted with stride N : the m-th permuted block is x̄m = 3) Real-time multi-core CPU implementation: A wave⊤ front pipeline realized with TBB flow graphs preserves [x[m], x[m+N ], . . . , x[m+(L−1)N ]] . Multi-block filtering FIFO order and achieves 3.95× scaling on the six describes the action of a single second-order section (1) on performance cores of a heterogeneous Meteor Lake a signal block group: as shown in [14], the permuted group satisfies the linear system processor, 2.4 GS/s for a 16th-order filter. 4) Batched GPU implementation: A single kernel carries A y = B x + B− x−1 − A− y−1 , (2) each group through the entire cascade in registers, resolving the state dependency with decoupled lookback; where x = [x̄0 , . . . , x̄N −1 ]⊤ , y is defined analogously, and a communication-based cost model reduces tuning to two x−1 = [x[−2], x[−1]]⊤ , y−1 = [y[−2], y[−1]]⊤ are the last parameters and predicts the measured behavior across two input and output samples of the preceding group. The two GPU generations. The best kernels reach 38.2 GS/s matrix A ∈ RN L×N L is block-Toeplitz with corner coupling, for a single second-order section on an RTX 3060, 85%   I a2 R a1 R of the memory-bandwidth roof. a1 I I a2 R 5) Cross-architecture evaluation: The cascade algorithms   a2 I a1 I  I A= (3) exceed the strongest published GPU recurrence en,   . . . . . . gine [20] at every filter order and remain numerically   . . . valid at order 16, where the direct form fails. a2 I a1 I I This paper is organized as follows. Section II reviews the multi-block filtering framework of Part I and extends it where every L × L entry is a scalar multiple of the identity with the two ingredients that parallel execution requires: the matrix I or the lower shift matrix R. The matrix B has theNsame L×2 divide-and-conquer form of cyclic reduction and the inter- structure with a1 , a2 replaced by b1 , b2 , and A− , B− ∈ R group state propagation algebra. Section III describes the real- couple the initial state to the first blocks. The computation of (2) can be separated into two stages. time multi-core CPU implementation: wavefront execution, its The non-recursive stage evaluates realization with TBB flow graphs, and measurements on a heterogeneous Meteor Lake processor. Section IV describes v = B x + B− x−1 (4) the batched GPU implementation: the thread-block mapping, the decoupled lookback protocol, a communication-based cost in parallel across all N blocks. The recursive stage complemodel, the four algorithms, and measurements on two GPU ments A− y−1 and solves generations. Significant throughput improvements over the A y = v − A− y−1 . (5) strongest published parallel recurrence baseline are shown. Section V presents a summary and conclusion as well as future Note that this differs slightly from the theory paper, where work. the initial-state term A− y−1 was folded into the right-hand side of (4). The reason is that the theory paper focuses on II. A RITHMETIC BACKGROUND single-core execution: serial group execution guarantees that This section summarizes the multi-block filtering framework y−1 is available when a group starts. This paper targets parallel developed in Part I [14] and extends it with the two ingredients group execution on multi-core CPUs and GPUs, and y−1 holds that parallel group execution requires: the divide-and-conquer the last two output samples of the preceding group. This means formulation of cyclic reduction and the principle of inter-group that A− y−1 carries the inter-group dependency. The non-recursive stage is computed trivially in parallel; the state propagation. recursive stage is the subject of the remainder of this section. A. The Cascaded Recursive Filter A general recursive filter of order 2K is realized as a cascade of K second-order sections (biquads). Each biquad computes y[n] = x[n] + b1 x[n−1] + b2 x[n−2] − a1 y[n−1] − a2 y[n−2], (1) where the scaling coefficients b0,k of all sections are split off into a common gain factor. The output of section k is the input

C. PH Factorization The PH factorization A = P H decomposes A into a banded block-Toeplitz factor P without corner coupling and a sparse  factor H [14], Eq. (23) . Substituting A = P H into (5) and moving A− y−1 to the left-hand side yields
P H y + P −1 A− y−1 = v. (6)

3

Defining

ye [7]

w =Hy+P

−1

A− y−1 ,

(7)

⊤

ye [4] C2

C

ye [2]

C

C

ye [1] C

we [7]

C

we [6]

we [5]

C

we [4]

we [3]

C

we [2]

we [1]

ye [0] C

m = 0, . . . , N −1, (9)

u1,N −1 e0

L

where e0 = [1, 0, . . . , 0] ∈ R and the coefficients ul,m follow the  same recurrence as those in H12 of H [14], Eq. (21) . Combining (7) with (10), the first N −2 blocks satisfy ȳm + u2,m ȳ−2 + u1,m ȳ−1 = w̄m ,

ye [5]

C2

(8)

with w̄−1 = w̄−2 = 0. Since each block requires the two preceding blocks, the forward substitution has sequential depth O(N ). Note that neither v nor (9) involves y−1 : the particular solution of every group can be computed independently. Recovering the output. The matrix P −1 A− in (7) can be derived as the explicit form   u2,0 e0 u1,0 e0  u2,1 e0 u1,1 e0    P −1 A− =  (10) , .. ..   . . u2,N −1 e0

C3

C2

for w, then recover y from (7). Particular solution. Since P is banded block Toeplitz, (8) is solved by the forward substitution w̄m = v̄m − a1 w̄m−1 − a2 w̄m−2 ,

ye [6]

ye [3]

the solution proceeds in two stages: first solve Pw=v

C4

m = 0, . . . , N −3, (11)

we [0]

y−1

Fig. 1. Sklansky construction [24] for the matrix scan (13) with L = 8. The initial state y−1 of the group enters the network at ye [0].

the homogeneous solution (11) follow the particular solution directly. Under parallel group execution, however, every group other than the first obtains y−1 from the recursive doubling of its predecessor: later groups cannot complete their homogeneous solution until the terminal states of all earlier groups are known. The terminal states must therefore propagate across the groups first; once this inter-group state propagation completes, the homogeneous solutions of all groups are independent and proceed fully in parallel. We will discuss the principle of inter-group state propagation later in Section II-F. D. Cyclic Reduction

Cyclic reduction in [14] solves (5) by eliminating alternating where ȳ−2 = [y[−2], y[N −2], . . . , y[(L−1)N −2]]⊤ and ⊤ blocks. It proceeds in three steps: reduction, terminal block ȳ−1 = [y[−1], y[N −1], . . . , y[(L−1)N −1]] , i.e., the termifiltering, and back substitution. nal blocks shifted by one position with the initial state in the Reduction. Each reduction level removes the even-indexed first entry. The last two blocks satisfy blocks of the current system and produces a system of half the        I ye [0] C we [0] size with the same second-order form. The elimination of the C I   ye [1]   0   we [1]  left-hand side A is developed in Part I [14]; here, the right-hand           +  ..  y−1 =  , .. .. .. .. side additionally carries the initial-state term A− y−1 of (5).         . . . . . After the same sequence of permutations and eliminations, the C I ye [L−1] 0 we [L−1] system at level i becomes (12) ⊤ where ye [l] = [y[(l+1)N −2], y[(l+1)N −1]] , we [l] is de(i) (i) A(i) = v(i) − A− y−1 , (14) CR y fined analogously, and the entries of the 2 × 2 matrix C are (i) precomputed from the filter coefficients [14], Eq. (30) . where A− is the initial-condition matrix at reduction level i, Recursive doubling. Each row of (12) is the first-order matrix given by  (i)  recurrence h e0 g (i) e0  d(i) e0    ye [l] = we [l] − C ye [l − 1], l = 0, . . . , L−1, (13)     (i) ,  (15) A− =  (i+1) with ye [−1] = y−1 . The recurrence is solved by recursive (i+1)  h e g e 0 0   doubling with the Sklansky construction [24], which takes (i+1)  e0  f   log2 L levels. Figure 1 illustrates this construction for L = 8, with the initial state y−1 entering the network at ye [0]. The efficient computation of the recursive doubling structure shown and d(i) is the updated coupling coefficient of the even-indexed in Figure 1 is developed in the Part I [14]. (i) . The lower half of A− carries over to Homogeneous correction. With the terminal blocks and blocks [14], Eq. (41) (i+1) (i) the initial state (y−1 ) known, the remaining N −2 blocks are the next level as A− , mirroring the reduction of ACR itself. This pattern continues until the final level i = log2 N − 1, recovered from (11) independently of signal groups. where Note that in Part I [14], since the focus is single group  (i)  algorithm computation, we assume the initial state y−1 is h e0 (g (i) + d(i) R)e0 (i) A− = (i+1) . (16) known when a group starts. Then, recursive doubling and h e0 (g (i+1) + f (i+1) R)e0

4

The initial-condition coefficients are updated alongside the reduced coefficients e(i) and f (i) , with e(0) = a1 and f (0) = a2   [14], Eq. (41) , as h(0) = a2 , h(i+1) = −e(i) h(i) ,

g (0) = a1 , g (i+1) = f (i) − e(i) g (i) .

(17)

All coefficients depend only on the filter parameters and are precomputed. The reduction proceeds over log2 N levels; all block updates within a level are independent. Terminal block filtering. At the deepest level, a single block remains: the coupling of ȳN −1 to all other blocks has been eliminated, and the lower row of (16) yields
(log N ) I+e(log2 N ) R+f (log2 N ) R2 ȳN −1 =v̄N −12 −h(log2 N ) e0 y[−2]
− g (log2 N ) +f (log2 N ) R e0 y[−1],

(18) (log N ) where v̄N −12 is the accumulated right-hand side of the last block after all log2 N reduction levels. Since the left hand side in (18) is a unit banded lower-triangular Toeplitz matrix, it can be solved by block filtering. Back substitution. With y−1 = [y[−2], y[−1]]⊤ known, the back substitution recovers the eliminated blocks level by level in reverse order of the reduction. All recoveries within a level are independent, and the even-indexed blocks of the original system—including ȳN −2 —are recovered at the last level. Cyclic reduction requires approximately 6N block FMAs per group, with sequential depth O(log2 N ) in both the reduction and the back substitution. Owing to this shorter dependency chain compared to the O(N ) particular solution of the PH factorization (9), the SIMD implementation of cyclic reduction achieved the best throughput on a single CPU core shown in [14]. E. Divide-and-Conquer Cyclic Reduction

retained input blocks and has the same block structure as A(1) CR , (1) (1) with coefficients eB and fB ; the cross matrix B(1) cross carries the contribution of the eliminated input blocks, with coefficients (1) (1) cB and dB :   (1) (1) I fB R eB R  (1)  (1)  eB I I fB R    (1) Bmain =  (1) . . ,  fB I . . . .    . . (1) . eB I I (20)  (1)  (1) dB I cB R  (1)  (1)  cB I dB I  (1)  . Bcross =   .. ..   . . (1)

(1)

cB I dB I Symmetrically, eliminating the odd-indexed blocks yields the even-indexed subsystem (1) e (1) A(1) CR y even = Bmain xeven + Bcross xodd ,

(21)

with yeven = [ȳ0 , ȳ2 , . . . , ȳN −2 ]⊤ . The left-hand-side matrix and the main matrix are identical to those in (19), since both reductions apply the same Schur complement to A and B. The (1) (1) cross matrix uses the same coefficients cB and dB , with the banding shifted by the reversed elimination direction,   (1) (1) cB R dB R  (1)  (1)  dB I cB R  (1)   e B (22) . cross =  (1)  cB I . . .    . . (1) . dB I The feedforward coefficients after the first reduction level are (0)

(0)

eB = b1 The cyclic reduction in [14] is one-sided: each level elimi- fB = b2 (1) (1) (0) (0) (0) (0) nates one parity and retains the other, so the eliminated blocks fB = f fB eB = fB + f (0) − e(0) eB are recovered only during back substitution. In particular, the (1) (0) (0) (1) (0) cB = f (0) eB − e(0) fB dB = eB − e(0) terminal state of the group is split across the two parities: (23) ȳN −1 emerges at the deepest reduction level, but its companion where e(0) = a and f (0) = a . The initial-condition and 1 2 ȳN −2 is recovered only at the last substitution level. Cyclic wrap-around terms reduce alongside, through initial-condition reduction is, at its core, a divide-and-conquer method [25], and matrices of the same form as (15), and are omitted here for the classical device for extracting more parallelism from it is clarity. to apply the elimination around both parities, trading redundant In this work, one level of divide-and-conquer is performed. arithmetic for a second, independent subproblem. The two subsystems (19) and (21) are independent, and each is Unlike the preceding algorithms, the divide-and-conquer solved by the one-sided cyclic reduction on N/2 blocks, starting variant does not evaluate the non-recursive stage (4) separately. from level i = 1. When two processors are available, the two Instead, the elimination is applied to both sides of (2), reducing subsystems proceed concurrently, which is the first benefit A and B together, so that the non-recursive computation is of the divide-and-conquer form: it doubles the parallelism of folded into the reduction steps. Grouping the blocks by parity every reduction and substitution level. and eliminating the even-indexed blocks yields the odd-indexed The second benefit concerns the terminal state. The odd subsystem subsystem exposes ȳN −1 at its deepest reduction level, and (1) (1) A(1) (19) the even subsystem exposes ȳN −2 . Both terminal blocks—and CR y odd = Bmain xodd + Bcross xeven , with them the terminal state [y[N L−2], y[N L−1]]⊤ of the ⊤ where yodd = [ȳ1 , ȳ3 , . . . , ȳN −1 ] and xodd , xeven collect the group—are therefore available before any back substitution (1) odd- and even-indexed input  blocks. The matrix ACR is the begins. This property is essential for parallel group execution: level-1 reduced system of [14], Eq. (40) . On the right-hand the terminal states can first be propagated across all groups, side, the main matrix B(1) after which the back substitution of every group depends only main carries the contribution of the

5

on data within the group, and all groups proceed independently in parallel. The added parallelism is paid for in arithmetic: reducing both parities requires approximately 4N block FMAs per direction and 8N in total, against 6N for the one-sided variant, with the block shuffle count roughly doubled as well. Divide-andconquer cyclic reduction is therefore worthwhile exactly when parallel resources can absorb the redundancy, as in the batched GPU execution of Section IV. F. Inter-Group State Propagation

The upper level admits the same zero-state treatment in turn. Unrolling (25) expresses each y−1 [g] as a combination of the zero-state terminal vectors t[0], . . . , t[g−1] and the initial state y−1 [0] of the first group alone. Consequently, once the zerostate terminal vectors of all groups and the initial condition of the first group are known, the true initial condition of every group can be computed independently—sequentially, or by one more recursive doubling across the groups. The zero-state terminal vector of a group is commonly referred to as its local carry, and the resolved initial condition as its global carry: local carries are computed fully in parallel, and global carries are resolved by a lightweight propagation whose cost is independent of the group size. The same strategy underlies the GPU implementation discussed in Section IV.

The algorithms above compute a single signal block group, and the initial state y−1 of each group is the terminal output of its predecessor. The straightforward execution is therefore sequential: the groups are processed in order, and the overall III. R EAL -T IME P ROCESSING ON M ULTI -C ORE CPU S rate is set by the per-group throughput of the chosen algorithm, The single-core benchmarks in Part I [14] established the best which is the basis of comparison in Part I [14]. Parallel achievable throughput for each computational kernel. This execution across groups requires breaking this chain: the section composes those optimized kernels into a multi-core initial conditions of all groups must be resolved before, or implementation targeting real-time signal processing, where independently of, the bulk of the per-group computation. We input data must be processed in FIFO order with minimal refer to this as inter-group state propagation. latency. We first present the wavefront execution strategy, then The scan (13) includes the initial state ye [−1] = y−1 , which its realization with TBB flow graphs, and finally the multi-core depends on the terminal output of the preceding group. Two measurements on Meteor Lake. strategies exist for incorporating this term, shown in Figure 1. The first includes y−1 at the initial level, so that the recursion propagates its contribution through all subsequent levels. The A. Wavefront Execution second sets ye [−1] = 0 and runs the recursive doubling We adopt a wavefront execution strategy for multi-core CPU independently of the initial state, producing the zero-state implementations of cascaded IIR filters. Figure 2 illustrates response. Since (13) is linear, the contribution of y−1 can be the execution pattern for a cascade of K second-order sections. superposed afterward as Each circle represents one section applied to one signal block group; each column is a group, and colors indicate thread ye [l] = ye,0 [l] − C l+1 y−1 , l = 0, 1, . . . , L−1, (24) assignment. Within a group, the signal flows vertically down where ye,0 [l] denotes the zero-state output. Both strategies have the cascade; between adjacent groups, each section passes its the same computational complexity, but the second is preferable output state along the gray diagonal to the next group. Time for parallel architectures: the zero-state recursive doubling advances downward, so all nodes at the same height execute requires no inter-group information and can be executed concurrently—one wavefront. Wavefront execution possesses several properties that make independently across all groups, with the correction (24) applied it well suited for real-time processing on multi-core CPUs. once the initial states are resolved. The same decomposition First, the pipelined cascade structure preserves FIFO ordering applies to cyclic reduction, where the initial state enters the by construction: each signal block group enters and exits the terminal solve (18) only through the explicit precomputed cascade in the same order. Second, the same signal block terms. group is processed by a single thread across all secondThe structure of Figure 1 makes the resolution of the initial order sections, maximizing data locality and minimizing interstates themselves equally simple. Setting l = L−1 in (24), core context switching. Third, within each section the groups the terminal state of a group is its zero-state terminal vector L execute strictly in order, so the initial state is always available corrected by −C y−1 . To track this across groups, index when a group arrives: the serial single-group functions of the groups by g and append the index as a second argument: [14] apply unchanged, and no inter-group state propagation is ye,0 [l, g] denotes the zero-state output of group g, and y−1 [g] required. A downstream section waits only for the state of its its initial state. Writing t[g] = ye,0 [L−1, g] for the zero-state predecessor, not for the entire cascade, so the pipeline fill depth terminal vector of group g, the initial states of consecutive is proportional to the number of sections and the steady-state groups therefore satisfy throughput scales with the number of available cores. y−1 [g+1] = t[g] − C L y−1 [g], (25) which is a first-order matrix recurrence of exactly the same form as (13), one level up: the groups themselves form an upper-level recursive doubling, with C L in place of C and the zero-state terminal vectors in place of we [l].

B. Implementation with TBB Flow Graphs The wavefront execution is implemented with the Threading Building Blocks (TBB) [23] flow graph, a task-based runtime for expressing pipelined parallelism without explicit thread

6

signal group / core group 0 F1 group 1

filter order / time

F2

group 4 FK

Fig. 4. Node-level execution timeline of a 16th-order IIR filter on 6 Meteor Lake P-core threads, captured from TBB Flow Graph Analyzer. Each color indicates a different thread. The steady-state pipelining pattern and the approximately 150 ns inter-node overhead are visible.

Fig. 2. Wavefront execution of a cascaded IIR filter with K second-order sections. Each circle is one section applied to one signal block group; colors indicate thread assignment, with one thread following each group through the entire cascade. The signal flows vertically through the cascade, section states pass diagonally to the next group, and nodes at the same height execute concurrently.

f ()

f (x)

f (x)

source

permute

iir0

···

f (x)

f (x)

f ()

iir7

deperm.

sink

Fig. 3. TBB flow graph for a 16th-order cascaded IIR filter. Each node has serial concurrency to enable wavefront pipelining.

management. The flow graph is a linear chain of second-order filter nodes, as shown in Figure 3, each configured with serial concurrency: only one thread can execute a given node at any time, while different nodes execute concurrently across threads, realizing the pipeline of Figure 2. The backbone of TBB is its task-queue-based scheduler: each worker thread maintains a local task queue, allowing it to carry the same signal block group through consecutive cascade stages whenever the next stage is immediately available. This behavior aligns naturally with the wavefront pattern. C. Function Size and Grain Size Figure 4 shows the steady-state execution timeline of a 16th-order filter using the CR algorithm under the best singlecore configuration (L = 8, N = 64) on six P-core threads, captured with the TBB Flow Graph Analyzer. A clear pipelining pattern emerges: each thread consistently processes the same signal block group, and approximately four filter nodes execute concurrently across threads. A key observation from Figure 4 is that the execution gap between adjacent nodes of the same thread—the task-generation overhead of the runtime—is approximately 150 ns. For CRbased multi-block filtering with a function size of N L = 512

samples, the per-section latency on a Meteor Lake P-core at 4.5 GHz is 7.3 cycles/block 512 samples · ≈ 103 ns, 8 samples/block 4.5 × 109 cycles/s where 7.3 cycles per block of L = 8 samples is the per-section single-core measurement of Part I [14]. Since this execution time is comparable to the 150 ns task-generation overhead, the optimized single-core function size is too small for efficient per-node execution in the task system. To address this, we introduce an additional parameter, the grain size—in TBB terminology [23]—which specifies how many signal block groups are aggregated and processed within a single node invocation. A grain size of 256 increases the per-node execution time to 103 × 256 ≈ 26,400 ns ≫ 150 ns, effectively amortizing the task-generation overhead. The grain size thus controls task granularity, at the cost of a proportional increase in per-group latency. Function size and grain size together form a two-level strategy: the function size determines the efficiency of processing each signal block group on a single core, while the grain size determines the task granularity needed to hide the overhead of the runtime scheduler. D. Heterogeneous Cores Meteor Lake-H features a hybrid architecture with six hyperthreaded performance cores (P-cores) and eight efficiency cores (E-cores) [26]. E-cores support only 128-bit vector registers, so an L = 8 single-precision block requires two vector instructions instead of one on a P-core; they also run at a lower clock, 3.8 GHz against the P-cores’ 4.5 GHz. Since the throughput of a pipeline is determined by its slowest stage, Ecore nodes become the bottleneck. Figure 5 shows the execution timeline when two E-cores are introduced alongside six P-cores: the slower E-core execution, visible as longer segments for iir4 – iir7 , creates pipeline bubbles that reduce aggregate throughput below the six-P-core configuration. E. Measurement and Results The multi-core experiments target the Meteor Lake-H platform of Part I [14]: six hyperthreaded P-cores at 4.5 GHz

7

TABLE I T HE TWO TEST DEVICES . T HE DERIVED ROWS FOLLOW FROM (28) AND (29). GTX 1070 Architecture SMs Threads / warps / TBs per SM Shared memory per SM / per TB Registers per SM / per thread FP32 cores per SM / FMA latency Memory bandwidth Throughput roof (28) Hiding floor Wmin (29)

Fig. 5. Execution timeline with 8 threads (6 P-cores + 2 E-cores). The slower E-core nodes (longer segments) become the pipeline bottleneck.

RTX 3060

Pascal GP104 Ampere GA106 15 28 2048 / 64 / 32 1536 / 48 / 16 96 / 48 kB 100 / 48 kB 65536 / 255 65536 / 255 128 / 6 cyc 128 / 4 cyc 256 GB/s 360 GB/s 32 GS/s 8 warps

45 GS/s 6 warps

Third, throughput scales with the thread count and saturates at six P-cores, where the implementation reaches 3.95× at grain size 512. Fourth, adding E-cores degrades performance: the eight-thread configuration (one E-core per cluster) and the tenthread configuration (two per cluster) fall below the six-P-core curve at every grain size, as the slower E-core stages bottleneck the pipeline (Figure 5). The ten-thread configuration degrades at smaller grain sizes than the eight-thread one because pairs of E-cores sharing the cluster L2 reach the effective per-core cache capacity earlier. In absolute terms, the six-P-core configuration at grain size 512 achieves 3.95 × 616 MS/s ≈ 2.4 GS/s for the 16th-order system, where 616 MS/s is the single-core CR throughput of Part I—itself an 8× improvement over scipy.signal.sosfilt [15]. IV. BATCHED P ROCESSING ON GPU S

Fig. 6. Normalized multi-core throughput of the 16th-order CR-based IIR filter on Meteor Lake versus grain size, for thread counts from 1 to 10. Throughput is the ratio to the single-core throughput (616 MS/s) of Part I.

and eight E-cores at 3.8 GHz [26]. Each P-core has 2 MB of private L2 cache; each cluster of four E-cores shares 2 MB. The single-core optimized kernels from Part I, implemented in C++ with the Vector Class Library [27], are composed into the TBB flow graph of Figure 3 with the CR-based multi-block filtering algorithm (L = 8, N = 64). Threads are pinned to P-cores first; E-cores are used only when the thread count exceeds six. Throughput is measured by processing a sufficiently large number of signal block groups to reach steady state, and is reported normalized to the single-core baseline of Part I. Figure 6 presents the normalized throughput as a function of grain size for one to ten threads, with the largest value of each configuration marked next to its curve. Four trends emerge. First, the single-thread curve approaches but does not reach the baseline, saturating at 0.87; the remaining gap is the TBB runtime overhead. Second, throughput increases with grain size up to 512 as the task-generation overhead is amortized, then drops sharply at 1024: the working set of 1024 groups of N L = 512 single-precision samples occupies 1024 × 512 × 4 B = 2 MB, reaching the private L2 capacity.

The multi-core implementation of Section III achieves lowlatency streaming by pipelining signal block groups across a small number of powerful cores. A GPU offers the opposite resource: thousands of simple execution lanes whose collective throughput is large, but only for work that decomposes into many independent pieces. This section maps the multi-block filtering algorithms onto that architecture and reports the measured results. An NVIDIA GPU is an array of streaming multiprocessors (SMs). Threads are scheduled in warps of 32 that execute in lockstep; warps are grouped into thread blocks (TBs) that share a fast on-chip memory (shared memory) and synchronize through barriers; TBs are launched as a grid and cannot synchronize with each other except through global (DRAMbacked) memory. When a warp stalls, the SM switches to another resident warp at zero cost: latency is not avoided but hidden by oversubscription [16]. Arithmetic is abundant; what costs is communication, at three price levels—register shuffles inside a warp, barrier-synchronized shared memory across warps, and global-memory round trips across TBs. The experiments use the GTX 1070 (Pascal) and the RTX 3060 (Ampere); Table I collects the hardware parameters used throughout. A. Thread-Block Mapping and Shared-Memory Staging The input stream is divided into groups of C = N L samples, and each group is assigned to one thread block. Within its TB,

8

L the group is carried through all K sections of an order-2K filter: the output of section s, resident in registers, is the input Global memory: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 · · · of section s+1, so the group is read from global memory once, (i) coalesced load N +1 filtered with K second-order sections on chip, and written Shared memory: 0 1 2 3 4 5 6 7 once. This single-kernel, register-resident cascade contrasts 8 9 10 11 12 13 14 15 L 16 17 18 19 20 21 22 23 with prior GPU work on multi-block recursive filtering [22], 24 25 26 27 28 29 30 31 which launches separate kernels per stage—at least 2K+1 (ii) row read (padded banks) launches for a K-section cascade—with every intermediate result passing through global memory. Each TB learns which Register file: T0 : 0 1 2 3 4 5 6 7 group it owns by drawing a ticket from a global atomic counter (iii) all K sections T1 : 8 9 10 11 12 13 14 15 T2 : 16 17 18 19 20 21 22 23 in registers at launch; tickets ascend, so a TB’s predecessors are always T3 : 24 25 26 27 28 29 30 31 scheduled no later than itself, a property the state protocol of (iv) row write-back inter-block lookback (Section IV-B) relies on. The multi-block algorithms operate on stride-N perShared memory: 0 1 2 3 4 5 6 7 muted data: thread i holds the N consecutive samples 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 x[iN ], . . . , x[iN +N −1], so that at block step m the warp’s 24 25 26 27 28 29 30 31 L lanes collectively hold the permuted block {x[m+iN ]}i . (v) coalesced store Issued directly against global memory, such strided access Global memory: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 · · · would fragment into up to 32 separate transactions, because the memory system serves contiguous 32-thread transactions only. Fig. 7. Memory-access pattern at group entry and exit (L=4 threads, N =8 The implementation therefore never lets a strided address reach blocks shown); steps (i)–(v) are described in the text. DRAM. Figure 7 shows the pattern at every group boundary: (i) the TB loads its C samples contiguously into a sharedmemory array of L rows and N +1 columns; (ii) each thread (the gray arrows in Figure 8); if predecessor c−j has posted reads its row into registers, kept free of bank conflicts by the its global carry (its resolved terminal state), the walk stops, one-word row padding; (iii) all K sections execute in registers; and the TB composes that global carry with the local carries (iv) results are staged back row-wise; and (v) committed to of the intervening groups through the 2 × 2 matrix action of a global memory contiguously. With single-precision (4-byte) full group, resolving its own global carry (the small rectangles samples, the buffer costs in Figure 8). The lookback distance is therefore variable: the TB waits only for its most recent finished ancestor, however Sbuf = 4L(N + 1) bytes per TB, (26) far back, not for its immediate predecessor to finish. Finally, which will later be shown to be the limiting resource of the the TB applies the state correction—for PH, the homogeneous configuration space (Section IV-D). Because the cascade is solution (11)—publishes the resolved carry, and proceeds to fused, the staging is paid twice per group—at entry and exit— section s+1. The protocol’s overhead is a few operations per regardless of K, so the per-section staging cost falls as 1/K; predecessor inspected against O(C) useful work per group, and this is the same amortization the CPU implementation obtains the wait itself is hidden by the SM scheduling other resident groups. from its shuffle-based permutation. The CPU and GPU split on this choice for structural reasons. The wavefront of Section III parallelizes across sections: at B. Decoupled Lookback for State Propagation most K groups are in flight, the state hand-off is free, FIFO Groups are independent except for the filter state: section s order is preserved, and the cheaper one-sided cyclic reduction of group c needs the last two outputs of section s of group c−1 applies. This matches a streaming CPU, where the core count as its initial conditions. Left untreated, this dependency chains is close to K and output order matters. The GPU must instead every group to its predecessor and the grid degenerates into a fill 15 × 2048 = 30,720 (GTX 1070) or 28 × 1536 = 43,008 serial pipeline. The decoupled lookback protocol [19], extended (RTX 3060) resident threads to hide latency; the section axis to higher-order recurrences in [20], removes the serialization offers few units (usually K ≤ 8), so only the group axis is by exploiting linearity, exactly as in Section II-F: because the long enough, and batched processing supplies the required initial state enters additively, a group computes almost all of backlog while imposing no output order. Lookback is therefore its output before its state is known. necessary on the GPU rather than merely preferred. Its price Figure 8 shows the execution. The protocol runs once per is paid in Section II-E’s currency: publishing the complete section. For section s, the TB first performs the local solve: terminal state early forces cyclic reduction into its divide-andit filters the section under zero initial conditions, which is conquer form. the overwhelming majority of the arithmetic and requires no communication; the byproduct is the section’s local carry, its terminal output pair under zero state. The TB publishes the local C. The GPU Cost Model Part I [14] ranked the algorithms by operation count, the carry to a global descriptor array, issues a memory fence, and sets a status flag. It then performs the lookback: one designated right measure for a processor whose bottleneck is arithmetic warp walks backward through the predecessors’ descriptors latency along a dependency chain. The GPU prices the same

9

signal group / thread block group 0

group 1

group 4

F1 DRAM

filter order / time

F2

FK

Fig. 8. Decoupled lookback execution of the same cascade as Figure 2; one thread block carries one group (color) through all sections. Local carries are written to, and read from, a descriptor array in DRAM shared by all thread blocks (gray arrows); the small rectangles resolve each group’s global carry before the next section.

dependency graphs differently: operations are nearly free, because each warp’s arithmetic latency is hidden by switching to another resident warp, while every dependency whose producer and consumer sit in different threads must execute as an explicit exchange. A register shuffle moves a value between lanes of one warp and is issued and hidden like arithmetic. A cross-warp exchange is a store to shared memory, a block-wide barrier, and a load; the barrier stalls every warp of the TB, and warp switching cannot hide it, because the waiting warps are the ones that would provide the hiding. A cross-TB exchange is a global-memory round trip. Following the standard cost models of communication-avoiding linear algebra [28], BSP synchronization [29], the roofline bound [30], and latencyhiding analysis [16], [31], we model the steady-state time per sample as F t ≈ tmem + + β σ + εhide . (27) C Memory floor. Every algorithm reads each input sample once and writes each output once; intermediates stay on chip and the carries add a negligible O(K/C) bytes per sample. At single precision this is bs = 8 bytes of irreducible DRAM traffic per sample, so tmem = bs /BW, (28) giving throughput roofs of 256/8 = 32 GS/s (GTX 1070) and 360/8 = 45 GS/s (RTX 3060). No kernel can exceed its device’s roof, and the experimental measurements (Section IV-F) normalize every plateau by it. No arithmetic term. A cascade performs roughly 6K FMAs per sample (PH; 8K for STCR) against the fixed 8 bytes, an arithmetic intensity of 1.5K–2K FLOP per byte. The devices’ balance points (peak FLOP/s over bandwidth) sit near 25 and 35 FLOP per byte, so the workload is memory-bound everywhere and the arithmetic term is dropped. The measurements test this

deletion directly, and locate its boundary where the margin is thinnest—order 16 on the GTX 1070. Fixed per-group overhead. F collects the constant work a group performs because it has boundaries: the ticket draw, reading and injecting carries, publishing with fences and flags, and the entry/exit staging synchronizations. It is independent of C, hence the dilution F/C: negligible at the plateau, visible at small batch sizes. Barrier term. σ counts TB barriers per sample—an algorithm and configuration property—and β is the machine’s cost per barrier: an architecture property. The separation lets the measurements attribute a gap to the algorithm on one device and to the architecture between devices. Exposed latency. εhide is the latency that resident parallelism fails to cover [16]. Keeping a pipeline of latency λ cycles busy at an issue rate of I warp instructions per cycle requires Iλ instructions in flight per SM. The supply has two factors. Across warps, every resident warp (resident TBs × warps per TB) is an independent instruction stream. Within a warp, instructions issue back to back until one needs a result still in the pipeline, so the number a warp keeps in flight—its instruction-level parallelism (ILP)—is set by the dependency structure of its inner loop: how many instructions the loop issues between consecutive dependent ones. ILP is therefore a property of the kernel’s code, and it differs across kernels. A single floor follows from the least-parallel loop in the suite, the per-thread particular solution (9): of its three instructions per iteration—two FMAs and one load—only one FMA depends on the preceding iteration, so three instructions issue per dependent step and ILP ≈ 3. Every other loop offers at least this much. Setting supply equal to demand, the minimum resident warps per SM are Wmin = ⌈I λ/ILP⌉, (29) and with I = 4 and a dependent FMA latency of λ = 6 cycles on Pascal [32] and 4 cycles on Ampere [33], (29) gives Wmin = 8 and 6 warps (Table I). Memory and synchronization latency add further demand, so Wmin is a floor: configurations at or below it should underperform, which the 32×128 row of Table II confirms. The remainder of the section is organized around this model. Section IV-D applies (27), together with the sharedmemory budget, to reduce the configuration space to the two parameters L and N and four candidate configurations, with (29) flagging one of them as hiding-starved. Section IV-E prices each algorithm’s exchange schedule and barrier count in the model’s terms (Table III), yielding testable predictions: where PH and STCR coincide, where their gap opens, and when DTCR’s trade of barriers for warps pays. Section IV-F then confronts these predictions with measurement; no term of (27) is fitted to data—the model supplies signs, floors, and ratios, and the experiments check them. D. Configuration: The Roles of L and N The multi-block algorithms expose two tuning parameters: L threads per group and N samples per thread, with group size C = N L. Their roles on the GPU differ sharply from the

10

TABLE II C ASCADE CONFIGURATIONS AGAINST THE PER -SM BUDGETS OF TABLE I; DEVICE PAIRS READ GTX 1070 / RTX 3060. S HARED MEMORY LIMITS EVERY TB S /SM ENTRY EXCEPT 32×32 ON THE RTX 3060, WHERE THE TB- SLOT CAP LIMITS FIRST. R EGS / THREAD IS THE PTXAS - MEASURED COUNT, SHARED BY PH AND STCR; DTCR USES ABOUT HALF. † BELOW THE HIDING FLOOR OF TABLE I. σsec TBs/SM warps/SM regs/thread PH/STCR/DTCR

L×N

Sbuf

32×32 32×64 32×128 64×64

4.2 kB 23 / 16 8.3 kB 11 / 12 16.5 kB 5 / 6 16.6 kB 5 / 6

23 / 16 11 / 12 5† / 6 10 / 12

72 107 168 105

0/0/2 0/0/2 0/0/2 1/2/2

TABLE III C OMMUNICATION COST PER SECTION PER GROUP ( PER SAMPLE FOR PLR). H ANDOFF : HOW SECTION s’ S OUTPUT REACHES SECTION s+1.

FMAs / sample Exchange rounds Barriers, L≤32 Barriers, 64×64 Section handoff

PH

STCR

DTCR

PLR

∼6 log2 L 0 1 registers

∼8 ∼ 4×2 k log2 B 2 log2 L 2 log2 L split log2 B B 0 2 log2 32 +stg. 2 2 — registers shared mem. none

E. The Four Algorithms

Each algorithm is described by what changes relative to Section II; Table III collects the costs. CPU, where L was the SIMD width (amortizing instructions) PH factorization. PH transfers unchanged: L threads, each and N set the function size. holding N permuted samples. The particular solution (9) is The role of L. Solving a length-L coupled recurrence across per-thread arithmetic with no exchange at all—the O(N ) L threads requires log2 L pairwise exchange rounds at distances dependency chain that limits it on the CPU is hidden by 1, 2, . . . , L/2; this is a lower bound on parallel evaluation [34], residency here. The only inter-thread communication is the not an algorithmic choice. PH traverses this exchange schedule terminal recursive doubling, one traversal of the exchange once per section, cyclic reduction twice. Every round at distance schedule (log2 L rounds, Sklansky-style). At L = 32 every below 32 is an intra-warp shuffle, priced at zero by (27); any round is a shuffle and the barrier count is zero; at L = 64, round at distance 32 or more crosses warps and becomes a exactly one round crosses warps—the last, whose exchange barrier-synchronized exchange. Under SIMT, each sample owns distance of 32 exceeds the warp width, while the earlier rounds its own lane, so the per-sample operation count is invariant at distances 1, 2, . . . , 16 remain in-warp. in L: unlike the CPU, L buys no arithmetic. Below 32, lanes Single-thread-group cyclic reduction (STCR). The lookback idle, so L = 32 is the smallest full-utilization width—all protocol needs both terminal blocks of a group published before rounds shuffles, zero barriers. Each doubling beyond 32 adds the state-dependent remainder, so the GPU implementation one barrier-priced round but enlarges the group, diluting F/C must use the divide-and-conquer formulation of Section II-E. without spending warps. A thread group is the set of L threads that jointly holds one The role of N . Raising N dilutes F/C and increases persignal block group; STCR uses a single thread group, which thread ILP, and never touches the exchange schedule. Its costs: reduces both parities in sequence. Against PH this costs roughly at fixed batch size the group count n/C falls, starving the 8 FMAs per sample instead of 6, and twice the exchange rounds. grid; and the staging buffer (26) grows linearly in N . Shared On the CPU this operation gap decides the ranking; the model memory is the hard limit on the configuration. When it is the prices it at zero wherever the rounds stay in-warp, predicting limiting budget item—and Table II shows it is—the resident that PH and STCR coincide in every L = 32 configuration. TBs per SM are SSM /(4L(N +1)) and each TB holds L/32 Double-thread-group cyclic reduction (DTCR). DTCR has no warps, so the factors of L cancel: CPU counterpart. It assigns the even- and odd-parity reductions to two thread groups of L threads each, running concurrently SSM warps per SM ≈ , (30) within one TB. This doubles the resident warps per group 128 (N + 1) at the same staging footprint—additional hiding supply for a function of N alone. Raising N costs resident warps roughly εhide . The cost is the cascade handoff: the two thread groups one-for-one; raising L merely repackages the same warps into hold disjoint halves of the samples, so every section boundary fewer, larger TBs. requires a shared-memory exchange (two barriers per section, The tested configurations follow from one intrinsic constraint versus the register handoff of PH and STCR). The model and three hardware budgets. The proposed multi-block filtering therefore predicts a conditional outcome: DTCR wins where requires N = kL for integer k. Beyond this, the 48-kB per- hiding is scarce and loses where the added barrier cost exceeds TB shared-memory cap excludes L ≥ 128; warp utilization the hiding it buys, with the boundary set by the device’s cost excludes L < 32; and (30) excludes large N (at 32×256 of barrier β. residency collapses to 2–3 TBs per SM). The four survivors Parallel linear recurrence (PLR). PLR [20] is the highestare 32×32, 32×64, 32×128, and 64×64 (Table II). The last throughput published general linear-recurrence engine and two reach the same C = 4096 by different routes—32×128 serves as the literature baseline, in its native realization: the pays in resident warps, 64×64 pays one cross-warp barrier K sections are combined offline via the z-transform into one round—so their head-to-head is a controlled comparison of the direct-form recurrence of order 2K, computed by hierarchical two currencies. Registers never bind before shared memory: chunk merging with precomputed corrections; the published the cascade footprint is set by the configuration (N resident kernels are transplanted unchanged except for the shared samples plus a fixed overhead) and does not grow with the descriptor flags and the device geometry. The direct form filter order, since the K sections reuse the same registers. pays for its generality twice as the order grows. Its per-thread

11

RTX 3060 (Ampere) throughput (GS/s)

50

Order 2 (K = 1)

Order 4 (K = 2)

Order 8 (K = 4)

Order 16 (K = 8)

roof 45 GS/s

40 30 20 10

GTX 1070 (Pascal) throughput (GS/s)

0 30

roof 32 GS/s

20 10 0

216

218

220

222

224

216

218

220 32×32

222

224

216

218

batch size (samples) 32×64

32×128

220

222

224

216

218

220

222

224

64×64

Fig. 9. PH throughput across the four L×N configurations, against each device’s memory roof (28) (dashed).

carry footprint grows with the order: block sizes must halve SMs at smaller batches. Past saturation the curves rise slowly per order doubling to stay ahead of the register budget, and as the per-launch costs amortize; on the 15-SM GTX 1070 at order 16 no spill-free configuration exists on either device. the large-group configurations also show a local peak at 218 , And its coefficients degrade numerically with the order, leading where the group count fits a single residency wave, followed to numerical failure at order 16 (Section IV-F). by a dip when a second, partial wave appears. At the plateau, the configuration analysis of Section IV-D is visible in full: the order-2 leaders approach the memory roof, and 32×128, F. Measurement and Results the configuration at the hiding floor of Table II, trails on both All experiments use Butterworth low-pass designs of order devices. As the order grows, every curve falls away from the 2K, K ∈ {1, 2, 4, 8}, cutoff 0.2π, designed in second-order roof—the per-section cost grows while the memory floor does sections and b0 -normalized in float32; for PLR the same design not—but at sharply different rates: 64×64, with the largest is combined offline into the equivalent direct form, so the two per-section barrier cost, falls fastest, from 35.5 to 15.2 GS/s on realizations share the transfer function. Batch sizes span 216 to the RTX 3060 between orders 2 and 16, while 32×32 barely 225 samples. Every timed cell runs an accuracy gate: relative moves, from 13.4 to 11.8 GS/s. The latter is dominated by error below 10−4 of the float64 reference peak. The cascade the fixed per-group cost F/C rather than by per-section work, algorithms pass at every order on both devices (errors of order so additional sections are absorbed by stall slots that were 10−7 ). Direct-form PLR passes at orders 2 and 4, is rounding- idle anyway: filter order is nearly free exactly where the fixed limited at order 8 (the float32 direct-form coefficients alone overhead dominates. exceed the gate; its throughput is reported under annotation), Operations are free. PH and STCR carry a factor-of-1.5 and fails entirely at order 16, where rounding pushes poles operation gap into every L = 32 configuration, and the gap outside the unit circle and the output diverges—the original buys nothing: at the plateau, the twelve single-warp cells (three argument for the cascaded realization. Each cell launches its configurations × four orders) coincide inside the ±5% sessionkernel 4000 times as 20 batches of 200 after a 1000-launch noise band (Figure 10), with both signs. On the RTX 3060 warm-up, with a keep-alive kernel holding the boost clock this holds at every order (worst case −4.6%). The GTX 1070 and rotating input buffers defeating the L2 cache; the reported locates the boundary of the claim: the band holds through statistic is the 20th-percentile batch mean, reproducible across order 4, then opens to −9% at order 8 and −13% at order 16. sessions within ±3%. The plateau of a batch curve is the mean The model itself predicts the opening: at order 16 the kernels 24 25 of its 2 and 2 points. perform 12–16 FLOP per byte against the 1070’s balance PH across the configurations. Figure 9 shows PH across the point near 25, so the premise that priced arithmetic at zero is four configurations and all orders. Below saturation, throughput expiring—first on the device with the lower balance point. 16 is governed by the group count n/C: at 2 one of the two small-group configurations leads every cell, and neither C = The barrier price. Where the exchange schedule crosses 4096 configuration ever does, because more groups fill the warps (64×64), STCR’s doubled traversal becomes visible,

12

Fig. 10. (STCR−PH)/PH at the plateau. Gray: the twelve single-warp cells (three L=32 configurations × four orders), inside the ±5% session-noise band—the operation gap has no measurable price. Color: the 64×64 cells, where one exchange round crosses warps; the gap grows with the order and is priced three to five times higher on Pascal. The gray cells leaving the band on the GTX 1070 at orders 8 and 16 mark the balance-point boundary discussed in the text.

and the gap grows with K: from −2.1% at order 2 to −10.6% at order 16 on the RTX 3060, and from −10.2% to −32.6% on the GTX 1070 (Figure 10). The structural gap is the same on both devices; the older architecture prices it three to five times higher. In the terms of (27), σ is the algorithm’s property and β the device’s. DTCR’s conditional trade. The DTCR-versus-PH map at the plateau (Figure 11) follows (27) on both devices. DTCR wins the entire 32×128 column—the hiding-starved configuration— at every order: +7% to +23% on the RTX 3060 and +15% to +30% on the GTX 1070. It never wins where hiding is most ample (32×64 on the RTX 3060, 32×32 on the GTX 1070), because the added warps hide nothing while the two barriers per section still cost. Between these poles the boundary moves with the device, including sign flips (32×32 at order 2: +8% on Ampere, −10% on Pascal). On the RTX 3060 every column degrades toward order 16 as the per-section barrier cost compounds; on the GTX 1070 the hiding-starved 32×128 column instead widens in DTCR’s favor, to +30% at order 16, where PH’s exposed latency deepens faster than DTCR’s barrier bill. Below saturation the trade is no longer conditional (Figure 12): at 216 , with too few resident groups to hide latency, DTCR’s doubled warps supply the missing hiding and it wins every cell on the RTX 3060 (+4% to +60%) and ten of sixteen on the GTX 1070. DTCR is the algorithm of the latency-bound regime; the throughput-bound regime prices its barriers. Two routes to C = 4096, and the best configurations. In the controlled comparison of Section IV-D (Figure 9), the Lroute wins on both devices: 64×64 beats 32×128 by +22% (RTX 3060) and +17% (GTX 1070) for PH at order 2, and still by +7% and +4% for DTCR. One cross-warp barrier round costs less than half the resident warps. Whether the barrier is worth paying at all, however, depends on the device: against the zero-barrier 32×64, 64×64 gains +22% on the RTX 3060 and loses −13% on the GTX 1070. Accordingly, on the GTX 1070 the zero-barrier 32×64 is the best configuration in eleven of twelve algorithm–order cells, while on the RTX 3060 64×64 leads every order-2 and order-4 cell. The best cells are DTCR 64×64 at 38.2 GS/s on the RTX 3060 (85% of its roof) and

Fig. 11. DTCR against PH at the plateau, per order and L×N configuration; blue: DTCR wins. The 32×128 column—the hiding-starved configuration—is a win on both devices, the columns with the most hiding headroom are losses, and cells flipping sign between the two devices are the barrier price β repricing the same trade.

Fig. 12. The same map at the smallest batch (216 ): with too few groups resident to hide latency, DTCR’s doubled warps win almost unconditionally— every cell on the RTX 3060, and the losses on the GTX 1070 confined to the 32×32 column and order 16.

PH 32×64 at 21.1 GS/s on the GTX 1070 (66%). The best configuration is a property of the device, not of the algorithm. PLR. Figure 13 compares the four algorithms at their best configuration per order, as a fraction of each device’s roof. At the order-2 anchor—where cascade and direct form realize the same recurrence and PLR runs exactly as published— PLR reaches 63% of the roof on the RTX 3060 and 56% on the GTX 1070, behind every cascade leader. The deficit widens with order, tracking its k log2 B work law quantitatively: the predicted order-8/order-4 throughput ratio is 0.57 on the RTX 3060 (measured 0.573) and 0.5625 on the GTX 1070 (measured 0.538). At order 16 PLR contributes no valid output at any geometry. The order-2 anchor also bounds the hypothetical cascade of PLR stages without a further measurement: it inherits at best that per-stage throughput and then adds K inter-section handoffs through global memory, so the single-kernel cascade algorithms beat it at every order. Cross-architecture summary. Raw throughput conflates algorithm quality with device hardware: the two devices differ by 1.41× in bandwidth but by 1.7–1.8× in measured leaders. Normalized by each device’s roof, the order-2 leaders read 85% (Ampere) against 66% (Pascal). The dimensionless structure transfers: the PH–STCR coincidence inside the noise band, the growth of the barrier gap with K, the sign map of DTCR’s trade, and PLR’s order ratios all reproduce on both devices. What does not transfer is attributable to three named parameters: the barrier price β falls by a factor of three to five from Pascal to Ampere, moving

13

RTX 3060 (Ampere)

plateau throughput (% of roof)

100

PH STCR DTCR PLR

80 60 40 20

output non-finite

0

GTX 1070 (Pascal)

100

plateau throughput (% of roof)

best configuration per order

38.2 GS/s

marker = crowning config (L×N) 32×32 32×64 32×128 64×64

80 21.1 GS/s

60 40 20 0

output non-finite

2

4

filter order

8

16

Fig. 13. The four algorithms at their best configuration per order (marker: the L×N configuration), as a fraction of each device’s roof (28). The hollow marker at order 16 is zero valid PLR throughput: the direct-form output is non-finite there.

the best configuration from zero-barrier to barrier-tolerant. The hiding budget grows, visible in every roof fraction. And the order-16 boundary of the free-operations claim belongs to the balance point. The dependency graphs of Part I are invariant; the two architectures price them differently, and this section has measured the prices. V. C ONCLUSION This paper extended the theoretical multi-block filtering framework [14] from a single SIMD core to the two architectures where modern throughput resides: multi-core CPUs and GPUs. The obstacle to that extension is the inter-group state dependency—the terminal outputs of each signal block group are the initial conditions of the next—and its resolution is superposition: each group’s output splits into a zero-state response, computable before the group’s state is known, and a homogeneous correction applied when the state arrives. Integrating this decomposition into the multi-block algorithms, and reformulating cyclic reduction in a divide-and-conquer form that exposes both terminal blocks before back substitution, allows groups to execute concurrently on either architecture. The two implementations pair the two deployment scenarios with opposite treatments of the same dependency. For realtime streaming, a wavefront pipeline realized with TBB flow graphs parallelizes across cascade sections, preserving FIFO order and requiring no state protocol at all; a two-level

granularity strategy—function size within a core, grain size across the task system—yields 3.95× scaling on the six performance cores of a heterogeneous Meteor Lake processor, about 2.4 GS/s for a 16th-order filter, and the measurements show that adding efficiency cores degrades the pipeline rather than extending it. For batched processing, a single-kernel GPU implementation carries each group through the entire cascade in registers and parallelizes across groups with a per-section decoupled lookback protocol. A communication-based cost model—memory roof, barrier price, and latency-hiding floor— reduces configuration tuning to the two parameters L and N and predicts the measured behavior: operation counts are free while cross-warp barriers are not, the double-thread-group variant of cyclic reduction wins exactly where latency hiding is scarce, and the best kernels reach 38.2 GS/s on an RTX 3060, 85% of the device’s memory-bandwidth roof; for the 16thorder filter, the best kernel sustains 15.8 GS/s, about 6.5× the multi-core CPU rate. Against the strongest published GPU recurrence engine, evaluated in its native direct form, the cascade algorithms are faster at every filter order and remain numerically valid at order 16, where the direct form produces no usable output. Normalized by each device’s memory roof, the qualitative structure of the results transfers across two GPU generations, with the differences attributable to three named hardware parameters: the barrier price, the latencyhiding budget, and the arithmetic balance point. All GPU kernels presented in this paper are released as an open-source CUDA library, along with the raw measurement data, available at https://github.com/Haotian-RA/cuda recursive filtering. Together with Part I, these results establish cascaded second-order recursive filtering—long considered inherently sequential—as a first-class workload on every tier of commodity parallel hardware, from a single SIMD core to a saturated GPU. Future work will pursue persistent-kernel execution to bring the GPU design into the low-latency streaming regime, extension of the state propagation framework to time-varying and adaptive coefficients, and automated configuration selection driven directly by the cost model. R EFERENCES [1] M. Valkama, M. Renfors, and V. Koivunen, “Advanced methods for I/Q imbalance compensation in communication receivers,” IEEE Transactions on Signal Processing, vol. 49, no. 10, pp. 2335–2344, 2001. [2] J. A. Belloch, B. Bank, L. Savioja, A. Gonzalez, and V. Välimäki, “Multi-channel IIR filtering of audio signals using a GPU,” in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2014, pp. 6692–6696. [3] D. Nehab, A. Maximo, R. S. Lima, and H. Hoppe, “GPU-efficient recursive filtering and summed-area tables,” ACM Trans. Graph., vol. 30, no. 6, p. 1–12, Dec. 2011. [4] P. R. P. Tammana and B. Penumutchi, “An exploration on competent video processing architectures,” in 2023 International Conference for Advancement in Technology (ICONAT), 2023, pp. 1–6. [5] H. Helms, “Fast fourier transform method of computing difference equations and simulating filters,” IEEE Transactions on Audio and Electroacoustics, vol. 15, no. 2, pp. 85–90, 1967. [6] H. Voelcker and E. Hartquist, “Digital filtering via block recursion,” IEEE Transactions on Audio and Electroacoustics, vol. 18, no. 2, pp. 169–176, 1970. [7] C. Burrus, “Block implementation of digital filters,” IEEE Transactions on Circuit Theory, vol. 18, no. 6, pp. 697–701, 1971. [8] ——, “Block realization of digital filters,” IEEE Transactions on Audio and Electroacoustics, vol. 20, no. 4, pp. 230–235, 1972.

14

[9] W. Sung, S. Mitra, and B. Jeren, “Multiprocessor implementation of digital filtering algorithms using a parallel block processing method,” IEEE Transactions on Parallel and Distributed Systems, vol. 3, no. 1, pp. 110–120, 1992. [10] R. Kutil, “Parallelization of IIR filters using SIMD extensions,” in 2008 15th International Conference on Systems, Signals and Image Processing, 2008, pp. 65–68. [11] J. Ahn, H. Chang, J. Cho, and W. Sung, “SIMD processor based implementation of recursive filtering equations,” in 2009 IEEE Workshop on Signal Processing Systems, 2009, pp. 087–092. [12] D.-h. Lee and W. Sung, “Multi-core and SIMD architecture based implementation of recursive digital filtering algorithms,” in 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, 2010, pp. 1550–1553. [13] H. Zhai and B.-P. Paris, “Accurate and efficient implementations of recursive filtering based on SIMD and cascaded form,” in 2024 IEEE 14th Annual Computing and Communication Workshop and Conference (CCWC), 2024, pp. 0001–0007. [14] ——, “Fast cascaded recursive filtering via a block-matrix reformulation,” 2026. [Online]. Available: https://arxiv.org/abs/2607.14054 [15] P. Virtanen et al., “SciPy 1.0: Fundamental algorithms for scientific computing in Python,” Nature Methods, vol. 17, no. 3, pp. 261–272, 2020. [16] V. Volkov, “Understanding latency hiding on GPUs,” Ph.D. dissertation, University of California, Berkeley, 2016. [17] G. E. Blelloch, “Programming parallel algorithms,” Commun. ACM, vol. 39, no. 3, p. 85–97, Mar. 1996. [18] J. JáJá, An Introduction to Parallel Algorithms. Addison-Wesley, 1992. [19] D. Merrill and M. Garland, “Single-pass parallel prefix scan with decoupled look-back,” NVIDIA Corporation, Tech. Rep. NVR-2016002, Mar. 2016. [20] S. Maleki and M. Burtscher, “Automatic hierarchical parallelization of linear recurrences,” in Proceedings of the Twenty-Third International Conference on Architectural Support for Programming Languages and Operating Systems, ser. ASPLOS ’18. Williamsburg, VA, USA: ACM, 2018, pp. 128–138. [21] A. V. Oppenheim and R. W. Schafer, Discrete-time signal processing, 3rd ed. Upper Saddle River, NJ: Prentice Hall, 2010. [22] D.-h. Lee and W. Sung, “GPU based implementation of recursive digital filtering algorithms,” in 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, 2013, pp. 2684–2687. [23] M. Voss, R. Asenjo, and J. Reinders, Pro TBB: C++ Parallel Programming with Threading Building Blocks. New York: APress, 2019. [24] J. Sklansky, “Conditional-sum addition logic,” IRE Trans. Electron. Comput., vol. EC-9, no. 2, pp. 226–231, 1960. [25] W. Gander and G. H. Golub, “Cyclic reduction – history and applications,” in Proceedings of the Workshop on Scientific Computing, 1997, pp. 1–15. [26] C. Lam, “Meteor Lake’s E-Cores: Crestmont Makes Incremental Progress,” May 2024. [Online]. Available: https://chipsandcheese.com/p/ meteor-lakes-e-cores-crestmont-makes-incremental-progress [27] A. Fog, “Vector class library,” https://github.com/vectorclass/version2.git, 2004. [28] G. Ballard, J. Demmel, O. Holtz, and O. Schwartz, “Minimizing communication in numerical linear algebra,” SIAM J. Matrix Anal. Appl., vol. 32, no. 3, pp. 866–901, 2011. [29] L. G. Valiant, “A bridging model for parallel computation,” Commun. ACM, vol. 33, no. 8, pp. 103–111, 1990. [30] S. Williams, A. Waterman, and D. Patterson, “Roofline: An insightful visual performance model for multicore architectures,” Commun. ACM, vol. 52, no. 4, pp. 65–76, 2009. [31] S. Hong and H. Kim, “An analytical model for a GPU architecture with memory-level and thread-level parallelism awareness,” in Proc. ISCA, 2009. [32] Z. Jia, M. Maggioni, J. Smith, and D. P. Scarpazza, “Dissecting the NVidia Turing T4 GPU via microbenchmarking,” arXiv:1903.07486, 2019. [33] H. Abdelkhalik, Y. Arafa, N. Santhi, and A.-H. A. Badawy, “Demystifying the Nvidia Ampere architecture through microbenchmarking and instruction-level analysis,” in Proc. IEEE High Performance Extreme Computing Conf. (HPEC), 2022. [34] L. Hyafil and H. T. Kung, “The complexity of parallel evaluation of linear recurrences,” J. ACM, vol. 24, no. 3, pp. 513–521, 1977.