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!

This is part 2 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 (← you are here)
3. Optimizations to the FFT
4. The Inverse FFT
5. The FFT-20XX Complex-Input Algorithm
6. The FFT-20XX Real-Input Algorithm(s)

Part 1 went through the core of how the discrete Fourier transform (DFT) works to turn a series of samples in time into a set of samples in frequency, and ended up with the following algoritm:

void CalculateDFT(
  complex[] inputs,
  complex[] outputs,
  int length)
{
  for (outIdx = 0; outIdx < length; outIdx++)
  {
    // Start the sum at 0.
    float sum = 0

    // This is the angular frequency (in radians/sample) for
    // the given output sample.
    float angularFreq  = -2 * pi * outIdx / length

    for (inIdx = 0; inIdx < length; inIdx++)
    {
      // Get the rotation value for the current input index
      // then use that to calculate the resulting rotation.
      float radians = inIdx * angularFreq
      complex rotation = cos(radians) + i*sin(radians)

      // A complex multiply with a unit complex vector is
      // a 2D rotation.
      sum += inputs[inIdx] * rotation
    }

    outputs[outIdx] = sum
  }
}

Note

A reminder: we’re going to be talking about these values as complex numbers (of the form $a + i b$), but if you prefer, just think of them as a 2D vector instead ($[a, b]$ where $a$ is the x coordinate and $b$ is the y coordinate).

We’ll be doing complex multiplication of these vectors with others, but in all cases a multiply like $P \cdot R$ will be a rotation of point $P$ by the angle that $R$ represents (where $R = \cos(\theta) + i\sin(\theta)$).

This algorithm is $O(N^2)$, but there’s a much more efficient, $O(N\,\log\,N)$ way: the fast Fourier transform (FFT) – and we’re going to derive it.

Note

The algorithm we end up deriving here is going to work specifically for FFT lengths that are powers of two, because it’s the easiest to implement (plus, to be honest, power-of-two lengths were sufficient for FFT-20XX so that’s all I implemented). The following techniques can be adapted for other lengths (i.e. if you need a length that is a mutliple of 3 or 5), plus there are ways to do this for arbitrary lengths (for instance: Bluestein’s algorithm). This series isn’t going to touch on any of those.

Halving the Work

The first step of deriving the FFT algoritm is figuring out where we can halve the amount of work. Since the DFT boils down to “every output takes every input multiplied by a rotation”, we can visualize all of the rotations as a 2D grid, where each column is an input and each row is an output.

Here’s a visualization of the rotations for an FFT of length 8 (adding a dividing line between the top and bottom halves, which should make sense in a moment):

Each rotation (as per the above pseudocode) for a given input index $n$ and output index $k$ is $n k/N$ turns (where $N$ is 8, the FFT length). But rotations are cyclical, and, for instance, turning a full turn (360 degrees) is equivalent to not turning at all.

We can use this cyclicity to our advantage! If you look at the even inputs (columns indexed 0, 2, 4, and 6), you’ll note that the rotations in the bottom half are identical to the top half. For instance, with input 2, outputs 0 and 4 are both unrotated, outputs 1 and 5 both rotate a quarter turn, and so on.

For the odd inputs (columns 1, 3, 5, 7), the angles aren’t the same. However, they’re still related: the bottom half rotations are always 180 degrees from the corresponding top half rotations. As an example, looking at the column for input 1, we see that output 4’s rotation is a half turn, where output 0 is unrotated – they’re facing away from each other. The other odd outputs are similar, with the bottom half rotations always pointing the opposite direction as their top-half counterparts.

In other words, the bottom half rotations are negated versions of the top half (because with 2D vector A, Rotate(A, 180 degrees) is the same as the 2D vector -A).

Let’s rearrange the columns of the diagram to group the even and odd elements together (evens to the left, odds to the right):

We can then rewrite the bottom-half odd inputs as negations of the top-half ones (note the minus signs in place of the usual plus signs):

Now all of the corresponding top and bottom half outputs use the same corresponding rotations per input; the only difference is that in the bottom half we subtract the rotated odd inputs instead of adding them. Our code, then, can switch to summing up all of the even and odd elements separately, then computing two outputs at once by doing evenSum + oddSum for the first (top) half of the outputs and evenSum - oddSum for the second (bottom) half, which means we’re now doing half the work as before!

