Fast Fourier Transforms

Mathematical / Symbolic terms

Divide and Conquer

The FFT works because some folks figured out you can split a DFT into two DFTs half the size, and then DFT them together, and it works! What does it mean to 'split a DFT' or 'DFT them together'? It's pretty hard to get a genuine answer to that without learning what the parts that a DFT uses to do its job are.

Decimation In Frequency vs Decimation In Time

We are doing a fourier transform - so taking a sequence of values over time, and using that to develop a sequence of values indicating the strength of various frequencies in the signal. This means we take a 'time domain' signal, and use it to make a 'frequency domain' signal. When we split up our data to do a FFT, we can split either the output Frequency domain data or the input time domain data. This gets us the two main types of FFT: Decimation in Frequency (DIF) and Decimation in Time (DIT). In this example, where not otherwise stated, I will describe the DIT variant. Cooley and Tukey, the guys who did the original FFT, described a 'radix-2 DIT' FFT. The DIF one is kind of like a backwards version of DIT (but it's not an iFFT)!

'Decimation' is another way of saying 'factorisation' - i.e. finding the common factors. In our case, we are factorising the transforms that make up the DFT, pulling out all the repeated calculations we can. These transforms are probably best represented as matrices, so we are factorising matrices. This gets us 'sparce matrices' - matrices who have very few non-zero values (sometimes along a diagonal, sometimes in little blocks, etc); the minimum amount of information needed for a 2D transform. This process happens through the reordering and twiddling of the data. The data is re-organised because the order of 'time' is now irrelevant as we look for periodic patterns across time. Across, here, is very literal - the data from different time points has to 'cross over' the ones in-between to be used for a minimal calculation to get the combination and difference of those points. Those sums and differences are the core of the FFT, and they are called 'butterflies'.

The Butterfly

2-point - smallest

Twiddle the odd value first. Then sum both values and put the answer in the even slot, and take their difference and put it in the odd slot. All other butterflies can be composed of this basic structure; it is the core of our correlation-detection.

yk0 ------------>    |+|-> yk0 + yk1 * twiddle
				  \ /
				   X
				  / \
yk1  -->  * ω -->    |-|-> yk0 - yk1 * twiddle

When an FFT is described by its radix (e.g. "radix-16 DIF FFT"), this tells us the smallest butterfly found in the algorithm. HOWEVER! This will still usually be composed of pairwise calculations, even in highly parallel implementations. But then, if all these pairwise operations simply merge two values, how do we merge values from across the whole input sequence- as we were saying, to treat time as irrelevant to find the similarities and differences across time. Well, other than going across the butterfly, the data also goes across the set of sub-transforms using stride.

Stride

How do we compose a larger butterfly of 2-point calculations? Using symmetry.

We assign the second variable to the 'other side of the circle' - the opposite string/slice of cake. Remember, data in the FFT goes across.

We can do this by considering some of our symbolic notation (maths letters) from the start. We have this quantity of nth-size portions, called m. E.g., in a cake of n=6, when m=2, we cut three thirds A[k + j + m/2]

When combined with the butterfly, this results in the following movement of information:

stage 1) a0,a2 = (x0,x1), a1,a3 = (x2,x3)
stage 2) a0,a2 = (a0,a1), a1,a3 = (a2,a3)

Twiddle factors

Phase

Detecting aliases of higher frequencies than the base case (cos(0) and sin(0)) requires a phase 'kickback' for every increase in octave. You really need to see the graphs for this one - Mark Newman is the GOAT https://www.dsprelated.com/showarticle/1716.php The simple description is: the more wobbly you want to check for, the more phase offset you just apply! The larger a DFT, the more frequencies it can check for: 2 observations isn't much to establish a trend from. 4 is okay for simple curves. 8 helps see if there's something sine-shaped in the data you've collected. 65536 observations per second might be enough to model sound waves for human hearing! So, following; the larger a DFT, the bigger the 'twist' applied to rotate it back through the cycle of the (co)sine wave. Why 'twist' and not 'offset'?

Picture cosine weaving along the timeline. Now place sine face-down on the third spacial axis; the complex plane is perpendicular to the timeline. If you know euler's identity, you can picture a spiral that those curves project onto their plane, turning through time. If we want to slide these curves along the axis, we can imagine how if you pulled this great long drill-bit, it would naturally turn in a spiral. You have to twist back to arrange the original function such that it looks just as it did at the start, but now 0 on the time-axis is different. Phase adjustments are 'twisty' like this; it naturally digresses along a complex exponential, and requires work to displace the whole curve across that phase difference. So, this is just a little bit of why all this trigonometry is everywhere (cosine being a wave is part of something bigger). And if you don't like trigonometry, that's okay! The FFT favours matrices, which can be arranged with complex terms, or other ways to represent orthogonality. There are, however, some concepts that unite the different approaches to representing the transformations, such as roots of unity.

