The DFT by hand and the FFT as its fast algorithm

← Fourier as a bank of correlators

What we do in this lesson

There is no new physics today — the DFT formula is the same as in the previous lesson. The goal of the lesson is to make the formula less abstract and to answer the practical questions that come up when you actually calculate a spectrum in code:

1. We go through the DFT formula by hand on a toy example of four numbers, step by step, with no steps skipped.
2. We find out what the output index $k$ means and in which order np.fft.fft returns the results — this is the first thing that confuses people who work with NumPy.
3. We see why nobody calculates the DFT directly from the formula for long records.
4. We understand the main idea of the FFT — it is not new mathematics, but a method to reuse results that are already calculated instead of calculating from zero.

Back to the correlator formula

In the previous lesson, one correlator tested one frequency: it multiplied the signal samples by the reference rotation and added the results. The DFT is simply a bank of such correlators that work at the same time, one for each frequency of the grid:

Here $n$ is the index of the input sample in time, and $k$ is the index of the output bin in frequency. Not one symbol of the formula changed — we only repeat it for each $k$ from $0$ to $N-1$.

We calculate the DFT by hand on four samples

We take the smallest example in which we can still see each term. Let $N=4$, and let the input signal be a usual real cosine, with exactly one period in four samples:

Note: the signal $x[n]$ itself is only four usual numbers, not vectors on a plane. But the reference signals inside each correlator are still complex — the same as in the previous lesson.

To find $X[k]$ for a specific $k$, multiply each sample $x[n]$ by the related value of the reference rotation $e^{-j2\pi kn/4}$ and add the four results. Here the reference rotation makes an integer number of turns in 4 steps, so its values are the four points on the circle that we already know: $1,\,-j,\,-1,\,j$. We write them for all pairs $(k,n)$ in a table — these are all the reference signals that we need:

The table row with index $k$ is the reference signal of bin $k$. Now, for each $k$, we multiply this row by the input $x=[1,0,-1,0]$ element by element and add. Note: we have $x[1]=x[3]=0$, so these two terms become zero on their own, and we only need to calculate the contribution from $n=0$ and $n=2$:

• $k=0$: $1\cdot1 + 1\cdot(-1) = 0$
• $k=1$: $1\cdot1 + (-1)\cdot(-1) = 1+1=2$
• $k=2$: $1\cdot1 + 1\cdot(-1) = 1-1=0$
• $k=3$: $1\cdot1 + (-1)\cdot(-1) = 1+1=2$

We get $X=[0,\,2,\,0,\,2]$. All four numbers are real, with no imaginary part — this is not a coincidence: the input signal was real. And $2$ is exactly what the rule from the previous lesson predicts for a cosine with amplitude $1$: each of its two components (the clockwise and the counterclockwise rotation) gives the contribution $NA/2 = 4\cdot1/2=2$.

How to read the bin index k as a frequency

The four numbers $X[0..3]$ alone tell nothing about frequencies in hertz — the formula $f_k=kF_s/N$ gives the link to $F_s$. We take $F_s=8$ kS/s, as before:

The last number looks wrong: at $F_s=8$ kS/s, the Nyquist frequency is $\pm4000$ Hz, and $+6000$ Hz is above it. Here the same aliasing effect operates that we studied in Part 1: on a grid of 4 samples, a rotation of $+6000$ Hz is the same as a rotation of $-2000$ Hz. Thus np.fft.fft does not show frequencies above Nyquist — instead, it gives the second half of the bins negative frequency labels:

This is the order in which np.fft.fft actually returns the array: first zero, then the positive frequencies in increasing order, then the negative frequencies — also in increasing order, starting from the Nyquist frequency. The Nyquist frequency itself ($k=2$) has the label $-4000$ Hz by convention, but it could equally be $+4000$ Hz — it is the same rotation.

This order is not convenient for a plot — we want the frequencies to increase from left to right. This is what np.fft.fftshift() does: it changes the order of the array and calculates nothing again:

Now we see the same as in the previous lesson: a cosine gives two equal peaks — at $+2000$ and at $-2000$ Hz. A correlator cannot tell a positive rotation from a negative one while the signal is real, so the energy always divides equally between the mirror frequencies.

