a fast CPU FFT

8 Sep 2026

The fast Fourier transform (FFT) is one of the most important algorithms in computer science. Modern telephony, image compression, signal analysis, and many other applications rely on the FFT at their core. It’s correspondingly well researched, and has been implemented many times over. I will be reimplementing it myself, (1) because it’s fun, and (2) because I’ll need an unusual version of it for my ocean water simulation. This document serves as a follow-along derivation of my optimized CPU implementation, which exceeds the performance of rustfft’s scalar code at the input sizes I care about.

Background

The FFT is just an efficient way of computing something called the “discrete Fourier transform.” What is that, and why do we care? Essentially, a Fourier transform decomposes an input “signal”, such as a sound wave or an image, into a bunch of sines and cosines. When you add those sines and cosines together, you get back the original signal. In the real world, signals are typically considered to be “continuous”, meaning they aren’t composed of blocks of a minimum size (if you ignore quantum mechanics). The classical Fourier transform deals with those kinds of signals. However, our digital lives are composed of blocks of minimum size: 1s and 0s, or bits. The discrete Fourier transform (DFT) deals with this kind of signal.

This representation is extremely useful in a number of applications. For example, to compress an image, you can simply strip away all the high-frequency components of the signal. It turns out the human eye can’t really tell, and you can make an image much, much smaller before it becomes obvious that it’s been compressed. The same goes for music, telephony, and video. I will be using it to implement a realtime ocean water simulation - it turns out that the sum-of-sines representation is actually a highly accurate way to represent the dynamics of so-called “fully developed” oceans. More on that in a followup article - for now, we focus on the FFT.

The DFT

The DFT is defined as follows¹

$$X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}$$

Let’s unpack that.

Our input signal $x$ is comprised of $N$ samples: $x = [x_0, x_1, ... x_{N-1}]$. A typical 1-second sound wave would have 44,100 samples, where each sample represents the air pressure that a microphone measured at that point in time. For that reason, we say that the input signal is in the time domain.

The DFT version of our signal, $X$, is also comprised of $N$ samples, but we index them with $k$ instead: $X_k = [X_0, X_1, ... X_{N-1}]$.

We can simplify our expression a bit to get a sense of what’s happening:

$$X_k = \sum_{n=0}^{N-1} x_n W^{nk}_N$$

To get the $k$th term of the DFT, we have to add up every component of the input signal multiplied by some term $W^{nk}_N$. Since there are $N$ terms in the DFT, we must do (N additions of the input signal) * (N times for the output signal). Thus this is an $O(N^2)$ algorithm. We will get back to this later!

The inner term, $e^{-i ...}$ might be a head scratcher. What does it mean to exponentiate by an imaginary number? Where are the sines and cosines? Well, there’s a famous formula² from calculus which tells us that:

$$e^{it} = \cos{t} + i \sin{t}$$

(This formula drops out of the Maclaurin series for $e^x$, $\cos x$, and $\sin x$. By rearranging terms you wind up with this identity.)

So although our expression is expressed in the form $e^{i \dots}$, it is really representing a sum of sines and cosines. Neat!

Finally, there is something unintuitive to reflect on. In the real numbers, there are at most two solutions to this equation, assuming that $k$ is a natural number (1, 2, 3…):

$$1 = x^k$$

If $k$ is even, the only solution is $x = 1$; if $k$ is even, there is also the solution $x = -1$.

In the complex numbers, we can have more than one solution. In general, for any natural number $k$, there are $k$ solutions, and they are of the form:

$$1 = e^{\frac{2 \pi p}{k} i}$$

… where $p \in [1, k]$

(Read as “$p$ is an element of the range of numbers starting at 1 and ending at $k$”).

For $k=1$:

$$e^{\frac{2 \pi i}{1}} = \cos{2\pi} + i \sin{2\pi} = 1 + 0 = 1$$

For $k=2$:

$$\begin{array}{rclclcl} e^{\frac{2 \pi i}{2} 1} &=& \cos{\pi} + i \sin{\pi} &=& -1 + 0 &=& -1 \\ e^{\frac{2 \pi i}{2} 2} &=& \cos{2\pi} + i \sin{2\pi} &=& 1 + 0 &=& 1 \end{array}$$

So for a given $k$, the set of complex numbers that satisfy the relationship $1 = x^k$ are called the $k$th roots of unity, and there are $k$ of them. (“Unity” is just another word for 1.) We can visualize them as simply dividing a circle in the complex plane:

In summary:

• $x_n$ is the $n$th term of the input signal. There are $N$ total input terms.
• $X_k$ is the $k$th term of the DFT. There are $N$ total terms.
• $W^{nk}_N = e^{-\frac{2 \pi i}{N} nk}$. We observe that this is an $N$th root of unity.
• Each term of the DFT multiplies each term of the input by $W^{nk}_N$.
• Naively evaluating a DFT takes O(N^2) time.

Implementing the DFT

The DFT is fairly straightforward to implement in code. Here it is in Rust:

# Converts a `usize` to a float.
fn usize_to_float<T: Float>(value: usize) -> T {
    num::cast(value).unwrap()
}

# Evaluates the DFT of `data`.
fn naive_dft<T: Float + FloatConst>(data: &mut [Complex<T>]) {
    let big_n = data.len();
    let mut result = vec![Complex::new(T::zero(), T::zero()); big_n];
    for k in 0..big_n {
        for n in 0..big_n {
            let k_t     = usize_to_float::<T>(k);
            let n_t     = usize_to_float::<T>(n);
            let big_n_t = usize_to_float::<T>(big_n);
            let phase = -T::TAU() * k_t * n_t / big_n_t;
            let factor = Complex::<T>::cis(phase);
            result[k] = result[k] + data[n] * factor;
        }
    }
    data.copy_from_slice(&result);
}

This is technically correct, but there are many problems with this code:

1. The factors $W_N^{nk}$ are recomputed each time we call this function, even though they do not change with respect to data. We should hoist that computation out.
2. $W_N^{nk}$, also called twiddles, are computed with type T, which may be a low-precision float. We should compute them in high precision, then cast to low-precision at the end. Hoisting them out of this function also justifies running that computation in high precision, since it’s no longer on the hot path.
3. We accumulate floating point adds sequentially, which accumulates more error than if we accumulated them via a binary tree.
4. The phase calculation does several floating point multiplications and divisions in the hot path. Had we hoisted our twiddles out, we could get away with no divisons and a single multiply. More on that later.
5. The copy at the end is expensive, and we’d like to avoid it if possible.
6. Converting floats to ints in the hot path is not free.
7. We allocate and initialize an array, result, on the hot path. It’s better than doing it on the heap, but it’s still slow. The allocation should be hoisted out.

We won’t be addressing those until we get into our fast Fourier transform, but I want to start pointing out the kinds of issues we need to think about. The name of the game is doing as little work as possible in the hot path.

DFT Evaluation

Let’s take a look at the DFT’s numerical accuracy and speed. We will be comparing against rust’s rustfft crate as our speed of light. We will also be using a 4096-element array of randomized elements to measure both our numeric accuracy and speed. When measuring performance, we use Criterion to minimize the effects of cache hotness, scheduling noise, etc. To measure error, we use rustfft on a 64-bit signal as our source of truth. Finally, we will disable all vectorization (AVX/SSE) when measuring performance, since our end goal is a GPU-friendly algorithm which won’t have access to those intrinsics.

The results are as follows:

(Input size 4096, type f32.)

The speed-of-light implementation is not only ~5,690x faster, it’s ~2,190x more accurate in the worst case, and ~2,353x more accurate on average.

So, how are we going to bridge this gap?

The fast Fourier transform

As highlighted above, naively evaluating a DFT takes $O(N^2)$ time, where $N$ is the length of the input signal. There is an algorithm appropriately named the fast Fourier transform (FFT) which evaluates the same result in $O(N \log N)$ time. It works by dividing the input into two parts, evaluating the FFT on each part (which is now half as big), then using some clever math to efficiently combine the results. Let’s get into it.

Recall the definition of the DFT:

$$X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n}$$

We can split this by even and odd indices $n$:

$$X_k = \sum_{n=0}^{N/2-1} x_{2n} e^{-\frac{2 \pi i}{N} k (2n)} + \sum_{n=0}^{N/2-1} x_{2n+1} e^{-\frac{2 \pi i}{N} k (2n+1)}$$

Next, factor out $e^{-\frac{2\pi i}{N}k}$ from the second sum:

$$X_k = \sum_{n=0}^{N/2-1} x_{2n} e^{-\frac{2 \pi i}{N} k 2n} + e^{-\frac{2\pi i}{N}k} \sum_{n=0}^{N/2-1} x_{2n+1} e^{-\frac{2 \pi i}{N} k 2n}$$

(This factoring follows from the fact that, in general, $a^{b+1} = a a^b$.)

Inside the sum, multiply the exponent by $\frac{1/2}{1/2}$, i.e. 1:

$$\begin{align*} X_k &= \sum_{n=0}^{N/2-1} x_{2n} e^{-\frac{2 \pi i}{N/2} k n} + e^{-\frac{2\pi i}{N}k} \sum_{n=0}^{N/2-1} x_{2n+1} e^{-\frac{2 \pi i}{N/2} k n} \\ &= E_k + e^{-\frac{2\pi i}{N}k} O_k \end{align*}$$

Note what just happened: we have represented the $k$th term of the DFT in terms of the sums of two DFT’s with half as many terms! That is the essence of how the FFT runs in $O(N \log N)$ time. The only lurking issue is that this only holds for $k$ in the range $[0, N/2)$. To get $k$ in the range $[N/2, N)$, we have to do some analysis. We will replace every instance of $k$ with $k+N/2$, then attempt to refactor the expression to get a result that only deals with indices of $k$:

In conclusion:

$$\begin{align*} X_k &= E_k + W_N^K O_k \\ X_{k+N/2} &= E_k - W_N^K O_k \end{align}$$

Let’s reflect on a couple things.

First, we divide the input into evens and odds. This only works if the input is divisible by 2. Since we’re going to be doing this recursively, we actually need it to be a power of 2. We can relax this by dividing the input into thirds, fourths, fifths, etc., which we’ll have to get into later. If at all possible, you should try to FFT an input signal with a length whose prime factors are small. This lets us apply various analytic tricks to make it fast. It’s common to pad with 0s, although that can create artifacts in the frequency-domain spectrum.

Second, splitting the input into even and odd terms isn’t the only choice. This approach is called decimation in time, because you still have samples near the beginning and end, but half as many overall. Your sample rate has halved, but the time interval is about the same. We might instead split it into a lower and upper half. This approach is called decimation in frequency: your time intervals halve, but the frequency rate in each half is the same.

Implementing the FFT

The FFT is far less trivial to implement than the DFT. Here is the simplest code I could come up with:

// Checks that `n = k^p`, for some natural number `p`.
fn is_power_of_k(n: usize, k: usize) -> bool {
    match n {
        0 => false,
        1 => true,
        _ => n % k == 0 && is_power_of_k(n / k, k),
    }
}

// Helper to naive_fft. Takes `data` along with 3 numbers that let us recreate an even-odd subset:
//  - `start_idx` tells us where the subset begins;
//  - `big_n` is the number of elements in the subset;
//  - `stride` is the distance between elements.
// We also use a double buffer, `scratch`, to avoid clobbering data while merging results.
#[rustfmt::skip]
fn _naive_fft<T: Float + FloatConst>(data: &mut [Complex<T>], start_idx: usize, big_n: usize, stride: usize, scratch: &mut [Complex<T>]) {
    if big_n == 1 {
        return;
    }
    // Compute DFT of even elements.
    _naive_fft(data, start_idx,        big_n/2, stride*2, scratch);
    // Odd elements.
    _naive_fft(data, start_idx+stride, big_n/2, stride*2, scratch);
    for k in 0..(big_n/2) {
        let p = data[start_idx + 2 * k       * stride];
        let q = data[start_idx + (2 * k + 1) * stride];
        let k_t     = usize_to_float::<T>(k);
        let big_n_t = usize_to_float::<T>(big_n);
        let phase = -T::TAU() * k_t / big_n_t;
        let factor = Complex::<T>::cis(phase);
        scratch[start_idx + k               * stride] = p + q * factor;
        scratch[start_idx + (k + big_n / 2) * stride] = p - q * factor;
    }
    data.copy_from_slice(scratch);
}

// Naive implementation of Cooley-Tukey FFT. Modifies `data`in place. Panics if data.len() is not a power of two.
#[allow(dead_code)]
fn naive_fft<T: Float + FloatConst>(data: &mut [Complex<T>]) {
    assert!(is_power_of_k(data.len(), 2));
    let mut scratch = Vec::from(data.as_ref());
    _naive_fft(data, 0, data.len(), 1, &mut scratch);
}

This code is obviously highly suboptimal, for many of the same reasons as the DFT code. In addition, we also copy the entire array once per recursive call. There are $O(N)$ recursive calls, so this is extremely wasteful. We’ll fix that later by double-buffering.

Inefficiencies aside, this code still performs vastly better than the naive DFT:

(Input size 4096, type f32.)

We get a nice 62x speedup, and 1335x improvement on average error. However, the speed-of-light implementation is still ~91x faster than ours, and 1.7x more accurate. Most of our work is going to focus on bridging these two gaps, while retaining as simple an implementation as possible.

Note for a moment the impact of sequential adds. The naive DFT performed 4096 sequential adds for each term, and wound up ~2000x less accurate than the speed-of-light. Due to the FFT’s recursive structure, we use log(4096) = 12 sequential adds, and that brings our accuracy within a factor of 2 of optimal. Quite the stark difference!

Opt. 1: Precompute twiddles

The twiddle factors, $W_N^{nk}$, do not depend on the input to the FFT, so we can (and should!) hoist them out of the hot path. Production FFT libraries like fftw and rustfft do this, and we’ll follow in their footsteps. This also lets us precompute the twiddles in high precision before casting to low precision, which as we’ll see, improves the precision of the end result.

First, let’s precompute our twiddle factors:

// Calculates the "twiddle factors" for an n-element FFT, aka all of the nth roots of unity.
fn precompute_twiddles<T: Float + FloatConst>(n: usize) -> Vec<Complex<T>> {
    let mut result = vec![Complex::<T>::new(T::zero(), T::zero()); n];

    let n_f64 = usize_to_float::<f64>(n);
    for i in 0..n {
        let tw_f64 = Complex::<f64>::cis(-f64::TAU() * usize_to_float::<f64>(i) / (n_f64));
        result[i] = Complex::new(T::from(tw_f64.re).unwrap(), T::from(tw_f64.im).unwrap());
    }

    result
}

