Fourier Transform and the FFT

Time, frequency, and the basis.

On this page

1 A signal and its spectrum 2 The basis waves 3 The DFT 4 Magnitude and phase 5 Which bin is which frequency 6 Frequency resolution 7 Back to time 8 Why the FFT is fast 9 Convolution becomes multiplication 10 In practice 11 When a spectrum looks wrong 12 References

The Fourier transform turns a list of samples into a list of frequencies: for each one, how much of it the signal holds and where its wave starts. The fast Fourier transform (FFT) computes it once per clip, or once per block of samples for a long signal, and the inverse turns the frequencies back into the same samples.

1 A signal and its spectrum

A recording is a list of samples x[n], where n counts the samples and the sample rate f_s is how many arrive each second. The FFT of N samples gives one value per frequency, called a bin; X[k] is the value of bin k. Each bin is a complex number, a size and an angle: how much of that frequency the signal holds, and where its wave starts (Section 4). The size is written |X[k]|. Figure 1 plots the sizes scaled as 2|X[k]|/N, so each tone reads its own amplitude.

Figure 1. 1100 ms of x[n]. 250 Hz reads 0.7, 3120 Hz 1.0. 4Axis ends at f_s/2 = 500 Hz.

The FFT also lists negative frequencies, waves that turn the other way round the circle. For a real signal they mirror the positive ones, so a tone at frequency f puts half its size at +f and half at −f. np.fft.rfft keeps only 0 to f_s/2, N/2 + 1 = 501 bins, the one-sided spectrum, and the plot adds the two halves back by doubling each size, 2|X[k]|/N. Here f_s = 1,000 Hz and N = 1,000, so bin k is k Hz; Section 5 shows why.

Figure 2. 1Each tone, half at −f, half at +f. 2rfft keeps 0 to f_s/2, halves added.

The NumPy code that draws Figure 1.

import numpy as np
import matplotlib.pyplot as plt

fs, N = 1000, 1000                       # 1 s at 1,000 samples per second
t = np.arange(N) / fs
tones = 0.7 * np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t)
x = tones + 0.2 * np.random.default_rng(0).standard_normal(N)
X = np.fft.rfft(x)                       # bins 0 to N/2
f_k = np.fft.rfftfreq(N, d=1 / fs)       # k * fs / N, in Hz
fig, (top, bottom) = plt.subplots(2, 1)
top.plot(t[:100], x[:100])
bottom.plot(f_k, 2 * np.abs(X) / N)
bottom.set_xlabel("frequency (Hz)")
plt.show()

1Makes the signal of Figure 1, marker 1: two sines and a little noise, f_s = 1,000 Hz and N = 1,000.
2Takes the FFT, keeping bins 0 to N/2 because the rest mirror them, and the frequency k f_s/N of each bin.
3Draws the time view and 2|X[k]|/N, which reads each tone’s amplitude, markers 2 to 4.

2 The basis waves

How do you find which tones are inside? Compare the signal with a set of basis waves, one at a time, and measure how much of each one it contains. For N samples there are N basis waves.

Figure 3. Wave k makes k turns in N = 8 samples (dots), as 1a cosine and 2a sine.

Each basis wave has a cosine part and a sine part, because a tone can start at any point of its cycle. One complex number holds both: a real part plus a part times j, the square root of −1. The complex number written e^jθ is the point at angle θ on a circle of radius 1.

Figure 4. 1The red dot e^jθ at angle θ. 2Its horizontal shadow is cos θ. 3Its height is sin θ.

One number that holds a cosine and a sine.

e^jθ = cos θ + j sin θ

j is the square root of −1. As θ grows, e^jθ walks round the unit circle; its real part, cos θ, is the dot's shadow on the horizontal axis and its j part, sin θ, is its height. One complex number therefore carries a cosine and a sine at once. It builds on the circle of Figure 4, and the DFT card in Section 3 uses it for every basis wave.

This is a change of basis, the same data described along new axes. PCA learns its axes from the data; the Fourier basis is the same fixed set of waves for every signal.

3 The DFT

The discrete Fourier transform (DFT) measures this. Four samples of a 4 kHz tone taken at 16 kHz, x = [1, 0, −1, 0], give X = [0, 2, 0, 2]: all of it sits in bin 1 and its mirror, bin 3.