Note

If you want this (and the next section) done algebraically instead of visually, the wiki page for the Cooley-Tukey algorithm has a fairly clear breakdown of the derivation of the whole thing.

Here’s what the code looks like if we do that:

void CalculateDFT(
  complex[] inputs,
  complex[] outputs,
  int length)
{
  // NOTE: Now looping over half the length
  for (outIdx = 0; outIdx < length / 2; outIdx++)
  {
    // Start the sums at 0.
    float evenSum = 0
    float oddSum = 0

    // angular frequency for the given output
    // (in radians/sample)
    float angularFreq  = -2 * pi * outIdx / length

    // NOTE: now incrementing by 2 to do odds and evens
    // separately.
    for (inIdx = 0; inIdx < length; inIdx += 2)
    {
      // Get the rotation value for the current input index
      // then use that to calculate the resulting rotation.
      float evenRads = inIdx * angularFreq
      complex evenRot = cos(radians) + i*sin(radians)
      evenSum += inputs[inIdx + 0] * evenRot

      // Do the same with the next index (inIdx + 1), for
      // the odd sum.
      float oddRads = (inIdx + 1) * angularFreq
      complex oddRot = cos(radians) + i*sin(radians)
      oddSum += inputs[inIdx + 1] * oddRot
    }

    // Computing 2 outputs in the same loop!
    outputs[outIdx]            = evenSum + oddSum
    outputs[outIdx + length/2] = evenSum - oddSum
  }
}

If you compare this code to what we started with, it does half the rotations and half the adds (counting subtractions) as the original: the inner loop does effectively the same amount of work in both, but the outer loop runs half as much. Unfortunately, it’s not quite to $O(N\,\log\,N)$ yet, so we need a way to effectively do this halving of work recursively.

Halving Work Again (and Again and…)

Let’s look again at just the first half of the outputs from our 8-length example:

If we look at just the even inputs for a moment, the pattern of rotations in that 4×4 block of the diagram ends up being the exact same pattern of rotations as for a standard 4-element FFT (half the size of our 8-element example):

Note

If you care about the algebra: this is because for each even rotation ($\frac{n k}{N}$), $k$ is a multiple of 2. So if we say $2m = k$, then the rotation becomes $\frac{2 n m}{N}$ which is equivalent to $\frac{n m}{N/2}$, which is exactly how the angles for an FFT of length $N/2$ would be specified.

This means that we can calculate the even sums recursively, computing the $N/2$ even sums as if they were an $N/2$-length DFT, which means we can do the same even/odd dance to halve its work as well.

That’s great for the even inputs, but what about the odd inputs? They have a similar pattern, but it’s not quite the same:

The first row is the same as in the 4-element DFT: no rotation. The second row, however, is different: each angle in the row is 1/8th of a turn rotated relative to the corresponding row of the 4-element DFT. The third and fourth rows are also different, with an extra 2/8th and 3/8th turn added to the angles in each row, respectively.

In other words, any given output $k$’s odd inputs will have their rotations be $k/N$ turns more than doing a standard $N/2$-length FFT with those inputs.

Since “rotating every vector by t then summing” is the same as “summing then rotating by t” (that is, Rotate(A, t) + Rotate(B, t) equals Rotate(A + B, t)), we can calculate each odd sum by doing the work as a half-length DFT (like we can with the evens), then adjust each sum by rotating by its corresponding $k/N$-turn angle (called a twiddle factor, which is a fun name to say), like so:

Using that, we can now recursively calculate the even and odd sums using half-length FFTs, apply the twiddles, then do the add and subtract to get the final outputs.

Let’s write this as code again, this time recursively! To do this, for now we’ll use some temporary storage for simplicity:

// Helper function to build a complex
// rotation value
complex Twiddle(int k, int len)
{
  float angularFreq = -2 * pi * k / len
  return cos(angularFreq) + i*sin(angularFreq)
}

