Stockham FFT Butterfly Algorithm in Parallel GLSL

This article evaluates the implementation and performance of the Stockham butterfly algorithm for parallel Fast Fourier Transforms (FFT) within the OpenGL Shading Language (GLSL). While classical FFT algorithms like Cooley-Tukey remain standard on single-threaded CPUs, the Stockham formulation is uniquely optimized for massive parallelization on GPUs. By removing the expensive global bit-reversal permutation step and naturally facilitating out-of-place ping-pong rendering, the Stockham algorithm achieves superior data locality, uniform execution branching, and high computational throughput across shader execution pipelines.

The Bit-Reversal Bottleneck and the Stockham Advantage

The standard Cooley-Tukey decimation-in-time or decimation-in-frequency algorithms require reordering the data array using a bit-reversal permutation, either before the computation stages begin or after they conclude. On massively parallel architectures such as GPUs, this bit-reversal operation introduces severe memory-access bottlenecks. It requires scattered global reads and writes, divergent memory access patterns, and execution stalls that undermine parallel hardware performance.

The Stockham formulation resolves this issue by performing an automatic, step-by-step reordering across successive butterfly passes. Because the algorithm is inherently out-of-place, data read from a source array is combined and written directly into its correct sequential order in a destination array at each stage. In GLSL, this eliminates the need to execute a distinct bit-reversal shader pass.

Execution via GLSL Compute Shaders and Textures

In GLSL, the Stockham algorithm is evaluated either via compute shaders using uniform/shared memory or via fragment shaders using texture ping-ponging:

  1. Ping-Pong Buffer Strategy: Two Framebuffer Objects (FBOs) or texture storage targets alternate between acting as the read source and the write destination. For an FFT of size \(N = 2^M\), exactly \(M\) sequential passes are dispatched.
  2. Index Computation: In each stage \(s \in [0, M-1]\), every shader invocation maps directly to a target coordinate. A fragment shader computes its unique output index \(x\). The corresponding pair of inputs are retrieved from memory at indices: \[j = \lfloor x / 2^{s+1} \rfloor \cdot 2^s + (x \pmod{2^s})\] \[k = j + N / 2\]
  3. Twiddle Factor Calculation: The twiddle factor \(W_N^r = e^{-i 2\pi r / N}\) is computed dynamically in the shader using built-in GLSL trigonometric functions (cos and sin), where \(r = (x \pmod{2^{s+1}}) \cdot (N / 2^{s+1})\). Evaluating twiddle factors dynamically avoids the memory overhead of reading from precomputed lookup textures, effectively trading negligible ALU load for critical memory bandwidth savings.

Multi-Dimensional Transforms and Separability

For 2D or 3D datasets—such as spectral ocean wave simulations or image processing—GLSL evaluates the Stockham FFT separably:

  • Horizontal Pass: The 1D Stockham transform is executed along the rows of the texture grid across \(\log_2(Width)\) stages.
  • Vertical Pass: The same process is repeated along the columns across \(\log_2(Height)\) stages.

Because texture memory is natively 2D, the vertical pass incurs no structural indexing overhead compared to the horizontal pass. GPUs benefit from hardware texture caches that accelerate localized spatial lookups during these transpositions.

Performance and Pipeline Evaluation

When evaluating the Stockham butterfly in GLSL, execution performance is governed by memory throughput rather than arithmetic capability:

  • Branching Consistency: Every parallel thread in a GLSL invocation executes the exact same mathematical operations per stage, preventing SIMD warp or wave divergence.
  • Vector Optimization: GLSL naturally processes complex numbers using vec2 types (where x represents real and y represents imaginary components). Modern graphics hardware performs complex addition and multiplication in minimal clock cycles using fused multiply-add (FMA) instructions.
  • Cache Efficiency: When using compute shaders with shared memory (local workgroup memory), batches of butterfly steps can be completed within fast on-chip registers before writing results back to global texture memory, substantially reducing overall bus traffic.

The Stockham butterfly algorithm provides an ideal balance of linear memory access patterns and predictable arithmetic steps, making it one of the most efficient choices for real-time Fourier analysis in GLSL environments.