Next, adjust our function to take these twiddles as input:

fn _fft_v1_hoist<T: Float + FloatConst>(
    data: &mut [Complex<T>],
    start_idx: usize,
    big_n: usize,
    stride: usize,
    scratch: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    if big_n == 1 {
        return;
    }
    // Compute DFT of even elements.
    _fft_v1_hoist(data, start_idx, big_n / 2, stride * 2, scratch, twiddles);
    // Odd elements.
    _fft_v1_hoist(
        data,
        start_idx + stride,
        big_n / 2,
        stride * 2,
        scratch,
        twiddles,
    );
    for k in 0..(big_n / 2) {
        let p = data[start_idx + 2 * k * stride];
        let q = data[start_idx + (2 * k + 1) * stride];
        let factor = twiddles[k * stride];
        scratch[start_idx + k * stride] = p + q * factor;
        scratch[start_idx + (k + big_n / 2) * stride] = p - q * factor;
    }
    data.copy_from_slice(scratch);
}

// Modification of fft_naive: hoist out and precompute twiddles.
pub fn fft_v1_hoist<T: Float + FloatConst>(data: &mut [Complex<T>], twiddles: &[Complex<T>]) {
    assert!(is_power_of_k(data.len(), 2));
    let mut scratch = Vec::from(data.as_ref());
    _fft_v1_hoist(data, 0, data.len(), 1, &mut scratch, &twiddles);
}

We see a modest performance uplift, but our average-case error is now within spitting distance of the speed-of-light, and our worst-case error matches exactly:

Opt. 2: Double buffering

Our FFT algorithms thus far have done a fully length-$N$ copy at each recurisive step. Because each recursive step divides the length of the array by 2, we make a total of $1 + 2 + 4 + \dots N/2$ function calls, which sums to $N-1$ total calls. Each one does a copy of length $N$, so if each copy takes $O(N)$ itme, we spend $O(N^2)$ time copying buffers overall. Not good!

We can fix this pretty easily with double buffering:

fn _fft_v2_double_buffer<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    start_idx: usize,
    big_n: usize,
    stride: usize,
    twiddles: &[Complex<T>],
) {
    if big_n == 1 {
        return;
    }
    // Compute DFT of even elements.
    _fft_v2_double_buffer(dst, src, start_idx, big_n / 2, stride * 2, twiddles);
    // Odd elements.
    _fft_v2_double_buffer(
        dst,
        src,
        start_idx + stride,
        big_n / 2,
        stride * 2,
        twiddles,
    );
    for k in 0..(big_n / 2) {
        let p = src[start_idx + 2 * k * stride];
        let q = src[start_idx + (2 * k + 1) * stride];
        let factor = twiddles[k * stride];
        dst[start_idx + k * stride] = p + q * factor;
        dst[start_idx + (k + big_n / 2) * stride] = p - q * factor;
    }
}

#[allow(dead_code)]
pub fn fft_v2_double_buffer<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    assert!(is_power_of_k(src.len(), 2));
    dst.copy_from_slice(src);
    // Switching `src` and `dst` means that at the end, the result is in `src` - which is actually
    // what we want! We will be hiding `dst` and `twiddles` in a struct later on :)
    _fft_v2_double_buffer(dst, src, 0, src.len(), 1, twiddles);
}

Note that we only hoist out the allocation of the double-buffer. Initialization still occurs in the hot path.

Accuracy numbers are identical to before, as expected, and performance is vastly improved:

Pretty remarkable result. Minimizing memory writes gets us within a factor of 3 of the state of the art.

Still… we can go faster!

Opt. 3: Iterative instead of recursive

The recursive implementation we’re using is good for the classroom, but bad for performance. If we switch to an iterative implementation, we’ll be able to share work each time we step down a layer of recursion. It will also make it much easier to map this algorithm to the GPU (more on that later).

Let’s do it:

pub fn fft_v3_iterative<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    assert!(is_power_of_k(src.len(), 2));
    dst.copy_from_slice(src);
    let n_iter = log_k_of::<2>(src.len());

    if n_iter % 2 != 0 {
        dst.copy_from_slice(src);
    }

    let (mut input, mut output) = if n_iter % 2 == 0 {
        (dst, src)
    } else {
        (src, dst)
    };
    let mut stride = input.len();
    let mut big_n = 1;
    for _ in 0..n_iter {
        stride /= 2;
        big_n *= 2;
        std::mem::swap(&mut input, &mut output);

        for start_idx in 0..stride {
            for k in 0..big_n / 2 {
                // Get odd and even elements.
                let p = input[start_idx + 2 * k * stride];
                let q = input[start_idx + (2 * k + 1) * stride];
                // Combine.
                let factor = twiddles[k * stride];
                output[start_idx + k * stride] = p + q * factor;
                output[start_idx + (k + big_n / 2) * stride] = p - q * factor;
            }
        }
    }
}