The DFT and its inverse.

X[k] = Σ_n x[n] e^−j2πkn/N

Basis wave k is e^−j2πkn/N, k full turns across the N samples. The sum runs over the samples, n = 0 to N − 1. Multiplying x[n] by the wave sample by sample and adding up, a dot product, gives a large total only when x turns at the same rate, so X[k] measures how much of wave k is in x, with its cosine and sine parts together.

x[n] = 1N Σ_k X[k] e^j2πkn/N

The inverse adds the N basis waves back, k = 0 to N − 1, each weighted by its amount X[k], with the exponent's sign flipped so each wave turns the other way. The 1N undoes the N-fold total that the forward sum built up. Both lines build on the basis waves of Section 2, and every later section reads this X[k].

The DFT as one matrix product. Row k of F is basis wave k, and the last line checks the product against np.fft.fft, which does not scale (np.fft.ifft divides by N).

import numpy as np

N = 8
n = np.arange(N)
k = n[:, None]                          # bins as a column
F = np.exp(-2j * np.pi * k * n / N)     # shape (N, N)
x = np.random.default_rng(0).standard_normal(N)
print(np.allclose(F @ x, np.fft.fft(x)))

1Row k of F is basis wave k, so F @ x is the DFT (eq card in Section 3, line 1).
2Checks the matrix product against np.fft.fft on random input.

Basis wave k = 5

Figure 5. Tones in bins 5 and 12. 1Signal times cosine. 2Bin k, large on a match.

4 Magnitude and phase

Each X[k] is a complex number, so it has a size and an angle.

Figure 6. Same magnitudes. 1All phases 0 add to a click. 2Random phases add to noise.

The two parts of each DFT value.

|X[k]|

|X[k]|, the magnitude, is the length of the complex number, how much of wave k there is; doubling a tone doubles it.

∠X[k]

∠X[k], the phase, is its angle, where in its cycle wave k starts; delaying the signal turns this angle and leaves the magnitude alone. Both build on the DFT card of Section 3, and Section 7 needs both to go back to time.

Because phase sets timing, How FIR and IIR Filters Shape Magnitude and Phase tracks what a filter does to it, and Spectral Subtraction for Denoising edits the magnitude and keeps the noisy phase.

5 Which bin is which frequency

A 32 ms speech frame at f_s = 16 kHz holds N = 512 samples. Its bins sit 31.25 Hz apart, so 1 kHz lands in bin 32 and the Nyquist limit, 8 kHz, in bin 256.

The frequency of bin k.

f_k = k f_sN

Bin k makes k turns in N samples, and N samples last Nf_s seconds, so its frequency is k f_sN. Bins sit f_sN apart, so more samples give finer spacing, and bin N2 lands on the Nyquist frequency, f_s2. It builds on the DFT card of Section 3, and Section 6 asks when two tones get bins of their own.

For a real signal, bins above N/2 mirror the ones below, the same fold as in Aliasing and Anti-Aliasing Filters.

6 Frequency resolution

Two tones 4 Hz apart that start in phase, at 100 and 104 Hz, show as one peak in a 0.15 s record and as two in a 0.20 s one.

How finely a record of N samples separates frequencies.

Δf = f_sN

Δf is the bin spacing, the sample rate f_s split into N equal steps. At f_s = 1,000 Hz, N = 1,000 gives 1 Hz bins and N = 100 gives 10 Hz bins. Cut plainly, with no fade at the record’s edges, two tones that start in phase show as two peaks once they are about 0.65Δf apart; other starting phases can need up to about 1.3Δf. The in-phase 4 Hz pair of Figure 7 is one peak at T = 0.15 s, where 4 Hz is 0.6Δf, and two peaks at 0.20 s, where it is 0.8Δf.

Δf = 1T

N samples last T = Nf_s seconds, so Δf is one over the record length, the frequency resolution: only a longer recording separates closer tones. Zero-padding appends zeros before the FFT, so bins sit closer than Δf and draw more points on the same curve, but T, the length of real data, does not change, so close peaks do not separate.

Record length T = 0.15 s