void CalculateDFT(
  complex[] inputs,
  complex[] outputs,
  int len)
{
  if (len == 1)
  {
    // To end the recursion, special-case
    // the 1-length FFT, which sums one input,
    // unrotated (i.e. it does nothing).
    outputs[0] = inputs[0]
    return
  }

  // Temp storage for the recursion
  complex[length / 2] evens
  complex[length / 2] odds

  // Separate the even and odd inputs
  for (k = 0; k < len / 2; k++)
  {
    evens[k] = inputs[2 * k]
    odds[k]  = inputs[2 * k + 1]
  }

  // Recurse into the even and odd FFTs
  // (each using the same array as input and
  // output, which is fine here)
  CalculateDFT(evens, evens, len/2)
  CalculateDFT(odds,  odds,  len/2)

  for (k = 0; k < len / 2; k++)
  {
    // Get the corresponding recursively-calcualted
    // sum from the relevant arrays.
    complex evenSum = evens[k]
    complex oddSum = odds[k]

    // The odd sum needs to be rotated by this
    // output's twiddle factor.
    oddSum *= Twiddle(k, len)

    // Now, as before, calculate two outputs from
    // these sums.
    outputs[k]         = evenSum + oddSum
    outputs[k + len/2] = evenSum - oddSum
  }
}

This is now truly, algorithmically $O(N\,\log\,N)$! We halve the work at each level of the FFT, which does way fewer rotations and adds than the initial DFT algorithm (for instance, the original algorithm for a 1024-length DFT would do around a million rotations, but this one does around five thousand).

But it’s not ideal: it’s recursive, and uses temporary storage per recursion. We can do better! But now we get into the fun that is…

Deinterleaving and Bit Reversal

To do all of this without recursion or scratch space, we need a way to do each level of FFT in-place. So let’s look at how data is flowing through the system.

Here’s a diagram of the routine (at a single recursion level) where, from left to right:

• The inputs are split into evens and odds.
• Each of those two groups has a half-length sub-DFT performed on it.
• Those results are then used in pairs by each step of the computation loop to write corresponding pairs of outputs.

If you look at the flow through this diagram, what’s interesting is that in the box on the right (representing the loop) each computation pulls in from the same row that it writes out to (for instance, DFT result indices ‘0’ and ‘1’ are lined up perfectly with final outputs ‘0’ and ‘8’).

This means that, theoretically, we could compute the entire FFT in-place (where the inputs and outputs are the same array), by doing the following:

1. Deinterleave the data in-place (separating the evens and odds).
2. Calculate the sub-DFTs recursively (again, in-place), using the first and last halves of the data as the even and odd sub-DFTs, respectively.
3. Then, run the computation loop over the results, each iteration of which reads from and writes to the same elements, so no rearranging needs to take place.

There’s one problem with this: it turns out, efficiently deinterleaving an array in-place is one of those problems that sounds like it would be super easy, but is instead kind of a nightmare!

To try to avoid that, let’s keep investigating the pattern of data moves. We know that, in terms of data flow, the data into the last step comes in from two half-length sub-FFTs:

What if we track this pattern of deinterleaves all the way back to the lowest level (the two-element FFTs), where the indices from right to left (from output to input) get deinterleaved at each step? Here’s the ordering we end up with:

At a glance, the pattern of inputs on the left isn’t necessarily obvious, but it turns out it’s straightforward: the values are the original indices with their bits reversed!

For a 16-element FFT (like the diagrams above), we have 4 bits worth of index (since $\log_2 16 \implies 4$), so here’s the table of reversals:

As you can see, the ordering on the right of this table (0, 8, 4, 12, etc) is the same as the input ordering from the previous diagram and, in their bit representations, you can see that each index is just the 4 bits of the other reversed.

Here is a simple routine that does a bit-reversed copy from one array to another:

void BitReverseCopy(
  complex[] inputs,
  complex[] outputs,
  int len)
{
  for (int k = 0; k < len; k++)
    outputs[k] = inputs[BitRev(k, len)]
}

Note

The above code doesn’t have an implementation of BitRev, because an efficient implementation on many systems is sorta unreadable; it turns into a sequence of swapping increasingly-large groups of bits until you reach the size of the int. For a C++20 implementaton, here’s a compiler explorer link.

Additionally, the above routine doesn’t work in-place. It is totally possible to do this bit reversal in place: it can be written efficiently as a square matrix transpose (or two) if you get creative with the indexing of the “matrix” rows. However, the final FFT-20XX algorithms do not have an explicit bit reversal step, so I’m not going to go into those details. Sorry!

