JoshJers' Ramblings

This is part 3 of a series of posts about FFTs and the specific algorithms that FFT-20XX uses.

1. The Discrete Fourier Transform
2. The Basics of the FFT
3. Optimizations to the FFT (← you are here)
4. The Inverse FFT
5. The FFT-20XX Complex-Input Algorithm
6. The FFT-20XX Real-Input Algorithm(s)

In the last post we got to a basic implementation of the radix 2 decimation-in-time Cooley-Tukey algorithm. But there are still other additions and optimizations we can make before getting into the way FFT-20XX works. We’ll start by diving into the “radix 2” part of the algorithm: we can make improvements there.

Butterflies

Before we go into the specifics, I’m finally going to bring in a classic FFT diagram: the butterfly, an FFT data flow diagram. The most basic one is for a FFT of length 2, which has only a single step over two inputs. Here’s a diagram for the radix-2 decimation-in-time 2-length FFT:

In this diagram, the inputs are on the left, each computation (just one, in this example) is represented as a small circle, and the outputs are on the right. The computation (for a decimation-in-time FFT) rotates the bottom input by the given rotation value; in this diagram it’s just a $0$, but in general the value there will be some $t$ where the actual rotation angle is $t/N$ clockwise turns ($-2\pi t/N$ radians).

Note

The formatting of these butterfly diagrams is a little different than a standard butterfly diagram, which tend to look more like the following:

These diagrams are more comprehensive (the $w^k_N$ value makes the whole turn fraction of $k/N$ clockwise turns clear, where in my diagrams the “$/ N$” is implied), it has a $-1$ multiplier on the lower input going into the sum circle (to signify that the bottom “sum” is a difference), and the sums are each on their own line as opposed to a central circle that represents the add/subtract pair, but in practice – especially on lower-resolution displays – my variant is, imo, less visually noisy.

Plus, once you get into the higher radixes, even the standard diagrams start to be drawn more simply.

Now here’s one for the decimation-in-time 4-length FFT (where the bit reverse is performed on the input):

There are two “phases” in this radix-4 FFT (two runs of the outer loop of the code), which are separated via the vertical dotted lines. Really, though, this is four of the previous (length 2 FFT) butterflies bolted together as two separate 2-element computations (the first pass) followed by a single 4-element computation loop (the second pass, itself just two 2-element computations stacked on top of each other wearing a trenchcoat).

Finally, here’s a butterfly diagram for a full 16-element FFT:

There are four phases in this one, and note that you can still see the structure discussed in the last post of 8 2-element computations in the first phase, then four 4-element computations in the next, then two 8-element computations, then a single 16-element computation. Also note the rotation angles in each phase: the last phase has rotations 0-7 (corresponding to $0/16$ to $7/16$ turns), the previous one has just the even numbers, the one before that just 0 and 4 (every other even value), down to no rotations in the first phase.

Note

The first radix-2 phase of any decimation-in-time FFT will always have no rotations, and the second phase will only have 0- or 90-degree rotations (a $4$ in the above diagram, as it’s $\frac{4}{16}$ turns, or $\frac{1}{4}$). When it comes time to start optimizing, these first two phases can be done with no multiplications at all (obvious when there’s no rotation, but a 90 degree (clockwise) rotation is just taking $a + ib$ and changing it to $b - ia$), and can typically be special- cased for a significant speed bost.

Doubling the Radix

We went with a radix 2 breakdown in the last post because it’s the most straightforward, but there are others. The one we’re going to focus on here is radix 4, which, in contrast to the above 16-element butterfly diagram, would instead look like the following:

Don’t worry about the specifics of the rotation values or what this four-in/four-out computation node is doing yet – we’ll get to that shortly. For now, the important thing is that there are now half as many passes over the data, and each computation handles four values at a time instead of two.

This will end up being more efficient in a couple different ways. The main efficiency gain is that having half as many passes means it halves the number of reads from and writes to memory, which as the FFT lengths get larger can make a big difference in how often there are cache misses. Second, it is possible to reduce the number of rotations needed by $\frac{1}{4}$, which can also make a difference on some CPUs (but not all, which we’ll get to).

To derive a radix 4 solution, you could go back to the original equation, split the inputs into four interleaved sections (where $k\mod 4$ is $0$, $1$, $2$, and $3$ respectively), then split the outputs into quarters, and redo all of the rotation shenanigans that we did when deriving radix-2.

Nah. I’m not doing all that.

Instead, we’re going to start with the radix-2 diagram and build up to what that four-in/four-out computation looks like. Every pair of radix-2 passes can be turned into a single radix-4 pass, so let’s start by highlighting a couple such groupings in our radix-2 diagram:

Each of the groupings of “two pairs radix-2 passes that all use the same 4 values” can be represented generically as:

The input/output indices are relative to the position in the FFT, but note that the middle two indices are in the opposite order as the output indices: this will always be true when we’ve done the bit reverse at the outset, and remains true no matter which intermediate step we’re on (i.e. it just automatically happens that way).

Also note that the angles themselves are relative, too: there will be some base angle $t$ (in the first highlighed group at input indices 0, 8, 4, and 12, $t$ is $0$. In the second at output indices “2, 6, 10, and 14”, $t$ is $2$) that the other angles are relative to, and we’ll set it as the first angle in the second radix-2 pass. Both angles in the first radix-2 pass are double this $t$ value (both labeled $2t$), and the second angle in the last pass is always a quarter turn more than $t$ (labeled $t + \frac{N}{4}$).

To get from here to a full radix-4, we can rearrange the angles a bit. If we separate the $t$ from the $\frac{N}{4}$ in the last angle (taking advantage of the fact that, since we’re using complex numbers, Rotate(a + b) is the same as Rotate(a) * Rotate(b)), there is a common $t$ factor to both of the second radix-2 angles, and we can distribute that back to the latter two inputs to the first radix-2 step (the first of which has no rotation initially so it just becomes a rotation of $t$, and the latter of which becomes effectively Rotate(2t) * Rotate(t) which, as just explained, becomes Rotate(2t + t) or Rotate(3t):

As mentioned in an above note, a quarter turn (90-degree clockwise rotation) can be done without any sine/cosine multiplications at all (it transforms $a + ib$ into $b - ia$), so this is mathematically more efficient: it does 3 rotations on 3 inputs at the start of the process instead of 4 rotations (split between the first and second inner radix-2 phases). The number of adds/subtractions remains the same, but that’s okay: there’s no way to reduce it any further.

This, then, is what becomes the computation node of the radix-4 butterfly diagram we made above! This, then, is the equivalent radix 4 diagram to the previous diagram (which hides the internal complexities, including the effectively-free quarter turn):

Note

An implementation detail worth mentioning here: on some systems using the first 4-element breakdown from above (the one with the $t + \frac{N}{4}$ term) can be more efficient! There are systems that have FMA instructions (the fused multiply-add) which have a cost close to doing either a multiply or an add by itself (x64 CPUs with the FMA extension – like most with AVX2 instructions – have this property). Doing it that way (which is sometimes referred to as radix-$2^2$), even with what seems like an extra complex multiply, can end up being cheaper due to less overall operations and one less required angle (the $3t$ angle isn’t required, using $t$ and $2t$ is sufficient).

There’s a breakdown of this trick on The ryg blog: Notes on FFTs: for implementers.

Here’s the 16-element radix-4 diagram we used earlier, with our same two radix-4 sections highlighted:

Note

This formulation of radix-4 is slightly different than the canonical one, which has a different “scrambled” input or output order from the bit reverse that we’re using. In that case, the reordered indices are reversed pairs of bits: for instance, an 8-bit binary value with bits 01234567 would becomes 67452301 instead of 76543210.

Instead, we’ve kept the radix-2 style bit reverse, the consequence of which is the reversal of the middle 2 inputs in the above diagrams (with reversed bit pairs, those would be in proper order, and the angles would instead be $t$, $2t$, and $3t$, respectively).

Now, hopefully, it’s clear why the angle values are the way they are in this diagram. The first phase has no rotations (except for the free 90-degree rotation), and the second phase has four computations with increasing rotation values of the form $2t$, $t$, $3t$.

The only issue with using the radix-4 algorithm is that it only fully works with FFTs with power-of-four lengths (because it splits in fours each time, so it works for length 4 but not 8, 16 but not 32, etc). The good news is that the way we’ve formulated it (by effectively smooshing together two radix-2 passes into a single radix-4 pass), you can just put a single radix-2 pass in there to make up the difference (and it doesn’t even matter where!)

Here is an 8-length FFT where the first pass is still radix 2 but the second is radix 4:

And here’s the same 8-length FFT the other way around: radix-4 then radix-2:

For longer FFTs, the radix-2 pass can be anywhere in the middle, as well - whatever is most efficient for your implementation! I found it easiest to place it right after the very first radix-4 pass (i.e. like that last diagram), mostly for simplicity.

Other Radix Flavors

So if radix 4 is more efficient than radix 2, that raises the question: are there other, even more efficient radixes? radix-8? radix-16? How far is too far?

Radix-4 is the largest power-of-two radix that has “free” rotations in the middle (that quarter turn). Radix-8 has some 45-degree-multiple rotations in the middle, which still require some multiplication. However, because both cosine and sine of 45 degrees (and related friends) are both $\frac{1}{\sqrt{2}}$, the multiply can be done as a common factor so it still is somewhat more efficient.

The problems with radix-8 tend to have to do with CPU register limits: you need to load 8 complex values (so 16 scalars), then (similar to how radix-2 needs 1 and radix-4 needs 3) it tends to need 7 complex twiddles plus the $\frac{1}{\sqrt{2}}$ internal value (so another 15 scalars) and you tend to also need some scratch registers as well to do computations. AVX and AVX2 only have 16 SIMD registers so even before you get to the angles you’ve already run out of register space. AARCH64 has 32, so you could maybe squeeze in there, but it would be tight. Plus you’re going to be up against CPU cache associativity for systems with eight-way associative caches, because you’re loading from 8 locations but also need to be loading the angles from somewhere.

Past radix-8, you’re out of register space, hitting CPU cache associativity limits, and also have more internal rotations, so it’s going to perform way worse under most circumstances.

Note

One place where it’s really tempting to use radix-8, however, is on the very first decimation-in-time step, which has no twiddles (all seven of the twiddle angles are 0), and thus only has those nice-to-compute internal 45-degree rotations. However, the register limits are still an issue and, depending on how you have to load your values (i.e. for SIMD purposes you may want to do a deinterleave to separate the real and imaginary components), you may still not have enough space. But it’s worth considering, if you can.

There is another radix option worth mentioning: the split-radix FFT algorithm, which I’m not going to detail here. In terms of required math operations, the split-radix algorithm gets you to some of the lowest-known possible operation counts. Unfortunately, in practice, I found that it breaks up the nice regularity of using just a standard radix-4 algorithm and ultimately – at least in my attempts at implementing it – had worse performance overall.

Okay, we’re through all of the radix shenanigans!

Lookup Tables

The last thing to touch on in this post is doing lookup tables for the angles. These tend to be exactly what they sound like: you precompute a table for your required FFT length that contains all of the sin/cos values needed for that length.

The simplest lookup table is one where you have $\frac{N}{2}$ cos/sin values ($0$ turns through $\frac{N}{2} - 1$ turns), in order. However, there are some potential ways to make the table a little more efficient (with size, cache-friendliness, or both):

• You don’t need both $\frac{N}{2}$ $\cos$ and $\frac{N}{2}$ $\sin$ values (unless you’re doing SIMD with interleaved real and imaginary components): because $\cos(t) = \sin(\frac{\pi}{2} + t)$, you can store $\frac{3N}{4}$ scalar values and start the cosine lookup partway into the table, and use $\frac{1}{4}$ less memory.
• if you’re doing radix-4, you can avoid storing the $3t$ angles and instead compute them using the $t$ and $2t$ rotations, which saves a bunch of table space. This may sound like you’re adding a complex multiply back in after eliminating it, but you can move that computation into an outer loop so it isn’t an issue.
• It can also be convenient – for SIMD loads and CPU cache-friendliness – to store not only the angles in order for size $N$, but to store all of the table sizes from some minimum length (like 4 or 8) up to your max length.
  • Doing this is nice because it means you can build a single table for your maximum required FFT length and it will perform just as well for any smaller FFT lengths as well).
  • Also, the earlier passes of the algorithm can take advantage of these smaller tables (the second radix-2 pass only needs the angles for a length-4 FFT, and the next only needs the angles for a length-8 FFT, etc). Loading from them can help with CPU cache coherence (as all of the required angles at that level are contiguous), plus SIMD loading if you need multiple angles within a single pass.

Note

It’s also worth noting, for some algorithmic variants (i.e. if you do the bit reverse at the end, or if you are doing a decimation-in-frequency implementation), it can be better to store the table in bit-reversed order (for example, for a 16-length FFT you’d store the rotations in order 0 4 2 6 1 5 3 7)! This has a neat property: if you make a table for a length-$N$ FFT, the first $N/2$ entries are also exactly what would be in a table for a length-$N/2$ FFT! This means the table can be built for the maximum required FFT length, but unlike the last bullet point above, with no extra memory footprint.

For FFT-20XX I found all of the above bullet points to be useful.

What Now?

No code this time, but using the above tricks (namely: using radix-4 instead of radix-2 where possible and switching to lookup tables for the rotations) gets you to a somewhat decent baseline FFT implementation! There are many more places to optimize: figuring out SIMD support and optimizing the bit reversal pass (or – spoiler alert – eliminating it altogether). We’ll get to those once we finally get into FFT-20XX’s algorithms.

Before that, next time, we’ll touch on the decimation-in-frequency version of the algorithm on our way to finally talking about the inverse FFT and how to compute it!