Why we need a faster algorithm

For each of the $N$ bins, the DFT formula must go through all $N$ samples — this is approximately $N^2$ multiplications and additions. For our example with $N=4$, this is $16$ operations, a very small number. But imagine a receiver that continuously makes the spectrum of an audio stream at $48\,000$ samples per second and updates the picture many times per second — at this scale, $N^2$ quickly becomes too large. For a record of length $N=1\,048\,576$ ($2^{20}$, approximately 1 second on some SDRs), the number of operations is $N^2\approx1.1\times10^{12}$ — more than a trillion for one single spectrum.

The FFT is a fast algorithm that calculates the same sum as the DFT, but not directly. It is not a different transform and not an approximation: the FFT returns the same numbers as the DFT formula, but it gets to them by a different path.

The idea of the FFT: we use what is already calculated

The key observation is simple: if you divide a signal of length $N$ into two halves — the even and the odd samples — then each half is itself a signal of length $N/2$, and you can calculate its DFT separately. When it is calculated, both outputs of the full DFT — $X[k]$ and its “pair” $X[k+N/2]$ — come from the same two half-length results, with a different sign on one term. Thus, instead of going through all $N$ samples again for each of the $N$ bins, we do the work only once on each half of the data and use it two times.

We write this as a formula. We write the DFT of the even samples as $E[k]$ and the DFT of the odd samples as $O[k]$ (both have length $N/2$):

Then you can make the full DFT of length $N$ from them like this:

This pair of equations is called a “butterfly” — both outputs come from the same pair of numbers $E[k]$, $O[k]$, only with a plus and a minus.

We check this on our example $x=[1,0,-1,0]$ to see that the numbers do not change. The even samples are $[1,-1]$, the odd samples are $[0,0]$. For two numbers, we calculate the DFT with the same logic as above, but with $N/2=2$:

• $E[0]=1+(-1)=0$, $\ E[1]=1\cdot1+(-1)\cdot(-1)=2$
• $O[0]=O[1]=0$ — the odd samples are zero, so their DFT is also zero

We make the full result with the butterfly formula:

• $X[0]=E[0]+O[0]=0$, $\quad X[2]=E[0]-O[0]=0$
• $X[1]=E[1]+e^{-j\pi/2}O[1]=2$, $\quad X[3]=E[1]-e^{-j\pi/2}O[1]=2$

Again $X=[0,2,0,2]$ — the same as the result of the direct calculation. But instead of $16$ multiplications, it was sufficient to calculate two DFTs of two samples each and combine them with a pair of additions.

On such a small example, you cannot see the saving, but the method scales recursively: you can divide each of the two halves of length $N/2$ in half again in the same way, and so on, until only pieces of one sample remain (the DFT of one number is the number itself). There are $\log_2 N$ levels of division, and at each level the butterflies do approximately $N$ additions and multiplications in total. The total cost is $O(N\log N)$ instead of $O(N^2)$. For $N=1\,048\,576$, this is $N\log_2 N\approx2\times10^7$ operations — approximately $50\,000$ times fewer than $N^2$.

Lab: we compare the DFT from the definition and the FFT

The script first calculates the example of this lesson again and compares it with np.fft.fft. Then it measures the time of the same operation with two methods at different $N$. Before you run it, it is useful to predict: by how many times must the time of the “manual” method increase when $N$ doubles, and by how many times for the FFT?

uv run labs/03/02-dft-manual-fft.py --show

The left plot repeats the table above after fftshift: two peaks of $2$ at $\pm2000$ Hz. The right plot is the run time as a function of $N$ on logarithmic axes. The “manual” method is a direct multiplication by the $N\times N$ matrix of reference signals — the same DFT formula, with no division into even and odd samples. Its time increases as $N^2$ and makes a straight line with a clear slope on the plot. The time of np.fft.fft increases as $N\log N$ — an almost horizontal line in the same range.

At $N=4096$, the manual method takes several thousand times more time than np.fft.fft, but both return the same numbers within rounding error ($\sim10^{-16}$). In the rest of the course, we will always calculate the spectrum with np.fft.fft, and use the DFT formula only for explanations and checks on small examples.

Self-check