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 (eventually) the specific algorithms that FFT-20XX uses:

1. The Discrete Fourier Transform (← you are here)
2. The Basics of the FFT
3. Some 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.

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

Note

The “periodic signal” bit is important! The algorithm assumes that the block of samples it is processing repeats forever:

However, this won’t be the case for all blocks, such as this one where there’s a hard discontinuity:

This can cause noise in the resulting DFT output. There are a few tricks that can be done for cases like this (like windowing the signal) to effectively “pretend” that it is periodic when it’s not, but that’s out of scope for this post.

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. We’ll get to the inverse DFT in a later post.

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 k (starting at 0 and ending at N - 1) there are k 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 (relative to the x axis). 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 manner (to me, at least): 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, we’ll display the average instead of the final sum (as the sum would be too large to fit in the diagram):

Note that when the frequencies match (signal and test cycles), 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 this sum 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.

Math Time!

If you’re interested in the math of how this works, here you go:

The sum of all samples of a repeating $\cos\,x$ signal is $0$, as it has the same area above and below the $x$ axis, which cancels out in the sum.

Now, if we multiply two such waves together: $\cos\,x \cdot \cos\,y$, this can be rewritten (via trig identities) as:
 $\frac{1}{2}[\cos (x - y) + \cos(x + y)]$

If $x$ and $y$ don’t match, these are just two repeating cosines, and each of them sums to $0$ over all of the samples (just like the single wave above).

However if $x = y$, this becomes:
 $\frac{1}{2}[\cos (x - x) + \cos(x + x)] \implies \frac{1}{2}[\cos\,0 + \cos\,2x]$
 $\implies \frac{1}{2} + \frac{1}{2}\cos\,2x$

That remaining $\cos$ term, over all of the samples, once again sums to 0, which leaves us with a constant $\frac{1}{2}$ in the average.

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 the cosine values for our example frequencies of 1 and 15 cycles:

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 (you need both to represent a wave). 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 cosine 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).

Note

This is just like the optical illusion where rapidly-spinning things like car wheels can appear to be rotating backwards.

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, the frequencies above Nyquist are referred to via their negative indices instead of as the positive aliased ones.

This explains why all of the frequencies in the above example (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, and getting an amplitude-relative result. 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 the resulting average to oscillate! While we get the expected 0.5 average when the phase is 0, it is -0.5 at a phase of 0.5 (a half-wavelength offset) and, even worse, 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 the average is at 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 (8 cycles), 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 (because there is literally no way to represent such a phase shift with the given sample rate).

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 of the test wave at that sample position to get a new 2D position. Those 2D positions are then summed 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 (its length), 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 \cdot 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).

Note

If you, like me, prefer 2D vectors to complex numbers, just pretend that $x + iy$ is the 2D vector $[\,x\,\,\,y\,]$ as you read all of this: the real component is the $x$ coordinate and the imaginary component is the $y$ coordinate.

To be fair, complex numbers do have a nice property that standard 2D vectors do not: they can be multiplied in a meaningful way.

The specific complex multiplications done by the DFT (an input sample [$a + ib$] multiplied by a wave value [$cos\,t + i\,sin\,t$]) ends up rotating the input by the angle $t$!

Note

If you’re familiar with matrices and linear algebra, this complex multiply (with $c$ and $s$ representing the cosine and sine values):

$(x + iy) \cdot (c + is)$

$\implies (xc - ys) \cdot i(xs + yc)$

…results in exactly the same operations as multiplying a 2D vector (x, y) with a 2x2 rotation matrix representing the same angle:

$\begin{bmatrix}x & y\end{bmatrix} \cdot \begin{bmatrix} c & s \\ -s & c \end{bmatrix}$

$\implies\begin{bmatrix} xc - ys && xs + yc \end{bmatrix}$

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 (i.e. 2D) instead of purely real-valued (1D), the contributions to these paired frequencies can be different, and both are needed to fully reconstruct the signal.

For real-valued signals, however, the paired results will be complex conjugates, where the result ($a + ib$) for frequency $k$ is the conjugate of the result for frequency ${-}k$: ($a - ib$); they have the same real (x) term but opposite imaginary (y) terms.

When summed, these values ($a + ib$ and $a - ib$) become $2a$ as the imaginary components cancel out, which is how a signal made up of real $\cos$ terms and imaginary $\sin$ terms can be a purely real-valued signal.

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 finally 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 this process 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

      // `i` here is the imaginary unit
      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 we can use 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!

(thanks to Neintonine, zerostride, and Chuck for giving me some great feedback on this before it went live!)