Fourier Transform

The Fourier transform rewrites an image as a sum of sinusoidal gratings: plane waves of brightness, each with its own spatial frequency, orientation, amplitude, and phase. The spectrum holds the same information as the image, but some operations are far simpler on it: convolution becomes multiplication, a translation becomes a phase ramp, and blurring, sampling, and periodic noise show up in where the energy sits. The discrete version is computed with the fast Fourier transform (FFT).

Definition

The two-dimensional Fourier transform expresses an image as a weighted superposition of complex exponentials $e^{i 2\pi (ux + vy)}$, one for every pair of spatial frequencies $(u, v)$. The complex weight of each, the spectrum $F(u, v)$, has a magnitude (the grating’s amplitude) and an angle (its phase). The inverse transform recovers the image exactly. The discrete Fourier transform (DFT) does the same for a digital image with as many frequencies as there are pixels.

Intuition

A single grating, such as $\cos(2\pi(ux + vy))$, looks like a set of parallel stripes. The vector $(u, v)$ points perpendicular to the stripes; its length is the number of cycles per unit length, the reciprocal of the stripe spacing. Low frequencies are slowly varying shading; high frequencies are fine texture, sharp edges, and noise. The amplitude sets the stripes’ contrast, and the phase slides them sideways without changing their spacing or direction.

Any image can be built from such gratings. A smooth sky needs mostly low frequencies. A sharp vertical edge needs gratings of many frequencies, all with horizontal wave vectors, whose crests line up at the edge and cancel elsewhere. The spectrum of a photograph therefore has most of its energy near the origin, with bright lines perpendicular to the dominant edges. Because a real cosine is the sum of two complex exponentials, a single grating appears as two dots, at $(u, v)$ and $(-u, -v)$.

Formal Definition

For a function of two continuous variables, with frequency in cycles per unit length, a negative sign in the forward exponent, and no normalization constant,

$$F(u, v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y)\, e^{-i 2\pi (ux + vy)}\, dx\, dy, \qquad f(x, y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(u, v)\, e^{i 2\pi (ux + vy)}\, du\, dv .$$

For an $M \times N$ digital image $f[m, n]$, with row index $m$ and column index $n$, the DFT and its inverse are

$$F[k, l] = \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f[m, n]\, e^{-i 2\pi \left( \frac{km}{M} + \frac{ln}{N} \right)}, \qquad f[m, n] = \frac{1}{MN} \sum_{k=0}^{M-1} \sum_{l=0}^{N-1} F[k, l]\, e^{i 2\pi \left( \frac{km}{M} + \frac{ln}{N} \right)} .$$

The index $l$ counts cycles across the width of the image and $k$ cycles down its height. Indices above half the size stand for negative frequencies: $F[M - k, l]$ is the coefficient at vertical frequency $-k$. This is NumPy’s convention for fft2 and ifft2; other texts move the $1/(MN)$ to the forward transform, split it as $1/\sqrt{MN}$, or use angular frequency $\omega = 2\pi u$. The magnitude $|F[k, l]|$ is the amplitude spectrum, its square the power spectrum, and $\arg F[k, l]$ the phase spectrum.

The DFT treats the image as one period of an infinite periodic pattern, since its basis functions repeat every $M$ rows and $N$ columns.

Properties

These hold for the DFT with indices modulo the image size, and for the continuous transform without the periodicity. Bracewell (2003) treats the two-dimensional theorems in detail.

• Linearity. The transform of $a f + b g$ is $a F + b G$.
• Shift theorem. A circular shift by $(m_0, n_0)$ multiplies the spectrum by a phase ramp, $$f[m - m_0,\, n - n_0] \;\longleftrightarrow\; F[k, l]\, e^{-i 2\pi \left( \frac{k m_0}{M} + \frac{l n_0}{N} \right)} .$$ The magnitude is unchanged: the amplitude spectrum says what gratings an image contains, not where.
• Convolution theorem. Circular convolution in space is multiplication in frequency, and multiplication in space is convolution in frequency. See convolution for its use in fast filtering.
• Parseval’s theorem. Energy is preserved up to the normalization: $\sum_{m,n} |f[m, n]|^2 = \frac{1}{MN} \sum_{k,l} |F[k, l]|^2$. The power spectrum distributes an image’s energy over frequencies.
• Separability. The kernel $e^{-i2\pi(km/M + ln/N)}$ factors into a row term and a column term, so the 2D DFT is a 1D DFT of every row followed by one of every column.
• Scaling. Stretching an image compresses its spectrum: $f(ax, by)$ has transform $\frac{1}{|ab|} F\!\left(\frac{u}{a}, \frac{v}{b}\right)$. A narrow kernel has a broad frequency response, and no function is narrow in both domains: the Gaussian minimizes the product of the two spreads (the uncertainty principle).
• Rotation. Rotating an image by an angle rotates its continuous spectrum by the same angle. On a finite grid this holds only approximately.
• Symmetry for real images. For real $f$, $F[-k, -l] = \overline{F[k, l]}$: the magnitude is symmetric about the origin and the phase antisymmetric, so half the spectrum suffices, which real-input transforms such as rfft2 exploit.

Magnitude and phase

The amplitude spectrum is the more intuitive half, but most of the structure lives in the phase. Oppenheim and Lim (1981) reviewed evidence that many important features of a signal survive when only its Fourier phase is kept, illustrated with phase-only reconstructions of images, speech, and crystallographic structures, and results showing that under certain conditions, such as finite extent, the phase alone determines a signal to within a scale factor. The shift theorem suggests why: phase encodes where the crests of every grating lie, and edges exist only where many gratings line up. Amplitude spectra, by contrast, are much alike: in a large set of natural images, van der Schaaf and van Hateren (1996) found power falling off roughly as $1/f^2$ (amplitude as $1/f$), with variation between images. In the example below, a hybrid of one image’s magnitude and another’s phase resembles the image that supplied the phase.

The Fast Fourier Transform

Evaluating the DFT sum directly costs $O(N^2)$ operations for $N$ samples, and $O(N^4)$ for an $N \times N$ image. Cooley and Tukey (1965) gave an algorithm for lengths that factor into small numbers, which splits the transform into shorter ones and needs a number of operations proportional to $N \log N$. By separability, a 2D FFT applies $2N$ such 1D transforms and costs $O(N^2 \log N)$. Libraries are fastest for sizes with only small prime factors.

Practical Issues

Periodic extension and windowing

Since the DFT assumes the image tiles the plane, its left border is adjacent to the right and the top to the bottom. Unless they match, the tiling has seams, and a seam needs energy at all frequencies perpendicular to it. This is why the spectrum of almost any photograph shows a bright cross along the horizontal and vertical frequency axes that has nothing to do with the content. The same mechanism causes spectral leakage: a grating that does not complete a whole number of cycles across the image does not land on a single DFT frequency, and its energy spreads over many.

The usual remedy is to multiply the image by a window, such as a 2D Hann window, that falls smoothly to zero at the borders, at the cost of slightly broader peaks and a weighting that favors the center. Mirror padding, which makes the borders continuous, is an alternative.

Display

The DFT stores zero frequency at index $[0, 0]$; for display, the quadrants are swapped to put it at the center (fftshift in NumPy). The amplitude spectrum spans many orders of magnitude, so it is shown as $\log(1 + |F|)$.

Frequency-Domain Filtering

A linear filter can be specified as a transfer function $H[k, l]$ that multiplies the spectrum, which is equivalent to convolving with its inverse transform. Nonlinear filters, covered in image filtering, have no such description.

• Ideal low-pass. Setting $H = 1$ inside a radius $D_0$ and $0$ outside passes low frequencies exactly and removes all others. Its spatial kernel is a sinc-like function with slowly decaying oscillating lobes, so the filtered image rings: edges acquire ripples parallel to them.
• Gaussian low-pass. The transform of a Gaussian is a Gaussian: $e^{-(x^2 + y^2)/(2\sigma^2)}$ transforms to $2\pi\sigma^2 e^{-2\pi^2 \sigma^2 (u^2 + v^2)}$. With no sharp cutoff and no negative lobes, the Gaussian filter smooths without ringing.
• High-pass. $1 - H$ for a low-pass $H$ keeps edges and texture; adding the result back to the image sharpens it.
• Notch filters. Periodic interference, such as scan lines or halftone screens, concentrates in a few isolated spectral spikes. Zeroing small neighborhoods around them and their symmetric partners removes the pattern while barely affecting the image; broadband image noise cannot be isolated this way.

Examples

The code checks grating spectra, leakage, the shift theorem, and a magnitude–phase swap on synthetic images, then recovers a known translation by phase correlation (described below).

import numpy as np

rng = np.random.default_rng(0)
N = 128
y, x = np.mgrid[0:N, 0:N]

# 1. A grating's spectrum is a pair of peaks at +/- its frequency.
def peaks(img):
    F = np.fft.fftshift(np.fft.fft2(img))
    F[N // 2, N // 2] = 0                       # ignore the mean (DC)
    mag = np.abs(F)
    idx = np.argsort(mag.ravel())[-2:]          # the two strongest bins
    return sorted((int(c) - N // 2, int(r) - N // 2) for r, c in
                  zip(*np.unravel_index(idx, F.shape)))

for u, v in [(8, 0), (0, 8), (12, 5)]:          # cycles per image along x and y
    g = 0.5 + 0.5 * np.cos(2 * np.pi * (u * x + v * y) / N)
    angle = np.degrees(np.arctan2(v, u))
    print(f"grating (u, v) = ({u:2d}, {v}): period {N / np.hypot(u, v):5.2f} px, "
          f"wave vector at {angle:4.1f} deg -> peaks at {peaks(g)}")

# A non-integer frequency is not periodic on the grid: its energy leaks across the spectrum.
k = np.fft.fftfreq(N, 1 / N)                    # integer frequencies, cycles per image
def far_leakage(img, u):
    """Share of energy more than 3 bins away from the true peaks at (+-u, 0)."""
    P = np.abs(np.fft.fft2(img)) ** 2
    near_x = np.minimum(abs(k - u), abs(k + u)) <= 3
    near = (abs(k)[:, None] <= 3) & near_x[None, :]
    return P[~near].sum() / P.sum()

hann = np.outer(np.hanning(N), np.hanning(N))
for u in [8, 8.5]:
    g = np.cos(2 * np.pi * u * x / N)
    print(f"u = {u}: energy leaked > 3 bins from the peaks = {far_leakage(g, u):.4f}, "
          f"with Hann window = {far_leakage(g * hann, u):.4f}")

# 2. Shift theorem: a circular shift changes only the phase.
img = rng.random((N, N))
dy, dx = 5, 12
shifted = np.roll(img, (dy, dx), axis=(0, 1))
F, G = np.fft.fft2(img), np.fft.fft2(shifted)
fy = np.fft.fftfreq(N)[:, None]                 # cycles per pixel
fx = np.fft.fftfreq(N)[None, :]
predicted = F * np.exp(-2j * np.pi * (fx * dx + fy * dy))
print(f"shift ({dy}, {dx}): max |magnitude change| = {np.abs(np.abs(G) - np.abs(F)).max():.1e}, "
      f"max phase change = {np.abs(np.angle(G / F)).max():.2f} rad, "
      f"max error vs predicted spectrum = {np.abs(G - predicted).max():.1e}")

# 3. Swap magnitude and phase between two images.
def shapes(seed, n_shapes=12):
    r = np.random.default_rng(seed)
    out = np.zeros((N, N))
    for _ in range(n_shapes):
        cx, cy, rad = r.uniform(10, N - 10, 2).tolist() + [r.uniform(4, 16)]
        if r.random() < 0.5:
            out[(x - cx) ** 2 + (y - cy) ** 2 < rad ** 2] += r.uniform(0.3, 1)
        else:
            out[(abs(x - cx) < rad) & (abs(y - cy) < rad / 2)] += r.uniform(0.3, 1)
    return out

A, B = shapes(1), shapes(2)
corr = lambda p, q: np.corrcoef(p.ravel(), q.ravel())[0, 1]
for (m, mname), (p, pname) in [((A, "A"), (B, "B")), ((B, "B"), (A, "A"))]:
    hybrid = np.real(np.fft.ifft2(np.abs(np.fft.fft2(m)) * np.exp(1j * np.angle(np.fft.fft2(p)))))
    print(f"magnitude of {mname}, phase of {pname}: corr with A = {corr(hybrid, A):+.2f}, "
          f"corr with B = {corr(hybrid, B):+.2f}")

# 4. Phase correlation recovers a translation.
scene = shapes(3, 30) + 0.05 * rng.standard_normal((N, N))
moved = np.roll(scene, (-9, 23), axis=(0, 1)) + 0.05 * rng.standard_normal((N, N))
R = np.fft.fft2(moved) * np.conj(np.fft.fft2(scene))
R /= np.abs(R) + 1e-12                          # keep only the phase difference
surface = np.real(np.fft.ifft2(R))
py, px = np.unravel_index(surface.argmax(), surface.shape)
py, px = (py + N // 2) % N - N // 2, (px + N // 2) % N - N // 2   # wrap to signed shifts
second = np.sort(surface.ravel())[-2]
print(f"phase correlation: true shift (-9, 23), estimated ({py}, {px}), "
      f"peak {surface.max():.2f} vs next highest {second:.2f}")

Output:

grating (u, v) = ( 8, 0): period 16.00 px, wave vector at  0.0 deg -> peaks at [(-8, 0), (8, 0)]
grating (u, v) = ( 0, 8): period 16.00 px, wave vector at 90.0 deg -> peaks at [(0, -8), (0, 8)]
grating (u, v) = (12, 5): period  9.85 px, wave vector at 22.6 deg -> peaks at [(-12, -5), (12, 5)]
u = 8: energy leaked > 3 bins from the peaks = 0.0000, with Hann window = 0.0000
u = 8.5: energy leaked > 3 bins from the peaks = 0.0683, with Hann window = 0.0001
shift (5, 12): max |magnitude change| = 4.3e-14, max phase change = 3.14 rad, max error vs predicted spectrum = 4.3e-13
magnitude of A, phase of B: corr with A = +0.13, corr with B = +0.75
magnitude of B, phase of A: corr with A = +0.75, corr with B = +0.09
phase correlation: true shift (-9, 23), estimated (-9, 23), peak 0.83 vs next highest 0.04

Each grating with a whole number of cycles produces two peaks, at $(u, v)$ and $(-u, -v)$. A grating with 8.5 cycles per image leaks about 7% of its energy more than three bins away from its peaks; the Hann window reduces that to 0.01%. A circular shift leaves the magnitude unchanged to rounding error while changing phases by up to $\pi$, matching the predicted phase ramp. Each hybrid correlates at 0.75 with its phase source and only 0.13 or 0.09 with its magnitude source. Phase correlation finds the shift exactly despite independent noise, with a peak about twenty times higher than any other value.

Common Misconceptions

• “The bright cross in a spectrum shows horizontal and vertical structure in the scene.” It usually comes from the seams of the periodic extension; windowing removes most of it.
• “FFT filtering is the same as spatial filtering.” Without padding, multiplying DFTs performs circular convolution, which wraps content around the borders.
• “Translating an image changes its amplitude spectrum.” A circular shift changes only the phase; an ordinary shift alters the magnitude only through content crossing the borders.

Where It Is Used

• Fast convolution. Large kernels are applied by multiplying spectra; see convolution.
• Image registration. Phase correlation, introduced by Kuglin and Hines in 1975, estimates the translation between two images from the normalized cross-power spectrum $$R[k, l] = \frac{G[k, l]\, \overline{F[k, l]}}{\left| G[k, l]\, \overline{F[k, l]} \right|} .$$ If $g$ is $f$ shifted by $(m_0, n_0)$, the shift theorem makes $R$ a pure phase ramp, whose inverse DFT is an impulse at $(m_0, n_0)$. Discarding the magnitude gives a sharp peak, insensitive to a global change of contrast. Reddy and Chatterji (1996) extended the method to rotation and scale by applying it to the Fourier magnitudes resampled on a log-polar grid, where rotation and scaling become translations. It is used for motion estimation, stabilization, and mosaicking.
• Sampling and aliasing. Sampling replicates the spectrum, and overlapping copies cause aliasing, as sampling and aliasing explains; the smoothing step of an image pyramid is the low-pass filter it calls for.
• Compression. JPEG and related codecs use the discrete cosine transform, a close relative of the DFT that implicitly extends a block by mirror reflection rather than tiling. This avoids border seams and concentrates the energy of smooth blocks in a few low-frequency coefficients.
• Texture analysis. Energy in rings and wedges of the power spectrum summarizes a texture’s dominant scales and orientations.
• Deconvolution. Blur by a known point spread function multiplies the spectrum by its transfer function; dividing by it naively amplifies noise where that function is small, so methods such as Wiener filtering regularize the division.
• Neural networks. Tancik et al. (2020) showed that a multilayer perceptron taking low-dimensional coordinates as input learns high-frequency components very slowly, and that first mapping the coordinates through sinusoids of several frequencies (Fourier features) lets it learn fine detail. Coordinate networks for images and 3D scenes use such mappings.