This is essentially identical to the v2 code, except that we use iteration instead of recursion. Regardless, the performance uplift is dramatic:

We’re well within a factor of 2 of SOTA now! No, we’re not done.

Aside: the radix-4 FFT

Let’s think, for a moment, what our FFT would look like if instead of splitting the input into 2 parts at each stage, we broke it into 4:

$$\begin{array}{rcll} X_k = & \sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & \sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n+1)k} & + \\ & \sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n+2)k} & + \\ & \sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n+3)k} & \end{array}$$

Apply the usual factoring trick:

$$\begin{array}{rclll} X_k = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & e^{-\frac{2\pi i}{N}k} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & e^{-\frac{2\pi i}{N}2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & e^{-\frac{2\pi i}{N}3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)k} & \end{array}$$

This is only valid for $k$ on $[0, N/4)$. To get the others we have to do the same analysis as before - replace every $k$ with $k + N/4$, then do some eliminations and factoring.

Calculating $k+N/4$

$$\begin{array}{rclll} X_{k+N/4} = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)(k+N/4)} & + \\ & e^{-\frac{2\pi i}{N}(k+N/4)} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)(k+N/4)} & + \\ & e^{-\frac{2\pi i}{N}2(k+N/4)} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)(k+N/4)} & + \\ & e^{-\frac{2\pi i}{N}3(k+N/4)} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)(k+N/4)} & \end{array}$$

Simplify the shared inner term:

$$\begin{align*} e^{-\frac{2\pi i}{N}(4n)(k+N/4)} &= e^{-\frac{2\pi i}{N}(4n)k} e^{-\frac{2\pi i}{N}(4n)N/4} \\ &= e^{-\frac{2\pi i}{N}(4n)k} e^{-2\pi i n} \\ &= e^{-\frac{2\pi i}{N}(4n)k} \end{align*}$$

Simplify the first outer term:

$$\begin{align*} e^{-\frac{2\pi i}{N}(k+N/4)} &= e^{-\frac{2\pi i}{N}k} e^{-\frac{2\pi i}{N}N/4} \\ &= e^{-\frac{2\pi i}{N}k} e^{-\frac{2\pi i}{4}} \\ &= e^{-\frac{2\pi i}{N}k} (-i) \\ \end{align*}$$

By inspection, we can see that the second and third terms will be of this form as well. We’re basically just multiplying by a vector that’s rotating 90 degrees clockwise in the complex plane:

$$\begin{align*} e^{-\frac{2\pi i}{N}(2k+2N/4)} &= e^{-\frac{2\pi i}{N}2k} e^{-\frac{2\pi i 2}{4}} \\ &= e^{-\frac{2\pi i}{N}2k} (-1) \\ e^{-\frac{2\pi i}{N}(3k+3N/4)} &= e^{-\frac{2\pi i}{N}3k} e^{-\frac{2\pi i 3}{4}} \\ &= e^{-\frac{2\pi i}{N}3k} (i) \\ \end{align*}$$

Plugging in:

$$\begin{array}{rrlll} X_{k+N/4} = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-i) e^{-\frac{2\pi i}{N}k} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-1) e^{-\frac{2\pi i}{N}2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (i) e^{-\frac{2\pi i}{N}3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)k} & \end{array}$$

Using $W$ syntax:

$$\begin{array}{rrlll} X_{k+N/4} = & &\sum_{n=0}^{N/4-1} x_{4n} & W_N^{4nk} & + \\ & (-i) W_N^k &\sum_{n=0}^{N/4-1} x_{4n+1} & W_N^{4nk} & + \\ & (-1) W_N^{2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & W_N^{4nk} & + \\ & (i) W_N^{3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & W_N^{4nk} & \end{array}$$

Calculating $k+N/2$

$$\begin{array}{rclll} X_{k+N/2} = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)(k+N/2)} & + \\ & e^{-\frac{2\pi i}{N}(k+N/2)} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)(k+N/2)} & + \\ & e^{-\frac{2\pi i}{N}2(k+N/2)} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)(k+N/2)} & + \\ & e^{-\frac{2\pi i}{N}3(k+N/2)} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)(k+N/2)} & \end{array}$$

Simplify the shared inner term:

$$\begin{align*} e^{-\frac{2\pi i}{N}(4n)(k+N/2)} &= e^{-\frac{2\pi i}{N}(4n)k} e^{-\frac{2\pi i}{N}(4n)N/2} \\ &= e^{-\frac{2\pi i}{N}(4n)k} e^{-2\pi i 2n} \\ &= e^{-\frac{2\pi i}{N}(4n)k} \end{align*}$$

(We can see from the above that the last quarter will also have the same simplification applied, so we will skip deriving it later.)

Simplify the first outer term:

$$\begin{align*} e^{-\frac{2\pi i}{N}(k+N/2)} &= e^{-\frac{2\pi i}{N}k} e^{-\frac{2\pi i}{N}N/2} \\ &= e^{-\frac{2\pi i}{N}k} e^{-\frac{2\pi i}{2}} \\ &= e^{-\frac{2\pi i}{N}k} (-1) \\ \end{align*}$$

