JoshJers' Ramblings

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!