Figure 7. 1Two tones (dashed). 2Spectrum of T seconds; dots are bins.

7 Back to time

In the signal of Figure 1, bin 50 alone brings back the 50 Hz sine and bin 120 the 120 Hz wave (Figure 8); all 501 bins together give every sample back.

A real signal rebuilt from its rfft bins.

x[n] = Σ_k A_k cos(ω_kn+φ_k)

Every sample is a sum of cosines, one per bin from 0 to N2. Bin k adds a cosine of height A_k that turns ω_k radians per sample and starts at angle φ_k. All N2 + 1 bins, 501 for the signal of Figure 1, give the samples back exactly; leaving a bin out removes that frequency from the result.

ω_k = 2πkN

Bin k makes k full turns in N samples, the basis wave of Section 2. In hertz it is f_k = k f_sN.

A_k = 2|X[k]|N

The height of wave k. The forward sum grows a tone N2 times, and the other half of its size sits in the mirrored negative bin that rfft drops, so 2N undoes both. Bins 0 and N2 have no mirror and use |X[k]|N.

φ_k = ∠X[k]

Where wave k starts, the angle of the complex bin. A sine starts a quarter turn after a cosine, so both sines come back with φ near −90°. Figure 8 reads −90° for bin 50 and −89° for bin 120, where the added noise nudges the angle.

Figure 8. 1Bin 50, a cosine 0.70 at −90°. 2Bin 120. 3Their sum fits the gray signal.

Keep every magnitude and scramble the phases, and the same cosines peak at different moments and add to a different signal, the noise of Figure 6 instead of its click. In Figure 9 each target repeats every 8 ms, a base frequency of 125 Hz, so wave k sits at 125k Hz. A square needs only odd harmonics, waves at odd multiples of its base; its even waves have size zero, so they count in K but add nothing.

Waves used K = 3

Figure 9. 1The sum of K waves. 2Each wave’s size, filled if used.

The round trip in NumPy, with the signal of Figure 1.

import numpy as np

fs, N = 1000, 1000
t = np.arange(N) / fs
rng = np.random.default_rng(0)
x = 0.7 * np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t) + 0.2 * rng.standard_normal(N)
X = np.fft.rfft(x)                       # 501 bins
x_back = np.fft.irfft(X, n=N)            # n=N keeps the original length
print(np.allclose(x, x_back), np.max(np.abs(x - x_back)))   # True 8.9e-16

1Makes the signal of Figure 1.
2Goes to the 501 bins (eq card in Section 3, line 1).
3Adds the waves back (eq card in Section 7, line 1); n=N keeps the length when N is odd.
4Prints True and the largest error, near machine precision.

8 Why the FFT is fast

Done directly, the DFT costs N multiply-adds for each of N bins. The fast Fourier transform (FFT) gets the same numbers with far less work. It is not an approximation of the DFT but exactly the DFT, regrouped (Anushiravani 2016). It splits the samples into even and odd, does two half-size DFTs and combines them. Each half splits again, until only pairs are left.

A butterfly joins two half DFTs.

X[k] = Σ_n x[n] W^kn

The DFT the FFT computes, with W = e^−j2π/N. Each of the N bins sums all N samples; the splits regroup this sum into even and odd samples.

X[k] = E[k] + W^k O[k]

E and O are the DFTs of the even and odd samples. The odd ones start one sample late, so the number W^k = e^−j2πk/N turns O[k] to line up. Adding gives the first half of X.

X[k+N2] = E[k] − W^kO[k]

Bin k + N2 flips the sign of the turned odd part, so the same E, O and W^k give the second half for one multiply per pair. It builds on the DFT card of Section 3, and Figure 11 counts the operations its stages need.

Figure 10. 1Even and odd splits down to pairs. 2Butterflies join them into 8 bins.

Each split halves the size, so there are log₂ N stages (log₂ N is how many times you can halve N before reaching 1), and each stage costs about N operations. At N = 1,024 that is 10 stages of about 1,024 operations each.

Operations for a direct DFT and for an FFT.

N² → N log₂ N

The direct DFT does N operations, each a multiply-add, for each of N bins. The FFT halves the size log₂ N times and each stage costs about N; Figure 11 plots both counts on log scales.