Let’s pause here to reflect. In the $[0, N/4)$, we rotated our outer terms by a quarter turn in the complex plane for each term. Now we’re rotating by a half turn. The next leg, we will rotate by 3/4 of a turn.

I will truncate the derivation there. The reader may do the rest as an exercise if needed.

Plugging in:

$$\begin{array}{rrlll} X_{k+N/2} = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-1) e^{-\frac{2\pi i}{N}k} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (+1) e^{-\frac{2\pi i}{N}2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-1) e^{-\frac{2\pi i}{N}3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)k} & \end{array}$$

Using $W$ syntax:

$$\begin{array}{rrlll} X_{k+N/2} = & &\sum_{n=0}^{N/4-1} x_{4n} & W_N^{4nk} & + \\ & (-1) W_N^k &\sum_{n=0}^{N/4-1} x_{4n+1} & W_N^{4nk} & + \\ & (+1) W_N^{2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & W_N^{4nk} & + \\ & (-1) W_N^{3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & W_N^{4nk} & \end{array}$$

Calculating $k+3N/4$

Per the lemmas in the last section, we can jump right to the result:

$$\begin{array}{rrlll} X_{k+3N/4} = & &\sum_{n=0}^{N/4-1} x_{4n} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (+i) e^{-\frac{2\pi i}{N}k} &\sum_{n=0}^{N/4-1} x_{4n+1} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-1) e^{-\frac{2\pi i}{N}2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & e^{-\frac{2\pi i}{N}(4n)k} & + \\ & (-i) e^{-\frac{2\pi i}{N}3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & e^{-\frac{2\pi i}{N}(4n)k} & \end{array}$$

Using $W$ syntax:

$$\begin{array}{rrlll} X_{k+3N/4} = & &\sum_{n=0}^{N/4-1} x_{4n} & W_N^{4nk} & + \\ & (+i) W_N^k &\sum_{n=0}^{N/4-1} x_{4n+1} & W_N^{4nk} & + \\ & (-1) W_N^{2k} &\sum_{n=0}^{N/4-1} x_{4n+2} & W_N^{4nk} & + \\ & (-i) W_N^{3k} &\sum_{n=0}^{N/4-1} x_{4n+3} & W_N^{4nk} & \end{array}$$

Summary

The radix-4 FFT’s merge step works as follows:

And for each stage, the twiddles are:

• Term 0: 1
• Term 1: $W_N^k$
• Term 2: $W_N^{2k}$
• Term 3: $W_N^{3k}$

Opt. 4: Radix-4

With the above in mind, we can now implement the radix-4 FFT:

#[inline(always)]
fn mul_ni<T: Float + FloatConst>(x: Complex<T>) -> Complex<T> {
    Complex::new(x.im, -x.re)
}

pub fn fft_v4_radix_4<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    assert!(is_power_of_k(src.len(), 4));
    let n_iter = log_k_of::<4>(src.len());

    dst.copy_from_slice(src);

    let (mut input, mut output) = if n_iter % 2 == 0 {
        (dst, src)
    } else {
        (src, dst)
    };
    let big_n = input.len();
    let mut stride = big_n;
    let mut big_n = 1;
    for _ in 0..n_iter {
        stride /= 4;
        big_n *= 4;
        std::mem::swap(&mut input, &mut output);

        for start_idx in 0..stride {
            for k in 0..big_n / 4 {
                // Collect inputs.
                let i0 = input[start_idx + 4 * k * stride];
                let i1 = input[start_idx + (4 * k + 1) * stride];
                let i2 = input[start_idx + (4 * k + 2) * stride];
                let i3 = input[start_idx + (4 * k + 3) * stride];
                // Collect relevant twiddles.
                let ot1 = twiddles[1 * k * stride];
                let ot2 = twiddles[2 * k * stride];
                let ot3 = twiddles[3 * k * stride];

                let a = i0;
                let b = ot1 * i1;
                let c = ot2 * i2;
                let d = ot3 * i3;

                // To derive this, write the expression below in terms of
                // a/b/c/d, then factor out!
                let ac_sum  = a + c;
                let ac_diff = a - c;
                let bd_sum = b + d;
                let bd_diff_ni = mul_ni(b - d);

                output[start_idx + k * stride] = ac_sum + bd_sum;
                output[start_idx + (k + big_n / 4) * stride] =
                    ac_diff + bd_diff_ni;
                output[start_idx + (k + big_n / 2) * stride] =
                    ac_sum - bd_sum;
                output[start_idx + (k + 3 * big_n / 4) * stride] =
                    ac_diff - bd_diff_ni;
            }
        }
    }
}

As a quick aside - note that we could simply multiply [a, b, c, d] by a 4x4 matrix holding the terms we derived in the previous section. Possibly useful for a GPU implementation!

With this we pick up another ~10% speedup, and actually beat the reference implementation’s average-case error!

A few remarks:

• Obviously, this can only be used on inputs with a size that’s a power of 4. We could get fancy and do some split-radix stuff (some stages are radix a, others radix b, etc.) but I’m not planning to.
• This is faster because we save on complex multiplies and adds.

Opt. 5: Special-case first stage

The twiddles we look up in our inner loop are just $W_N^{k s}$. For the first iteration, big_n $= 4$, so $k$ is always 0 – therefore, we use $W_N^0$, which is just 1. We can special-case this and save a few more complex multiplies:

fn fft_butterfly_radix_4<T: Float + FloatConst>(
    input: &mut [Complex<T>],
    output: &mut [Complex<T>],
    stride: usize,
    big_n: usize,
    twiddles: &[Complex<T>],
) {
    for start_idx in 0..stride {
        for k in 0..big_n / 4 {
            // Collect inputs.
            let i0 = input[start_idx + 4 * k * stride];
            let i1 = input[start_idx + (4 * k + 1) * stride];
            let i2 = input[start_idx + (4 * k + 2) * stride];
            let i3 = input[start_idx + (4 * k + 3) * stride];
            // Collect relevant twiddles.
            let ot1 = twiddles[1 * k * stride];
            let ot2 = twiddles[2 * k * stride];
            let ot3 = twiddles[3 * k * stride];

            let a = i0;
            let b = ot1 * i1;
            let c = ot2 * i2;
            let d = ot3 * i3;

            // To derive this, write the output assignments in terms of
            // a/b/c/d, then factor out!
            let ac_sum = a + c;
            let ac_diff = a - c;
            let bd_sum = b + d;
            let bd_diff_ni = mul_ni(b - d);

            output[start_idx + k * stride] = ac_sum + bd_sum;
            output[start_idx + (k + big_n / 4) * stride] = ac_diff + bd_diff_ni;
            output[start_idx + (k + big_n / 2) * stride] = ac_sum - bd_sum;
            output[start_idx + (k + 3 * big_n / 4) * stride] = ac_diff - bd_diff_ni;
        }
    }
}

fn fft_butterfly_radix_4_s0<T: Float + FloatConst>(
    input: &mut [Complex<T>],
    output: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    let stride = input.len() / 4;
    let big_n = 4;

    for start_idx in 0..stride {
        for k in 0..big_n / 4 {
            // Collect inputs.
            let i0 = input[start_idx + 4 * k * stride];
            let i1 = input[start_idx + (4 * k + 1) * stride];
            let i2 = input[start_idx + (4 * k + 2) * stride];
            let i3 = input[start_idx + (4 * k + 3) * stride];

            let a = i0;
            let b = i1;
            let c = i2;
            let d = i3;

            // To derive this, write the output assignments in terms of
            // a/b/c/d, then factor out!
            let ac_sum = a + c;
            let ac_diff = a - c;
            let bd_sum = b + d;
            let bd_diff_ni = mul_ni(b - d);

            output[start_idx + k * stride] = ac_sum + bd_sum;
            output[start_idx + (k + big_n / 4) * stride] = ac_diff + bd_diff_ni;
            output[start_idx + (k + big_n / 2) * stride] = ac_sum - bd_sum;
            output[start_idx + (k + 3 * big_n / 4) * stride] = ac_diff - bd_diff_ni;
        }
    }
}

pub fn fft_v5_s0_opt<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    assert!(is_power_of_k(src.len(), 4));
    let n_iter = log_k_of::<4>(src.len());

    dst.copy_from_slice(src);

    let (mut input, mut output) = if n_iter % 2 == 0 {
        (dst, src)
    } else {
        (src, dst)
    };
    let big_n = input.len();
    let mut stride = big_n;
    let mut big_n = 1;
    for stage in 0..n_iter {
        stride /= 4;
        big_n *= 4;
        std::mem::swap(&mut input, &mut output);

        if stage == 0 {
            fft_butterfly_radix_4_s0(input, output, twiddles);
        } else {
            fft_butterfly_radix_4(input, output, stride, big_n, twiddles);
        }
    }
}

Here I refactored the inner loop of our FFT - called a butterfly in FFT research parlance - and made a variant which avoids those complex multiplies in stage 1. We get a few more microseconds out of this, with no change to our accuracy:

Within 12% of our speed-of-light! No, we’re not done yet :)