Going Iterative

We now have all of the pieces we need to assemble the final algorithm. Let’s go back to this earlier diagram of the data flow for a 16-element FFT:

This diagram takes the inputs in what we now know to be bit-reversed order, and does four ($\log_2 16$) passes on the inputs:

• The first pass does eight 2-element computation loops
• The next pass does four 4-element loops
• Next is two 8-element loops
• Finally, a single 16-element loop to get the final outputs

Note

In each of the above passes, the number of computation loops multiplied by the number of elements in the loop is 16 ($8\times2$, $4\times4$, $2\times8$, and $1\times16$), which means each pass processes $N$ elements. Since there are $\log_2 N$ passes of $N$ elements each, this algorithm – like the recursive one – is our target $O(N\,\log\,N)$.

This process can (finally!) be done iteratively instead of recursively, by getting the elements into the expected bit-reversed order, then modifying each element in place, in passes:

void CalculateDFT(
  complex[] inputs,
  complex[] outputs,
  int len)
{
  if (len == 1)
  {
    // Special case length 1.
    outputs[0] = inputs[0]
    return
  }

  // Do the bit reversal so we can do the rest
  // iteratively in the outputs, in-place.
  BitReverseCopy(inputs, outputs, len)

  // Outer loop iterates all of the subFFT lengths:
  // starting at 2, then doubling until we reach the
  // full length ... this is log2(len) iterations.
  for (int subLen = 2; subLen <= len; subLen *= 2)
  {
    // This loop iterates over each subsection, as the starting
    // index into the data.
    for (int subStart = 0; subStart < len; subStart += subLen)
    {
      // The data for this subsection starts at subStart.
      var *data = &outputs[subStart]

      // Finally, do the computation loop for this subsection.
      // The logic here is unchanged from the previous code
      // examples, but the reads and writes are to the same
      // pairs of elements.
      for (int k = 0; k < subLen / 2; k++)
      {
        int kOdd = k + subLen / 2
        complex evenSum = data[k]
        complex oddSum  = data[kOdd]
        oddSum *= Twiddle(k, subLen)
        data[k]    = evenSum + oddSum
        data[kOdd] = evenSum - oddSum
      }
    }
  }
}

Note

The bit reversal does not necessarily have to come at the beginning: the algorithm could instead be written to take the inputs in natural order and end up with the outputs in bit-reversed order, and do the bit reversal in-place as the last step. The indexing within the passes would be different, but the algorithm would otherwise work the same.

In fact, some applications can just deal with the frequency-space outputs being out of order, and so can skip the bit-reversal step entirely! This is really nice when it works out.

What we’ve ended up with is the radix-2 decimation in time Cooley-Tukey algorithm. It’s radix-2 because each “recursion” splits into two sub-FFTs, and it’s decimation in time because each pass deinterleaves the input, time domain samples.

Note

If you go back up to the rotation diagrams again, you might note that the rotations are diagonally symmetrical, which means you could instead split on even/odd output for the first and last half of the inputs (instead of the other way around, like we did). This will end up as a similar but different algorithm, and is the decimation in frequency version of an FFT.

Room for Improvement

While we’ve now gotten to the core FFT algorithm, this still isn’t anywhere near an optimal implementation, and there are a number of things we’ll want to do to get this even faster. For instance:

• As mentioned above, this is a radix 2 algorithm, but there are other radixes (such as radix-4, where each FFT splits into 4 sub-FFTs) that are more efficient and worth investigating.
• The sines and cosines could be done via a lookup table instead of calculating them at runtime, since for a given FFT length there are only so many unique angles that are required.
• Doing the bit reverse pass to get things into the correct order (either at the start or end) can be a sizeable chunk of time, and depending on implementation can wreak havoc on your CPU’s caches.
• And then there’s doing SIMD optimizations and otherwise trying to make the memory accesses CPU-cache-friendly.

We’ll touch on the first couple of those next time!

Because the world clearly doesn’t already have enough of them, I recently released a new FFT library: FFT-20XX. To the best of my ability to tell, it uses an algorithm that is at least slightly novel, and I thought it might be interesting to write about it.

This is the first of a series of posts about FFTs and the specific algorithms that FFT-20XX uses:

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