Figure 11. At N = 1,024, 1the DFT needs 1,048,576 operations, 2the FFT 10,240.

James Cooley and John Tukey published this algorithm in 1965; their version needs fewer than 2N log₂ N operations. Every np.fft call uses an FFT.

9 Convolution becomes multiplication

A filter’s sliding weighted sum, the convolution x ∗ h, is also what a convolution layer computes (CNN Basics: Convolution, Stride and Pooling). In frequency the same operation is much simpler, a bin-by-bin product of the two spectra.

Figure 12. 1In time, 3 weights slide along the input. 2In frequency, X times H gives Y.

Filtering in time and in frequency.

x ∗ h ↔ X · H

The ∗ is convolution, the weighted sum that slides along the signal, and ↔ links a signal to its spectrum. A filter only scales and shifts each basis wave, so in frequency the whole slide becomes one multiply per bin, X times H. It builds on the convolution of LTI Systems and Convolution, and How FIR and IIR Filters Shape Magnitude and Phase computes H for real filters.

The DFT’s waves repeat every N samples, so the product gives a circular convolution, where the end of the output wraps onto the start. Zero-pad both signals to the full output length first. A convolution reverb, which adds a room’s echo to a dry recording, does exactly this; Figure 13 runs a tiny case.

Figure 13. [1, 2, 3] through [1, 0.5, 0.25]. 1True. 2Wrapped. 3Padded.

10 In practice

The FFT runs inside every spectrum, STFT, PSD and fast convolution later in the series, once per block of samples.

Step	How often, on what data	State it keeps	Cost
DFT, direct	once per block of samples, on one block	none	block length squared multiplies
FFT	once per block of samples, on one block	none	block length times its log₂

Setting	Default	Source
norm	"backward", forward transform unscaled, inverse divided by `N`	`numpy.fft.rfft`
`rfft` output	`N`//2 + 1 bins for `N` real samples	`numpy.fft.rfft`

Knob, where it enters	Raise it	Lower it	Default
FFT size `N`, the bin spacing `f_s`/`N`	Finer bins and more cost (about `N` log₂ `N`); zero padding only interpolates	Coarser bins; `N` = 1 leaves only the mean	Every sample, `N` = `len(x)` (`numpy.fft.rfft`)
Zero-padding, `n` in `np.fft.rfft(x, n)`	More points drawn on the same curve, more cost; close tones still merge	`n` = `N` plots only the bins, `Δf` apart	No padding, `n` = `len(x)` (numpy.fft.rfft)

NumPy documents numpy.fft.rfft and its scaling.

10.1 Gotchas

Magnitudes grow with N; divide by N/2 to read a tone’s amplitude (Section 7).
rfft drops the mirror half, so power sums double bins 1 to N/2 − 1.

11 When a spectrum looks wrong

What you see	Likely cause	What to do
Peaks at frequencies that make no sense	Reading bin `k` as `k` Hz	Bin `k` is `k` `f_s`/`N`; use `rfftfreq`
The rebuilt signal is `N` times too big	The 1/`N` went missing	`numpy` puts it on the inverse; keep fft and ifft as a pair
A mirror image in the top half	fft on a real signal	Use `rfft` and keep bins 0 to `N`/2
One tone spreads over many bins	The tone falls between bins (leakage)	Use a window; see The STFT and Spectrograms
Filtered output wraps around at the edges	Circular convolution	Zero-pad to the full output length

How do you keep some frequencies and drop others? That is How FIR and IIR Filters Shape Magnitude and Phase.

12 References

J. W. Cooley and J. W. Tukey, "An Algorithm for the Machine Calculation of Complex Fourier Series," Mathematics of Computation 19(90), pp. 297 to 301, 1965.
A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson, 2010, chapters 8 (the discrete Fourier transform) and 9 (computation of the DFT).
S. W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing, chapters 8 (the DFT), 10 (Fourier transform properties), 12 (the FFT) and 18 (FFT convolution).
NumPy documentation, numpy.fft (normalization, rfft, rfftfreq).
Ramin Anushiravani, Example-Based Audio Editing, M.S. thesis, University of Illinois at Urbana-Champaign, 2016, sections 2.6 (the FFT) and 2.7 (circular convolution).