Opt. 6: Unsafe

Our butterfly does 8 array lookups, each of which Rust will bounds-check for us. However, we know by inspection that they will never go out of bounds. So we can tell Rust this with the unsafe keyword, and enable more compiler optimizations.

fn fft_butterfly_radix_4_unsafe<T: Float + FloatConst>(
    input: &mut [Complex<T>],
    output: &mut [Complex<T>],
    stride: usize,
    big_n: usize,
    twiddles: &[Complex<T>],
) {
    let input_ptr = input.as_ptr();
    let output_ptr = output.as_mut_ptr();
    for start_idx in 0..stride {
        for k in 0..big_n / 4 {
            unsafe {
                // Collect inputs.
                let i0 = *input_ptr.add(start_idx + 4 * k * stride);
                let i1 = *input_ptr.add(start_idx + (4 * k + 1) * stride);
                let i2 = *input_ptr.add(start_idx + (4 * k + 2) * stride);
                let i3 = *input_ptr.add(start_idx + (4 * k + 3) * stride);
                // Collect relevant twiddles.
                let ot1 = twiddles.get_unchecked(1 * k * stride);
                let ot2 = twiddles.get_unchecked(2 * k * stride);
                let ot3 = twiddles.get_unchecked(3 * k * stride);

                let a = i0;
                let b = ot1 * i1;
                let c = ot2 * i2;
                let d = ot3 * i3;

                // To derive this, write the output assignments in terms of
                // a/b/c/d, then factor out!
                let ac_sum = a + c;
                let ac_diff = a - c;
                let bd_sum = b + d;
                let bd_diff_ni = mul_ni(b - d);

                *output_ptr.add(start_idx + k * stride) = ac_sum + bd_sum;
                *output_ptr.add(start_idx + (k + big_n / 4) * stride) = ac_diff + bd_diff_ni;
                *output_ptr.add(start_idx + (k + big_n / 2) * stride) = ac_sum - bd_sum;
                *output_ptr.add(start_idx + (k + 3 * big_n / 4) * stride) = ac_diff - bd_diff_ni;
            }
        }
    }
}