As a first step I thought it might be fun to start with “what even is the discrete Fourier transform?”

Fourier Transforms

Any explanation of a Fourier transform would be remiss to not include a link to this fabulous 3Blue1Brown video that nicely illustrates how the Fourier transform works. I highly recommend watching it!

This post, however, is specifically about the discrete Fourier transform (referred to as the DFT), which is a Fourier transform that operates on sampled data (like, say, digital audio or images) instead of continuous data. The FFT (fast Fourier transform) is the standard algorithm used to compute the result of a DFT – we’ll get to it in the next post.

What Does the DFT Even Do?

Any periodic signal (i.e. a signal that repeats, such as looped audio) can be represented as the sum of a set of sine waves at different frequencies. Here is an example of a simple signal that can be broken down into three waves:

This is true of any periodic signal whether it is analog or sampled – no matter how complex it may seem.

Note

The “periodic signal” bit is important, but there are many applications that use the DFT on non-periodic signals (like small blocks of a longer piece of audio) and there are a few tricks that can be done (like windowing the signal) to effectively “pretend” that it is periodic when it’s not.

The DFT answers a simple question: given a set of samples that represents a periodic signal, which frequencies are represented in that signal?

To put it another way: the DFT converts a signal that is made up of samples in time (where each sample represents a point on the signal at a given time) into the same signal represented in frequency (where each sample represents the amplitude (how tall the wave is) and phase (the side-to-side position of the wave) of its corresponding frequency). These two representations are frequently referred to as time domain and frequency domain, and both encode the same signal in different ways.

Note

There is also an inverse DFT (the IDFT) which converts the other way – from the frequency domain back to the time domain. Many signal processing applications convert from time to frequency space to do some operation, then use the inverse DFT to convert back.

For a given repeating signal that is N samples long, the DFT returns contribution values for N frequencies, each with one more cycle than the previous. To put it a different way: for each frequency with index i (starting at 0 and ending at N - 1) there are i cycles of a wave (so index 0 has 0 cycles, index 1 has 1 cycle, etc).

Here are the 8 waves for a DFT with 8 samples: (the red wave with circles is the cosine and the green wave with squares is the sine)

Ultimately, a DFT will give you a set of values that tell you how much each of these frequencies contributes to the signal. Specifically, each frequency will have a value that represents both its amplitude and its phase.

Note

The frequency at index 0 always has cosine values of 1 and sine values of 0, and this ends up measuring how off-center the signal is. This is known as the DC bias or DC offset.

Let’s figure out how to calculate the DFT for a signal!

Frequency Testing

The first step of computing the DFT is to figure out how to get a single frequency out of a signal. Ignoring phase for a moment, this can be done in a surprisingly straightforward (to me, at least) manner: if we multiply each sample in the signal by the cosine value of the test frequency and sum the results, you end up with a value that represents the contribution of that frequency to the signal (which is 0 if the frequency is not present in the signal).

Here is an interactive example where you can set the signal and test frequencies and see the result of this process. For purposes of illustration, the final sum in the diagram is instead an average (otherwise it would be too large to fit in the diagram):

Note that when the frequencies match, the resulting average scales with the amplitude of the signal, otherwise (with one exception we’ll touch on in a moment) it is 0. This makes the average effectively a frequency detector (again, ignoring phase for now)! This is still true even when there are multiple frequencies represented within the signal: this will result in a zero value for frequenices that are not in the signal and non-zero values for those that are.

Note

The mechanism by which this works is that the product of two waves with mismatched frequencies (over a span where they have no fractional cycles) will have the same amount of area above and beneath the curve, and thus the average area is zero.

This is hard to illustrate now, but it will be clearer in a later diagram.

As noted, there’s an exception here, where the multiply-and-sum seems to register a match when it shouldn’t. For example: if you set the Signal Cycles to 1 and the Test Cycles to 15 (or vice versa), it will have a non-zero average even though the frequencies don’t match. In fact, with the exception of frequencies 0 and 8, there are two matches for every other frequency. Why is that?

Well, now it’s time to talk about…

Aliasing and Negative Frequencies

Let’s take a closer look at those two frequencies’ cosine values:

Note that while the waves are different frequencies, the cosine samples are the same! This is because of aliasing, which is where a frequency above the maximum representable frequency ends up sampling as if it were a different, lower frequency.

