# How does the Stockham FFT work?

I'm trying to implement a Fourier transform in an OpenGL shader. I found some literature that explains the general principles, but is a bit sparse on some details. And those details are seriously confusing me.

This paper describes the use of the Stockham FFT on the GPU. It contains the following diagram:

I understand how to calculate the discrete fourier transform for any individual point, but I'm confused how exactly the Stockham FFT proceeds. My understanding was that the FFT in general divides the input into odd and even data points, and then recursively processes these. But the diagram I added above looks reversed to me. It starts by dividing the input into the smallest segments and then adds up to larger ones, I would have expected this to work exactly the other way around.

How exactly does the divide-and-conquer part of the Stockham FFT work?

Both algorithms are DIT, but the main difference is in memory access patern.
Cooley-Tukey - the most popular form of transformation from DFT to FFT comes from cache awareness, which is optimal for CPU pipeline.
Cooley-Tukey in first phase reverses bit order while Stockham changes order at each stage.
This is observable if you take look at edges - once you have half of the data it stays in place for Cooley-Tukey but moves for Stockham.

For CPU Stockham makes cache mispredictions while Cooley-Tukey makes thread serialization for GPU.

FFT - look at BFS vs DFS strategy.

FFT stage decomposition - very nice pdf showing butterfly explicitly for different FFT implementations.

There are two variants of FFT: Decimation in Time (DIT) and Decimation in Frequency (DIF). What you describe is DIT, whereas the diagram shows DIF. The only difference is the order in which the bits are processed – LSB to MSB in DIT, and MSB to LSB in DIF.

There are many online resources describing this – see for example these slides.

• Hmm, the paper explicitly mentions DIT in this particular figure Commented Jan 11, 2016 at 21:35
• I guess one of us made a mistake. Commented Jan 11, 2016 at 21:35

The division step involves no computation, it is a mere permutation of the elements. It can indeed be described as a recursive process, top down (consider two halves of the input array and rearrange them, then two halves of each half and son on). In the diagram, this step is only represented by its final result (element-wise permutation).

The conquer step recombines the elements in larger and larger subarrays, and this is explicit in the diagram.

The Stockam method combines subdivision and computation, so that the initial rearragement step is unnecessary.

I could not find a good source on the Stockham FFT either. So here is my attempt of a comprehensive answer. Disclaimer: I did not track down the original source of the Stockham FFT. Most of my information comes from the following non-encrypted japanese website, so make sure You at least use a script blocker if You visit:

http://wwwa.pikara.ne.jp/okojisan/otfft-en/stockham1.html


# Stockham vs. Cooley-Tukey

To unterstand the Stockham FFT, let's compare it to Cooley-Tukey, looking at very simple implementations of both in the Scala programming language, using Spire for the complex numbers. Here's the Stockham FFT:

import spire.implicits.*
import spire.math.Complex
import spire.math.Complex.rootOfUnity

def fft_stockham( u: Array[Complex[Double]] ): Array[Complex[Double]] =

require(Integer.bitCount(u.length) == 1)
val v = u.clone

def fft( off: Int, n: Int ): Unit =

if n <= 1 then return;

val half = n / 2
val  mid = off + half

for i <- 0 until half do
val x = v(off + i)
val y = v(mid + i)
v(off + i) =  x+y
v(mid + i) = (x-y) * Complex.rootOfUnity[Double](n,-i)
end for

fft(off,half)
fft(mid,half)

end fft

fft(0,v.length)
return bitReversalPermutation(v)

end fft_stockham


And here's the Cooley-Tukey FFT:

def fft_cooleyTukey( u: Array[Complex[Double]] ): Array[Complex[Double]] =

require(Integer.bitCount(u.length) == 1)
val v = bitReversalPermutation(u)

def fft( off: Int, n: Int ): Unit =

if n <= 1 then return;

val half = n / 2
val  mid = off + half

fft(off,half)
fft(mid,half)

for i <- 0 until half do
val x = v(off + i)
val y = v(mid + i) * Complex.rootOfUnity[Double](n,-i)
v(off + i) = x+y
v(mid + i) = x-y
end for

