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] $$