The frequency above which aliasing occurs is called the Nyquist frequency, which is the highest frequency where there can be both a high and low sample for every cycle. In a signal with N samples, that frequency has N / 2 cycles. For example, with 16 samples the Nyquist frequency has 8 cycles (which you can see have a clear up then down pattern):

Above this frequency, the samples end up aliasing to a lower frequency and cannot be differentiated anymore. Here are some more example pairs of aliased cosines from our 16-sample example:

Note that each pair has matching cosine samples, despite the different wave. Specifically: with a repeating sample count of N, any frequency with $N - k$ cycles has the same cosine values as the frequency with $k$ cycles (that is 1 matches 15, 2 matches 14, etc).

That’s neat! But it’s not the whole picture. We can’t just look at the cosine values, we need to also look at what happens to the sine values as well:

While the cosines are the same in these pairs of frequencies, the sines are negated and appear vertically flipped! This means that a given frequency with $N - k$ cycles has negated sine values from the frequency with $k$ cycles.

The relationship here is that the aliased frequencies in each pair are negated versions of the other, since $\cos(-a) = \cos(a)$ and $\sin(-a) = -\sin(a)$. Thus, frequency $N - k$ has the same angle as the hypothetical frequency at $-k$:

Cosines:

Sines:

In practice, when computing the DFT, these frequencies are referred to via their negative indices instead of as the positive aliased ones.

This explains why all of the frequencies (except DC and Nyquist) have two matches: the cosines match at both the positive and negative versions of the angle. This may seem redundant, but it isn’t always – we’ll get into why in a bit.

Note

It’s also worth noting that the multiply-then-average in the example diagram results in half the signal amplitude for each match, and when both are summed, the result represents the full amplitude!

However, DC and Nyquist – because they each only have a single match – have an average that matches the amplitude exactly.

Measuring Phase

So far we have been looking at a signal that contains a cosine, comparing it with a test wave that is also a cosine. But the original description of the DFT mentions both amplitude and phase as being represented, so what happens if we let the phase of the signal shift around relative to the test?

Oh no! Sliding the phase around causes resulting average to oscillate! While we get the expected 0.5 average when the phase is 0, it is -1 at a phase of 0.5 (a half-wavelength offset) and 0 at phases of 0.25 and 0.75 (quarter- and three-quarter-wavelength offsets)! Clearly there is something else we need to consider here…and it turns out it’s the sine of the test frequency.

Instead of just multiplying by the cosine of the test value, we can instead multiply by both cosine and sine to get two values per sample. If we then treat those two resulting values as a 2D coordinate (where cosine corresponds to x and sine to y), we can average them all to get our final result, which will be at the origin for non-matching frequencies and at some non-origin position for matches.

Note

“Multiply by both sine and cosine” is a simplification for purposes of illustration. It works for now but is not the general form, which we’ll get to in a bit.

Some observations from this diagram:

• If the cycle value magnitudes don’t match, the points are all evenly distributed around the center (and thus average to the origin). This makes it much easier to visualize why we end up with no measured amplitude for mismatched waves.

• If the cycle value magnitudes do match (i.e. a signal of 2 cycles with either 2 or -2 test cycles), the points are off-center, and thus their average is also off-center.

• Increasing the phase causes the resulting average to spin clockwise when the test cycle count is positive, but counter-clockwise when the test cycle count is negative. The magnitude of the average remains the same no matter the phase.

• At Nyquist, the resulting points always lie on the x axis, as the test wave’s sine values are all zeros. It turns out: a phase-shifted wave at Nyquist is indistinguishable from a non-shifted wave with scaled amplitude.

This gives us a new way to think about this computation: rather than thinking of the operation as a multiplication, instead we can think of each sample of the signal being rotated by the angle in the test wave at that sample position to get a new 2D position. Those 2D positions are then averaged to get a value that corresponds to the frequency’s contribution.

The computed amplitude (the strength of this frequency’s contribution) is the magnitude of the resulting average, and the phase is the angle of the average (computed, in this case, using atan2).

2D Coordinates and Complex Numbers

Each result of this multiplication with cosine and sine ends up as a 2D coordinate. Canonically, the Fourier transform represents these as a complex number (of the form x + i*y). Why complex numbers? Honestly I think it’s a hack to get a 2D coordinate that is simple to represent with standard math notation (mathemeticians, don’t @ me).

