JoshJers' Ramblings

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!