end fft

fft(0,v.length)
return v

end fft_cooleyTukey


For completeness' sake, here's a simple implementation of the Bit-reversal permutation:

def bitReversalPermutation( u: Array[Complex[Double]] ): Array[Complex[Double]] =
require( Integer.bitCount(u.length) == 1 )
val s = 1 + Integer.numberOfLeadingZeros(u.length)
return Array.tabulate(u.length)( i => u(Integer.reverse(i<<s)) )
end bitReversalPermutation


Note that there are endless variants of each algorithm:

• stackless implementations
• implementations that perform bit-reversal permutation during recursion
• implementations that use pre-computed roots
• ...

The simplified variants however are particularly well-suited for comparison. You can easily see two characteristic differences between Stockham and Cooley-Tukey:

• Stockham uses a top-down approach, first modifying the array and then calling FFT recursively. Cooley-Tukey uses a bottom-up approach, first applying FFT recursively and then modifying the array.
• Stockham applies bit-reversal permutation at the very end, whereas Cooley-Tukey applies it at the very beginning.

Other than that, You can also see that both implementations are incredibly similar. For that reason, I would tend to disagree with anyone who claims that one algorithm is inherently more efficient than the other. They are just to similar for that to be true. Put the same optimization effort into either one of them and You will very likely end up with the same performance.

# Stockham Use Cases

If both algorithms are so similar, why is the Stockham FFT even a thing? Why not always use Cooley-Tukey? One very good reason is the difference in where the bit-reversal permutation takes place. Let's say You want to compute the positive wrapped convolution (PWC) of two vectors $$u$$ and $$v$$ of the same power-of-2 length. The PWC is an important building block of fast polynomial and large integer multiplication and division. We can use the FFT to compute it:

$$\operatorname{pwc}\left(u,v\right) = \operatorname{ifft}\left( \operatorname{fft}\left(u\right) \circ \operatorname{fft}\left(v\right) \right)$$

Where $$\circ$$ is the element-wise product, which is agnostic to (bit-reversal) permutation in the sense that:

$$\operatorname{bitReverse}\left(\operatorname{bitReverse}\left(u\right) \circ \operatorname{bitReverse}\left(v\right) \right) = u \circ v$$

This means we can use Stockham for the FFT, Cooley-Tukey for the IFFT, drop the bit-reversal permutation from both and still get the correct result for the convolution.

Bit-reversal permutation can be implemented as an $$\mathcal{O}(n)$$ operation. It should be negligible compared to the $$\mathcal{O}\left(n \log{n}\right)$$ operations of the FFT, but due to the pseudo-random nature of bit-reverse indices, it produces many cache-misses, even on modern CPUs. In my experience it can make up 5-10% of the time of a vectorized, optimized FFT implementation, even for large inputs.

# Theory behind Stockham

If You've made it this far, or skipped to this part, You must be really interested in the math behind the Stockham FFT. So here is my best attempt to explain it. Let's say we have a polynomial $$P$$ defined by its coefficient vector $$p = [p_0, p_1, \dots, p_{n-1} ]$$ where $$n$$ is the number of coefficients and also a power of 2.

$$P\left(x\right) = p_0 + p_1 x + p_2 x^2 + \dots + p_{n-1} x^{n-1}$$

Applying the FFT to $$p$$ is equivalent to evaluating $$P$$ at the first $$n$$ powers (starting at zero) of the $$n$$-th root of unity $$\omega_n$$:

$$\operatorname{fft}(p) = \left[P\left(1\right), P\left(\omega_n\right), P\left(\omega_n^2\right), \dots, P\left(\omega_n^{n-1}\right)\right]$$

To split the FFT into two smaller FFT problems, let's look at the even index (starting at 0) and odd index entries of $$\operatorname{fft}(p)$$ separately.

## Even Indices

$$k \in \left\{0, 1, 2, \dots, \frac{n}{2}-1\right\}$$