But, to be fair, complex numbers do have a nice property that standard 2D vectors do not: multiplying a complex number by another is equivalent to scaling the first value by the magnitude of the second then rotating it by the angle of the second.

Thus, multiplying a sample from the signal (whether it is complex or real-valued) by cos(t) + i*sin(t) is equivalent to rotating that sample by the angle t (the rotation value has a magnitude of 1 so no scaling occurs in this case).

Note

This also finally gives us an answer to the earlier question of why we need both the positive and negative versions of a given frequency: if the input samples are complex instead of purely real-valued, the contributions to these paired frequencies can be different, and both are needed to fully reconstruct the signal.

For real-valued inputs, however, the paired results will be complex conjugates, where the result (a + i*b) for frequency k is the conjugate of the result for frequency -k: (a - i*b). These results are, in fact, redundant – and we will absolutely take that redundancy into account when implementing the FFT for real-valued inputs.

If you’re more familiar with matrices and linear algebra than complex numbers, the reason this works is that this complex multiply:

gives exactly the same result as multiplying a 2D vector (a, b) with a 2x2 rotation matrix representing the same angle:

Note

It’s worth pointing out something that shows up in the standard formula for the DFT: mathematical notation has an absolutely wild shorthand for $\cos(t) + i*\sin(t)$:

$e^{it}$

That’s right, taking the mathematical constant $e$ to an imaginary power is the same as a rotation around the origin. Why? Well, 3blue1brown has an explanation if you’re curious.

Personally I don’t like this shorthand because it obscures the meaning of the equation, especially for non-mathemeticians (which is, if I’m doing my calculations right, most people).

Calculating the DFT

Okay, we have all of the pieces of how to calculate a DFT, it’s time to put them together.

A DFT takes a set of complex input samples in the time domain and produces the same number of complex output samples in the frequency domain, where the first output sample (at index 0) represents the DC offset, the middle one (at index N/2) represents the Nyquist frequency, and the rest are either the positive or negative angles in between.

There are a couple minor (but important) differences in the value calculation done in the above interactive diagrams vs how the DFT does them:

• Typically the DFT rotates the opposite direction as the above diagrams (i.e. using $\sin(-t)$ instead of $\sin(t)$), so the rotations actually spin the opposite way (clockwise becomes counter-clockwise and vice versa).
• Rather than averaging the results of these rotations, the DFT simply sums them – no division by N is performed.

Otherwise, the process described above can be used to get the contribution of any single frequency within the signal, so then we just need to do that for every possible frequency in the signal.

Thus, here we are, finally at some pseudocode to calculate the DFT:

void CalculateDFT(
  complex[] inputs,
  complex[] outputs,
  int length)
{
  for (outIdx = 0; outIdx < length; outIdx++)
  {
    // Start the sum at 0.
    float sum = 0

    // This is the angular frequency (in radians/sample) for
    // the given output sample (the "test frequency" from the
    // earlier examples).
    float angularFreq  = -2 * pi * outIdx / length

    for (inIdx = 0; inIdx < length; inIdx++)
    {
      // Get the rotation value for the current input index
      // then use that to calculate the resulting rotation.
      float radians = inIdx * angularFreq
      complex rotation = cos(radians) + i*sin(radians)

      // A complex multiply with a unit complex vector is
      // a rotation.
      sum += inputs[inIdx] * rotation
    }

    outputs[outIdx] = sum;
  }
}

Note

As mentioned in a previous note, the standard formulation of the DFT uses an exponential shorthand (which I do not like) to represent the rotation. Thus, the DFT is written in standard mathematical notation as:

$X_k = \displaystyle\sum_{n = 0}^{N - 1} x_n \cdot e^{-i2\pi\frac{k}{N}n}$

…which is equivalent to the above pseudocode.

We’ve finally done it! We have some code that we can run to calculate the DFT for a given signal.

The above algorithm is $O(N^2)$, which may seem like the best that can be done. However, there’s a bit of algorithmic trickery that can be done to get that down to $O(N \log N)$, which is considerably more efficient. That trickery is the fast Fourier transform (FFT), and the next post will go into the details of how that works!