Fourier Filtering for Noise Removal - Image Processing in the Frequency Domain
Fourier Transform and Frequency Representation of Images
The Fourier Transform converts images from spatial domain to frequency domain. Images represented as pixel intensity values in spatial domain become superpositions of sinusoidal waves at various frequencies in frequency domain. This transformation enables direct manipulation of frequency components for noise removal and edge enhancement.
2D Discrete Fourier Transform (DFT):
The 2D DFT of an MxN image f(x,y) is F(u,v) = sum sum f(x,y) x exp(-j2pi(ux/M + vy/N)). F(u,v) is complex-valued at frequency (u,v) with amplitude |F(u,v)| and phase arg(F(u,v)). Amplitude spectrum represents each frequency component's strength while phase spectrum represents positional information.
Frequency meaning:
Low-frequency components correspond to gradual intensity changes (large structures, overall brightness). High-frequency components correspond to rapid changes (edges, textures, noise). DC component F(0,0) represents average image intensity. In spectrum images, center shows low frequencies while periphery shows high frequencies.
FFT fast computation:
Direct DFT computation is O(N^4) but FFT (Fast Fourier Transform) reduces this to O(N^2 log N). NumPy's np.fft.fft2() and OpenCV's cv2.dft() implement FFT. A 1024x1024 image FFT completes in milliseconds on modern CPUs.
Low-Pass Filters - Noise Removal Fundamentals
Low-Pass Filters pass low-frequency components while attenuating high frequencies. Since noise generally appears as high-frequency components, low-pass filtering removes noise. However, edges also contain high frequencies, creating a noise-removal vs. edge-preservation tradeoff.
Ideal low-pass filter:
Passes all components below cutoff frequency D0 and blocks all above. H(u,v) = 1 (D(u,v) <= D0), H(u,v) = 0 (D(u,v) > D0). However, abrupt cutoff causes ringing (Gibbs phenomenon), producing ring-shaped artifacts. Rarely used in practice.
Butterworth low-pass filter:
H(u,v) = 1 / (1 + (D(u,v)/D0)^(2n)) provides smooth cutoff characteristics. Higher order n means sharper cutoff - n=1 is smoothest, n approaching infinity approaches ideal filter. n=2 is widely used as practical balance. Minimal ringing with natural blur effect.
Gaussian low-pass filter:
H(u,v) = exp(-D(u,v)^2 / (2*D0^2)) provides the smoothest cutoff with zero ringing. Equivalent to spatial-domain Gaussian blur (cv2.GaussianBlur) with relationship sigma = M/(2*pi*D0). Noise removal quality comparable to Butterworth with ringing-free advantage.
Cutoff frequency selection:
Cutoff frequency D0 significantly affects results. Smaller D0 increases noise removal but also increases blur. Common approaches include selecting D0 containing 90-95% of power spectrum, or determining D0 from noise spectral characteristics.
High-Pass and Band-Pass Filters
High-pass filters pass high frequencies while attenuating low frequencies, used for edge detection and sharpening. Band-pass filters pass specific frequency bands for extracting particular textures or structures from images.
High-pass filter:
Defined as H_HP(u,v) = 1 - H_LP(u,v) where H_LP is a low-pass filter. Emphasizes edges and textures while removing overall brightness (DC component). Results appear dark, so typically used as unsharp mask added to original: enhanced = original + k x highpass (k is sharpening strength).
Band-pass filter:
Defined as H_BP(u,v) = H_LP(D_H) - H_LP(D_L), passing only frequencies between D_L and D_H. Used for extracting specific-scale textures or isolating frequency-band noise. Applied in medical imaging structure enhancement and remote sensing land classification.
Band-reject filter:
Defined as H_BR(u,v) = 1 - H_BP(u,v), removing specific frequency bands. Particularly effective for periodic noise removal (power hum stripe patterns). When noise frequency is known, removing only that band preserves other image information while eliminating noise.
Homomorphic filtering:
Separates images into illumination (low-frequency) and reflectance (high-frequency) components for independent processing. Applies Fourier transform after logarithmic transformation, filtering to suppress low frequencies and enhance high frequencies. Effective for contrast improvement in non-uniform illumination, used in face recognition preprocessing.
Notch Filters - Periodic Noise Removal
Notch Filters remove specific points (notches) in the frequency spectrum. Specialized for periodic noise removal (moire, electromagnetic interference stripes, scanning banding artifacts).
Periodic noise characteristics:
Periodic noise appears as distinct peaks (spikes) in frequency spectrum. Horizontal stripes produce peaks on the vertical axis at specific frequencies; diagonal stripes produce peaks at corresponding angular positions. Identifying and removing these peaks effectively eliminates periodic noise.
Notch filter design:
Notch filters remove small regions centered at specific spectrum points (u0, v0) and their symmetric counterparts (-u0, -v0). Due to Fourier transform symmetry, real image spectra are origin-symmetric, so notches are always set in pairs. Gaussian notch H(u,v) = 1 - exp(-D1^2/(2*sigma^2)) x exp(-D2^2/(2*sigma^2)) achieves smooth removal.
Notch position identification:
Periodic noise notch positions are identified by visualizing power spectrum (|F(u,v)|^2) manually or via automatic detection algorithms. Automatic detection searches for local maxima after excluding DC component, using peak detection or threshold-based methods.
Practical example - scan moire removal:
Scanned printed materials exhibit moire from interference between print dot patterns and scanner resolution. Spectrum observation reveals regular peaks corresponding to dot frequencies. Applying notch filters to these peaks removes moire while preserving image content. Notch width (sigma) is adjusted to match peak spread.
Wiener Filter and Inverse Filtering
Wiener Filter is an optimal filtering method considering noise and signal statistical properties. Unlike simple low-pass filters, it determines optimal attenuation per frequency based on signal-to-noise ratio (SNR).
Inverse filtering problem:
The simplest restoration from degraded image G(u,v) = H(u,v)F(u,v) + N(u,v) is inverse filter F_hat = G/H. However, where H(u,v) approaches zero, noise N/H is amplified catastrophically. This is inverse filtering's fundamental problem.
Wiener filter formulation:
Wiener filter minimizes mean squared error between restored and original images. W(u,v) = H*(u,v) / (|H(u,v)|^2 + S_n(u,v)/S_f(u,v)) where H* is complex conjugate, S_n is noise power spectrum, S_f is signal power spectrum. Approaches inverse filter at high-SNR frequencies while strongly attenuating at low-SNR frequencies.
Regularization parameter:
In practice, noise and signal power spectra are often unknown, so S_n/S_f is approximated by constant K. W(u,v) = H*(u,v) / (|H(u,v)|^2 + K) where K adjusts for noise level. Larger K means stronger noise suppression (more blur); smaller K means sharper restoration (more noise remains).
Application - motion blur removal:
Camera shake motion blur appears as directional sinc functions in frequency domain. Estimating blur direction and length, constructing corresponding degradation function H, and applying Wiener filter restores blurred images. Blur parameter estimation uses spectrum zero-crossing pattern analysis or cepstrum analysis.
Signal processing and Fourier analysis textbooks can be found on Amazon
Implementation Guide - Fourier Filtering in Python
Implementing Fourier filtering with NumPy and OpenCV. Covers FFT computation through filter design to inverse transformation for practical processing pipelines.
FFT and inverse FFT basics:
NumPy FFT computes via F = np.fft.fft2(image), centering low frequencies with F_shift = np.fft.fftshift(F). After filtering: F_ishift = np.fft.ifftshift(F_filtered) then result = np.abs(np.fft.ifft2(F_ishift)) for spatial domain inverse transformation.
Gaussian low-pass filter implementation:
Generate filter mask by computing distance map from image center and applying Gaussian function. rows, cols = image.shape; crow, ccol = rows//2, cols//2; D = distance_map; H = np.exp(-D**2/(2*D0**2)) then apply via F_filtered = F_shift * H.
Notch filter implementation:
Visualize spectrum to identify periodic noise peak positions, placing Gaussian notches at those locations. Visualize with magnitude = 20*np.log(np.abs(F_shift)+1), identify peak coordinates (u0, v0). Generate notch filter as H = 1 - np.exp(-D1**2/(2*sigma**2)) * np.exp(-D2**2/(2*sigma**2)).
Performance considerations:
Large images may have FFT computation time issues. Padding to power-of-2 sizes optimizes FFT efficiency. optimal_size = cv2.getOptimalDFTSize(n) finds optimal sizes. Real image FFT uses np.fft.rfft2() exploiting symmetry to halve computation. GPU acceleration via CuPy's cupyx.scipy.fft.fft2() provides 10-50x speedup over CPU.