The 2D DFT & FFT in Practice

"I do half the work and hand the other half to my twin, who does half the work and hands the other half to her twin. Somehow we finish in log time and management takes the credit."

An Overworked FFT Butterfly

Section 4.1 argued that images are sums of waves; this section hands you the machine that finds the weights. The good news is that the machine is two function calls: one forward, one inverse. The work is in understanding its conventions, because nearly every confusing spectrum plot, mysterious wraparound artifact, and "why is my filtered image complex?" bug traces back to a convention that was never stated. We state them all.

1. The 2D Discrete Fourier Transform, Defined Intermediate

For an $M \times N$ image $f(x, y)$, the 2D DFT produces an $M \times N$ array of complex coefficients:

$$F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y)\; e^{-j 2\pi \left( \frac{ux}{M} + \frac{vy}{N} \right)}$$

and the inverse transform reconstructs the image exactly:

$$f(x, y) = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v)\; e^{+j 2\pi \left( \frac{ux}{M} + \frac{vy}{N} \right)}$$

Each coefficient $F(u,v)$ is the inner product of the image with one complex sinusoid, exactly the gratings of Section 4.1 in complex form. Three structural facts follow directly from the definition, and all three matter daily in practice.

First, $F(0, 0) = \sum_{x,y} f(x,y)$: the zero-frequency or DC coefficient is the plain sum of all pixels, so $F(0,0)/MN$ is the image mean. Second, for a real-valued image the spectrum has conjugate symmetry, $F(-u, -v) = F^{*}(u, v)$: half the coefficients are determined by the other half, which is why the spectrum dots in Figure 4.1.1 came in mirrored pairs, and why rfft2 can store just over half the spectrum with no information loss. Third, the DFT treats the image as periodic: it analyzes an infinite tiling of your image, not the image alone. The right edge wraps to meet the left edge, and if those edges differ in brightness, the artificial seam injects a bright horizontal-plus-vertical cross into the spectrum. The same periodicity is why unpadded frequency-domain convolution wraps around the borders, the trap discussed in part 4 of this section.

2. Computing and Displaying a Spectrum Beginner

The raw output of np.fft.fft2 stores frequency $(0,0)$ at array index $[0, 0]$, the top-left corner, with positive frequencies counting up and negative frequencies stored in the upper half of each axis (an inheritance from the transform's periodicity). For viewing, everyone re-centers the spectrum with fftshift so DC sits in the middle, and compresses its enormous dynamic range with a log. The full display recipe:

import numpy as np
import matplotlib.pyplot as plt
from skimage import data

img = data.camera().astype(np.float64)         # 512 x 512 grayscale

F = np.fft.fft2(img)                           # complex128, DC at [0, 0]
F_shifted = np.fft.fftshift(F)                 # DC moved to the center
magnitude = np.log1p(np.abs(F_shifted))        # log(1 + |F|) tames the range

# Sanity checks worth internalizing
print(F.shape, F.dtype)                        # (512, 512) complex128
print(np.allclose(F[0, 0].real / img.size, img.mean()))   # True: DC = mean
recon = np.real(np.fft.ifft2(F))
print(np.allclose(recon, img))                 # True: perfectly invertible

plt.imshow(magnitude, cmap="gray"); plt.axis("off"); plt.show()

Without the log, the DC term and its low-frequency neighbors are so much larger than everything else (often by six orders of magnitude) that the display would be a single white pixel on black. Without the shift, the visually meaningful "low frequencies in the middle, high frequencies at the edge" layout is scrambled into four corner fragments. Figure 4.2.1 shows exactly what fftshift rearranges; commit it to memory, because forgetting the shift (or applying it twice, or filtering a shifted spectrum with an unshifted mask) is the most common frequency-domain bug in student code.

3. The FFT and the Library Landscape Intermediate

Evaluating the DFT definition directly costs $O(M^2 N^2)$: every output coefficient touches every pixel. The fast Fourier transform computes the identical result in $O(MN \log MN)$ by recursively splitting each 1D transform into even-indexed and odd-indexed halves and reusing the shared work (the "butterfly" pattern). A 2D transform is just 1D FFTs over all rows, then all columns. At $4096 \times 4096$ the speedup factor is roughly a million; the FFT is the reason frequency-domain methods are practical at all, and it is routinely listed among the most important algorithms of the twentieth century.

In the Python ecosystem you will meet four implementations, and it is worth knowing when each earns its keep. numpy.fft is always there and always correct; its rfft2 exploits conjugate symmetry to do roughly half the work on real images. scipy.fft is a faster drop-in with a workers= argument for multi-core transforms and better support for odd sizes. OpenCV's cv2.dft is fastest when you feed it sizes it likes; cv2.getOptimalDFTSize tells you the next "friendly" size (products of 2, 3, and 5) to pad to, a habit inherited from the array-centric stack of Chapter 0 that can double throughput for awkward dimensions. And torch.fft runs on the GPU and, crucially, is differentiable, so a frequency-domain term can sit inside a training loss.

import cv2
import numpy as np
import torch

img32 = img.astype(np.float32)
h, w = img32.shape

# --- OpenCV: pad to an FFT-friendly size, then transform -------------
H_opt, W_opt = cv2.getOptimalDFTSize(h), cv2.getOptimalDFTSize(w)
padded = cv2.copyMakeBorder(img32, 0, H_opt - h, 0, W_opt - w,
                            cv2.BORDER_CONSTANT, value=0)
dft = cv2.dft(padded, flags=cv2.DFT_COMPLEX_OUTPUT)   # (H_opt, W_opt, 2)
mag_cv = cv2.magnitude(dft[..., 0], dft[..., 1])

# --- PyTorch: half-spectrum of a real input, GPU and autograd ready ---
x = torch.from_numpy(img)                  # add .cuda() if a GPU is present
F_half = torch.fft.rfft2(x)                # shape (512, 257) complex
power = F_half.real**2 + F_half.imag**2    # differentiable spectral power

4. The Convolution Theorem in Practice Intermediate

The single most consequential theorem of this chapter connects it back to Chapter 3: convolution in the spatial domain equals element-wise multiplication in the frequency domain,

$$f * h \;\longleftrightarrow\; F \cdot H$$

Sliding a kernel across an image, the operation you implemented with quadruple loops in Chapter 3, is the same as transforming both image and kernel, multiplying the spectra point by point, and transforming back. For a $K \times K$ kernel on an $M \times N$ image, direct convolution costs $O(MNK^2)$ while the FFT route costs $O(MN \log MN)$ regardless of kernel size. Small kernels are cheaper directly; large kernels are dramatically cheaper through the FFT. The crossover on typical hardware arrives surprisingly early, around kernel widths of 11 to 15:

from time import perf_counter
from scipy.signal import convolve2d, fftconvolve

big = np.random.rand(1024, 1024)

for k in [3, 7, 15, 31, 63]:
    kernel = np.random.rand(k, k)
    t0 = perf_counter(); convolve2d(big, kernel, mode="same"); t1 = perf_counter()
    fftconvolve(big, kernel, mode="same");                     t2 = perf_counter()
    print(f"{k:2d}x{k:<2d}  direct {t1-t0:7.2f} s   fft {t2-t1:6.3f} s")

# Representative output (one laptop run; your numbers will differ,
# the shape of the trend will not):
#  3x3   direct    0.07 s   fft  0.090 s
#  7x7   direct    0.30 s   fft  0.092 s
# 15x15  direct    1.30 s   fft  0.095 s
# 31x31  direct    5.50 s   fft  0.098 s
# 63x63  direct   22.10 s   fft  0.103 s

out = fftconvolve(big, kernel, mode="same")     # padding, cropping handled
# scipy.signal.oaconvolve: same API, overlap-add, better for huge images

This theorem is also a lens for what comes later in the book: the learned convolution layers of Chapter 19 are, in the frequency view, learned spectral multipliers, and that equivalence is exactly what several of the architectures in this section's research frontier exploit.

5. Phase Correlation: Alignment in One Transform Advanced

Here is the frequency domain solving a real task end to end. If image $b$ is image $a$ shifted by $(\Delta x, \Delta y)$, the Fourier shift theorem says their spectra differ only by a phase ramp: $B(u,v) = A(u,v)\, e^{-j2\pi(u\Delta x/M + v\Delta y/N)}$. Dividing out the magnitudes leaves a pure complex exponential whose inverse transform is a single sharp spike located exactly at the displacement. This is phase correlation: translation estimation for the price of three FFTs, robust to global illumination changes (which live in magnitude, not phase) and indifferent to how feature-rich the scene is.

import cv2
import numpy as np
from skimage import data

a = data.camera().astype(np.float32)
b = np.roll(a, shift=(7, 15), axis=(0, 1))     # down 7 pixels, right 15

(shift_x, shift_y), response = cv2.phaseCorrelate(a, b)
print(f"shift = ({shift_x:.2f}, {shift_y:.2f}), response = {response:.3f}")
# Recovers the 15-pixel horizontal and 7-pixel vertical displacement to
# sub-pixel accuracy with a response near 1.0. The sign convention
# follows the argument order: verify it once on a synthetic roll like
# this one before trusting it in production.