fn fft_butterfly_radix_4_s0_unsafe<T: Float + FloatConst>(
    input: &mut [Complex<T>],
    output: &mut [Complex<T>],
) {
    let stride = input.len() / 4;
    let big_n = 4;
    let input_ptr = input.as_ptr();
    let output_ptr = output.as_mut_ptr();
    for start_idx in 0..stride {
        for k in 0..big_n / 4 {
            unsafe {
                // Collect inputs.
                let i0 = input[start_idx + 4 * k * stride];
                let i1 = input[start_idx + (4 * k + 1) * stride];
                let i2 = input[start_idx + (4 * k + 2) * stride];
                let i3 = input[start_idx + (4 * k + 3) * stride];

                let a = i0;
                let b = i1;
                let c = i2;
                let d = i3;

                // To derive this, write the output assignments in terms of
                // a/b/c/d, then factor out!
                let ac_sum = a + c;
                let ac_diff = a - c;
                let bd_sum = b + d;
                let bd_diff_ni = mul_ni(b - d);

                *output_ptr.add(start_idx + k * stride) = ac_sum + bd_sum;
                *output_ptr.add(start_idx + (k + big_n / 4) * stride) = ac_diff + bd_diff_ni;
                *output_ptr.add(start_idx + (k + big_n / 2) * stride) = ac_sum - bd_sum;
                *output_ptr.add(start_idx + (k + 3 * big_n / 4) * stride) = ac_diff - bd_diff_ni;
            }
        }
    }
}

pub fn fft_v6_unsafe<T: Float + FloatConst>(
    src: &mut [Complex<T>],
    dst: &mut [Complex<T>],
    twiddles: &[Complex<T>],
) {
    assert!(is_power_of_k(src.len(), 4));
    assert_eq!(src.len(), dst.len());
    assert_eq!(twiddles.len(), src.len());
    let n_iter = log_k_of::<4>(src.len());

    dst.copy_from_slice(src);

    let (mut input, mut output) = if n_iter % 2 == 0 {
        (dst, src)
    } else {
        (src, dst)
    };
    let big_n = input.len();
    let mut stride = big_n;
    let mut big_n = 1;
    for stage in 0..n_iter {
        stride /= 4;
        big_n *= 4;
        std::mem::swap(&mut input, &mut output);

        if stage == 0 {
            fft_butterfly_radix_4_s0_unsafe(input, output);
        } else {
            fft_butterfly_radix_4_unsafe(input, output, stride, big_n, twiddles);
        }
    }
}

Note that “add” just means “add a value to this pointer.” Seems to be the canonical way to do pointer arithmetic in Rust. With this, we have nearly reached the speed of light!

No, we’re not done.

Opt. 7: Radix-8

Why stop at radix-4? If we extend to radix-8, we still get the desirable analytic property of our factors not requiring complex multiplies, as they’re just 45-degree rotations, but we also reduce the number of stages.

For a length-4096 input, aka $2^12$, radix-4 requires 6 stages, where radix-8 requires only 4. If each stage does 3 and 7 complex multiplies respectively, we wind up with 18$s$ vs. $14$ total complex multiplies.

We also reduce the number of times that we need to read the full data buffer from 6 to 4.

I can derive the radix-8 twiddles by inspection - it’s left as an exercise to the reader if needed. (Tip: visualize the rotations through the complex plane.)

Let $p = \frac{1}{\sqrt{2}}$. Then: