Lesson 15: The Fourier Transform and Frequency-Domain Filtering¶
Every filter so far has worked directly on pixel neighborhoods. The Fourier transform offers a completely different view: any image can be rewritten as a sum of sinusoidal gratings at different frequencies and orientations. Filtering can then be done by directly boosting or suppressing specific frequencies — and, thanks to the convolution theorem, this turns out to be mathematically the same operation as spatial convolution, just seen from a different angle.
import numpy as np
import cv2
import matplotlib.pyplot as plt
1D warm-up: a signal as a sum of sinusoids¶
Any signal can be approximated by adding up sine and cosine waves of different frequencies. For example, a square wave is the sum of sine waves at increasing frequencies (odd harmonics, decreasing amplitude). The Fourier transform is the tool that goes the other way: given a signal, it tells you exactly which frequencies (and how much of each) are present.
t = np.linspace(0, 1, 500, endpoint=False)
square_wave = np.sign(np.sin(2 * np.pi * 5 * t))
approx = np.zeros_like(t)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
for n_harmonics in [1, 9]:
approx = sum((4 / (np.pi * k)) * np.sin(2 * np.pi * 5 * k * t) for k in range(1, n_harmonics + 1, 2))
axes[0 if n_harmonics == 1 else 1].plot(t, square_wave, '--', color='gray', label='true square wave')
axes[0 if n_harmonics == 1 else 1].plot(t, approx, label=f'{n_harmonics} harmonic(s)')
axes[0 if n_harmonics == 1 else 1].legend(fontsize=8)
axes[0 if n_harmonics == 1 else 1].set_title(f'{n_harmonics} harmonic(s)', fontsize=9)
plt.tight_layout()
plt.show()
The 2D Fourier transform of an image¶
The 2D discrete Fourier transform (DFT) of an image decomposes it into 2D sinusoidal gratings at every combination of horizontal and vertical frequency. The algorithm for computing the DFT is the fast Fourier transform (FFT) (np.fft.fft2), which results in a complex-valued (real + imaginary) 2D array. Similar to the image gradient, we usually view the Fourier transform's magnitude (how much of each frequency is present) and orientation. Two common tricks for visualization: show the magnitude on a log scale, and shift the zero-frequency (average brightness) term from the top-left corner to the center.
def show_spectrum(image, ax, title):
F = np.fft.fftshift(np.fft.fft2(image)) # compute discrete Fourier transform (DFT), and shift DC value to center
magnitude = np.log1p(np.abs(F)) # display logarithm of magnitude
ax.imshow(magnitude, cmap='gray')
ax.set_title(title, fontsize=9)
ax.axis('off')
return F
img = np.zeros((200, 200), dtype=np.float64)
cv2.rectangle(img, (60, 60), (140, 140), 200, -1)
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(img, cmap='gray')
axes[0].set_title('Image')
axes[0].axis('off')
show_spectrum(img, axes[1], 'Log-magnitude spectrum')
plt.tight_layout()
plt.show()
The bright cross through the center comes from the square's sharp horizontal and vertical edges — a hard step edge contains energy at every frequency along the direction perpendicular to it, which is why a sharp-edged shape has such a spread-out spectrum.
A grating reveals its frequency directly as a spectrum peak¶
A pure sinusoidal grating is the simplest possible test: its spectrum should be (almost) a single pair of bright dots, at a distance from the center equal to its frequency, in the direction perpendicular to its stripes.
Why a pair, not one dot? A real-valued sine wave is, by Euler's formula, a mix of two complex exponentials in equal amounts: $\sin(2\pi f x) = \frac{1}{2i}\left(e^{i2\pi fx} - e^{-i2\pi fx}\right)$. The Fourier transform of the first exponential is a spike at frequency $f$, and for the other it is at $-f$. Therefore, a real sinusoid produces two spikes at $+f$ and $-f$, symmetric about the zero-frequency (DC) center. In fact, for any real-valued image, its Fourier magnitude is symmetric: $|F(-f)|=|F(f)|$.
size = 100
x = np.arange(size)
cycles_across_image = 12
grating = np.tile(np.sin(2 * np.pi * cycles_across_image * x / size), (size, 1))
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(grating, cmap='gray')
axes[0].set_title(f'Grating, {cycles_across_image} cycles across')
axes[0].axis('off')
F_grating = show_spectrum(grating, axes[1], 'Spectrum')
plt.tight_layout()
plt.show()
peak = np.unravel_index(np.argmax(np.abs(F_grating)), F_grating.shape)
print(f'brightest spectrum point at offset {peak[1] - size // 2} from center '
f'(grating frequency was {cycles_across_image})')
Filtering by editing the spectrum¶
An ideal low-pass filter keeps only frequencies inside a circle around the center (blurs); an ideal high-pass filter keeps only frequencies outside it (edge-enhances).
def circular_mask(shape, radius, keep_inside=True):
h, w = shape
yy, xx = np.mgrid[:h, :w]
dist = np.sqrt((yy - h / 2)**2 + (xx - w / 2)**2)
inside = dist <= radius
return inside if keep_inside else ~inside
photo = np.zeros((200, 200), dtype=np.float64)
cv2.rectangle(photo, (40, 40), (120, 120), 200, -1)
cv2.circle(photo, (150, 150), 30, 120, -1)
low_pass_mask = circular_mask(photo.shape, radius=15, keep_inside=True)
high_pass_mask = circular_mask(photo.shape, radius=15, keep_inside=False)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
for ax, im, title in zip(axes, [low_pass_mask, high_pass_mask],
['Ideal low-pass mask\n(frequency domain)', 'Ideal high-pass mask\n(frequency domain)']):
ax.imshow(im, cmap='gray')
ax.set_title(title, fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
Once an image is transformed to the frequency domain, we can filter it by masking frequencies directly, then transforming back with the inverse transform. Since the mask is the same size as the image, we can simply elementwise-multiply in the frequency domain.
def apply_frequency_filter(image, mask):
F = np.fft.fftshift(np.fft.fft2(image)) # forward transform
filtered = np.fft.ifft2(np.fft.ifftshift(F * mask)) # filter, then inverse transform
return np.real(filtered)
low_pass = apply_frequency_filter(photo, low_pass_mask)
high_pass = apply_frequency_filter(photo, high_pass_mask)
fig, axes = plt.subplots(1, 3, figsize=(9, 3.5))
for ax, im, title in zip(axes, [photo, low_pass, high_pass],
['Original', 'Ideal low-pass\n(blurred)', 'Ideal high-pass\n(edges only)']):
ax.imshow(im, cmap='gray')
ax.set_title(title, fontsize=9)
ax.axis('off')
plt.tight_layout()
plt.show()
The cost of a sharp cutoff: ringing¶
Look closely at these results: faint ripples radiate from the shape edges. An ideal filter has a perfectly sharp cutoff in frequency, which corresponds (by the same convolution theorem) to convolving with a spatial kernel that has infinite ripples of its own (a sinc function) — this is the Gibbs phenomenon. A Gaussian mask, which falls off smoothly instead of cutting off sharply, avoids this.
def gaussian_mask(shape, sigma):
h, w = shape
yy, xx = np.mgrid[:h, :w]
dist2 = (yy - h / 2)**2 + (xx - w / 2)**2
return np.exp(-dist2 / (2 * sigma**2))
smooth_low_pass = apply_frequency_filter(photo, gaussian_mask(photo.shape, sigma=15))
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
axes[0].imshow(low_pass, cmap='gray')
axes[0].set_title('Ideal (sharp cutoff)\nvisible ringing near edges')
axes[1].imshow(smooth_low_pass, cmap='gray')
axes[1].set_title('Gaussian mask\nsmooth falloff, no ringing')
for ax in axes:
ax.axis('off')
plt.tight_layout()
plt.show()
The convolution theorem¶
Why does frequency-domain filtering work? The convolution theorem
$$\mathcal{F}\{f * g\} = \mathcal{F}\{f\} \cdot \mathcal{F}\{g\}$$
says that convolving two signals in the spatial domain is exactly the same as multiplying their Fourier transforms elementwise, then transforming back — convolution becomes multiplication, and vice versa.
A catch: the DFT gives circular convolution¶
The DFT implicitly treats the image as if it were infinitely repeating (imagine seamless tiles extending forever in both directions). There's no "outside the image", only "the next copy of the same image." So multiplying two DFTs and transforming back doesn't give the ordinary ("linear") convolution from Lesson 10, where the kernel eventually slides off the edge into some explicit border (zero, replicate, reflect); rather, it gives circular convolution, where the kernel that slides off the right edge wraps around and picks up pixels from the left edge instead. The two only agree if the border handling on the spatial side is also wraparound (cv2.BORDER_WRAP).
rng = np.random.default_rng(0)
test_img = rng.integers(0, 255, (64, 64)).astype(np.float64)
sigma, ksize = 3, 15
ax = np.arange(ksize) - ksize // 2
gx, gy = np.meshgrid(ax, ax)
kernel = np.exp(-(gx**2 + gy**2) / (2 * sigma**2))
kernel /= kernel.sum()
# Spatial convolution, wrap-around border to match the FFT's implicit circular convolution
spatial_result = cv2.filter2D(test_img, cv2.CV_64F, kernel, borderType=cv2.BORDER_WRAP)
# Frequency-domain: pad the kernel to image size, center it at (0,0) via np.roll, then multiply spectra
kernel_padded = np.zeros_like(test_img)
kernel_padded[:ksize, :ksize] = kernel
kernel_padded = np.roll(kernel_padded, -(ksize // 2), axis=0)
kernel_padded = np.roll(kernel_padded, -(ksize // 2), axis=1)
freq_result = np.real(np.fft.ifft2(np.fft.fft2(test_img) * np.fft.fft2(kernel_padded)))
print(f'max abs difference between spatial and frequency-domain convolution: '
f'{np.abs(spatial_result - freq_result).max():.2e}')
In practice, though, frequency-domain filtering is rarely used in everyday computer vision. Most kernels — box blur, Gaussian, Sobel, sharpen — are small (3x3 or 5x5), and direct spatial convolution is already fast at those sizes, so the fixed overhead of two forward FFTs and one inverse FFT usually isn't worth paying. The DFT earns its keep mainly once kernels get large.
Exercise¶
- Rotate the grating in the spectrum-peak experiment (make the stripes diagonal instead of vertical) by building it as a function of
x*cos(theta) + y*sin(theta). Predict, then verify, where the spectrum peak moves to. - Build a band-pass frequency mask (an annulus: keep frequencies between an inner and outer radius, block everything else) and apply it to
photo. What kind of image content survives? - Try
apply_frequency_filterwith a very small Gaussian mask sigma (e.g. 3) versus a large one (e.g. 60) on a photo-like test image, and describe the trend as sigma increases toward the image size. - Get linear (not circular) convolution out of an FFT: zero-pad both
test_imgandkernelto at least(image_size + kernel_size - 1)along each axis before transforming, multiply the spectra, inverse-transform, then crop back down. Compare the result tocv2.filter2D(..., borderType=cv2.BORDER_CONSTANT)(zero-padding border) instead ofBORDER_WRAP. - Using
%timeit, measure direct convolution (cv2.filter2D) against FFT-based convolution (as infreq_resultabove) for Gaussian kernel sizes 5, 15, 31, 63, and 127 on a 256x256 image. At roughly what kernel size does the FFT approach start winning?