Roots of Unity

If your brain works like mine, you will find it helpful that these phase kickbacks line up with the 'roots of unity' - complex numbers that correspond with an index to which you can power a complex number (z) to get 1 (z^r = 1 where r is the root of unity) they line up at fractional positions on the unit circle in the complex plane - i.e. they are at positions you slice when cutting a cake evenly if you have things arranged evenly in a circle, you can use roots of unity of the number of things to split them evenly e.g. If you are weaving two portions of hair into a braid, the only root of unity is opposition - two opposite parts. If you have 4, you can split in two sets of two- x and y, cartesian axes, four points on a compass; 90 degrees, i. The root from -1 (the next split after opposition) If you have 6, you can split into two sets of three , or one set of two. The common denominator, 6, tells us if we pulled our gracious braidee's hair taut and arranged their hair around a 6-spoked cog, we could evenly make the required divisions. Roots of unity help us divide up multidimensional things into the core parts of the transformation by fractional amounts. Remember we were trying to factorise matrices?:)

Twiddling in Practice

We work our way through the twiddle factors in a sub-transform/DFT through repeated application of:

ω = 1
for j in m/2:
	butterfly
	ωm ← exp(−2πi/m)

As the sub-transforms get larger, the inner loop runs more frequently before the outer loop comes round, resets ω, and changes the values of k, j, and m. Combined with the stride, this does most of the 'magic' of the FFT that isn't obvious at first glance.

Order - Bit Shuffling and more

You may notice that repeatedly doing the above butterflies and strides ends up with a considerably jumbled order. In fact, it has the same end result as if you had first split your initial pile of data into a stack of odd-indexed items and a stack of even-indexed ones. Most of the time that we do an FFT, we need to do this to get our frequency data output in a low to high order. Now here's a neat trick; reversing the order of the digits of the indices does exactly this sort.

000	- 0		000 - 0
001	- 1		100 - 4
010	- 2		010 - 2
011	- 3		110 - 6
100	- 4		001 - 1
101	- 5		101 - 5
110	- 6		011 - 3
111	- 7		111 - 7

Putting it Together: Wikipedia pseudocode

Okay, ready to see how it looks all put together? I've added some comments to this pseudocode from the wikipedia page for Cooley-Tukey to make it easier to follow:

input: Array a, of n complex values, where n is a power of 2.
output: Array A, the DFT of a.

n ← A.length 							// get n
A ← bit-reverse-copy(a) 				// shuffle via bit-reverse

// stages:
for s = 1 to log(n) do
    m ← 2^s 							// calculate m based on the stage
    ωm ← exp(−2πi/m) 					// calculate the factor used to work out twiddles (m is the only variable it depends on)
    // sub-transform:
    for k = 0 to n-1 by m do 			// note here, k goes from 0 to n-1 in steps of m. e.g. when m=2 and n=8: k = 0, then 2, then 4, then 6
        ω ← 1 							// reset twiddle to 1 at the start of every sub-transform
        // sub-sub-transform:
        for j = 0 to m/2 – 1 do
            t ← ω A[k + j + m/2] 		// twiddle one value
            u ← A[k + j] 				// get another value
            // butterfly:
            A[k + j] ← u + t 			// sum values together
            A[k + j + m/2] ← u – t 		// find the difference of the values
            ω ← ω ωm 					// modify the twiddle factor for the next round

return A

Variants on the FFT

Radix (2,3,4,6,8,16, split...)

You can start your FFT with a bigger transform than 2 if you want. But you should really keep it to one of the products of twos (2^n), since everything is square in this algorithm (geometrically, it's full of orthogonal pairs). What would you do with the odd value out? Just don't try to use a prime number; you'll have a really hard time splitting things into even parts.

Stride II: Attack of the Clones (Stockham)

Okay the stride is all well and good but what's this bit reversal nonsense? Seems like a waste of time. It absolutely is a waste of time, too, because we already move all the values around in the algorithm. Couldn't we just do that a little better, without breaking it? Yes! If we use the Stockham FFT. As you might have guessed, it uses the method the FFT already moves things about with- stride.

First; for the twiddle, we look up A[j + r * stride + k * L] Here L = logR(n), where R is the radix (smallest transform). So in this radix-16 example, L = log16(n), all the way through. We do this 16 times per invocation for the radix-16, so we use a for loop to step through them as values of r.

Then, once we have our 16 twiddles calculated, we can do the butterfly:

acc ← (0,0)
for m = 0 to 16 do
	acc ← acc + cmul(x[m], W16[(r * m) & 15u])
y[r] = acc;

Lastly, Stockham requires that we reshuffle the data after each round of processing:

base ← k * (L * 16);
for r = 0 to 16 do
    data_out[j + r * L + base] ← y[r];

Great, now we have our FFT with no time wasted reshuffling bits around! Just addition, subtraction, multiplication, and a little trig.

Happy transforming:)