$$P(\omega_n^{2k}) = p_0 + p_1 \omega_n^{2k} + p_2 \omega_n^{2k \cdot 2} + p_3 \omega_n^{2k \cdot 3} + \dots + \omega_n^{2k \frac{n}{2}} \left( p_{\frac{n}{2}} + p_{\frac{n}{2}+1} \omega_n^{2k} + p_{\frac{n}{2}+2} \omega_n^{2k \cdot 2} + p_{\frac{n}{2}+3} \omega_n^{2k \cdot 3} + \dots \right)$$

Due to the definition of the principal $$n$$-th root of unity we have:

$$\omega_n^{2k \frac{n}{2}} = \left(\omega_n^n\right)^k = 1$$ $$\omega_n^{2k} = \left(\omega_n^2\right)^k = \left(\omega_{\frac{n}{2}}\right)^k$$

Where $$\omega_{\frac{n}{2}}$$ is the $$\frac{n}{2}$$-th root of unity. This allows us to simplify the even index FFT components further:

$$P(\omega_n^{2k}) = p_0 + p_{\frac{n}{2}} + \left(p_1 + p_{\left(\frac{n}{2}+1\right)}\right) \left(\omega_{\frac{n}{2}}\right)^{k} + \left(p_2 + p_{\left(\frac{n}{2}+2\right)}\right) \left(\omega_{\frac{n}{2}}\right)^{k \cdot 2} + \left(p_3 + p_{\left(\frac{n}{2}+3\right)}\right) \left(\omega_{\frac{n}{2}}\right)^{k \cdot 3} + \dots + \left(p_{\left(\frac{n}{2}-1\right)} + p_{n-1}\right) \left(\omega_{\frac{n}{2}}\right)^{k \cdot \left(\frac{n}{2}-1\right)}$$

This means we can use a half sized FFT to compute the even indices of our full-sized FFT as follows:

$$\operatorname{fft}\left(\left[\left(p_0 + p_{\frac{n}{2}}\right), \left(p_1 + p_{\left(\frac{n}{2}+1\right)}\right), \dots, \left(p_{\left(\frac{n}{2}-1\right)} + p_{n-1}\right) \right]\right) = \left[P\left(1\right), P\left(\omega_n^2\right), \dots, P\left(\omega_n^{n-2}\right) \right]$$

## Odd Indices

$$P(\omega_n^{2k + 1}) = p_0 + p_1 \omega_n \omega_n^{2k} + p_2 \omega_n^2 \omega_n^{2k \cdot 2} + p_3 \omega_n^3 \omega_n^{2k \cdot 3} + \dots + \left(\omega_n\right)^{\frac{n}{2}} \omega_n^{2k\frac{n}{2}}\left(p_{\frac{n}{2}} + p_{\left(\frac{n}{2}+1\right)} \omega_n \omega_n^{2k} + p_{\left(\frac{n}{2}+2\right)} \omega_n^2 \omega_n^{2k \cdot 2} + p_{\left(\frac{n}{2}+3\right)} \omega_n^3 \omega_n^{2k \cdot 3} + \dots\right)$$

From the properties of principal roots of unity, we can derive:

$$\left(\omega_n\right)^{\frac{n}{2}} = -1$$

Analogous to the proof of the even index entries, we get the following half-sized FFT for the odd index entries:

$$\operatorname{fft}\left(\left[\left(p_0 - p_{\frac{n}{2}}\right), \left(p_1 - p_{\left(\frac{n}{2}+1\right)}\right) \cdot \omega_n, \left(p_2 - p_{\left(\frac{n}{2}+2\right)}\right) \cdot \omega_n^2, \dots, \left(p_{\left(\frac{n}{2}-1\right)} - p_{n-1}\right) \cdot \omega_n^{\left(\frac{n}{2}-1\right)} \right]\right) = \left[P\left(\omega_n\right), P\left(\omega_n^3\right), \dots, P\left(\omega_n^{n-1}\right) \right]$$

• Damn you, Dyslexia! Apologies for that! Commented Aug 12, 2023 at 8:51