Frequency Domain Restoration — Inverse and Wiener Filtering
Why Work in the Frequency Domain
In the previous post we built the image degradation model: the observed image $g(x,y)$ is the original image $f(x,y)$ convolved with a degradation function $h(x,y)$ and corrupted by additive noise $n(x,y)$:
$$g(x,y) = h(x,y) * f(x,y) + n(x,y)$$
Restoration means: given $g$ and $h$, recover $f$ as closely as possible.
In the spatial domain, this requires solving a deconvolution problem. Deconvolution is a massive linear system. For a $512 \times 512$ image, you face 260,000 unknowns. Direct solving is practically impossible.
But convolution has a beautiful property: the convolution theorem. Spatial convolution equals frequency multiplication:
$$g(x,y) = h(x,y) * f(x,y) \quad \Longleftrightarrow \quad G(u,v) = H(u,v) \cdot F(u,v)$$
where $G, H, F$ are the Fourier transforms of $g, h, f$.
This tells us something crucial: in the frequency domain, convolution degradation becomes pointwise multiplication. Pointwise multiplication can be reversed by pointwise division, and division is far simpler than solving a linear system. This is the fundamental motivation for frequency domain restoration.
The noise term carries over to the frequency domain as well:
$$G(u,v) = H(u,v) \cdot F(u,v) + N(u,v)$$
Every difficulty in frequency domain restoration comes from this $N(u,v)$.
Inverse Filtering: The Direct Approach
The Idea
Since degradation is multiplication, restoration is division. The most straightforward approach divides out the degradation function $H(u,v)$:
$$\hat{F}(u,v) = \frac{G(u,v)}{H(u,v)}$$
With zero noise, $G = H \cdot F$, so $\hat{F} = F$, perfect recovery. This formula is beautifully simple.
The Fatal Flaw
Reality has noise. Substituting $G = H \cdot F + N$ into the inverse filter:
$$\hat{F}(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)}$$
The first term is the original image, no problem. The second term is noise divided by $H$. That is where things break.
Degradation functions $H(u,v)$ are typically very small at high frequencies, often approaching zero. The reason is that most blur processes (defocus, motion blur) act as low-pass filters that suppress high-frequency components. When $H(u,v) \to 0$, the ratio $N(u,v) / H(u,v) \to \infty$, and noise explodes.
Image energy concentrates at low frequencies, while noise energy spreads across all frequencies. Inverse filtering amplifies weak high-frequency noise until it drowns the entire image, making the result unusable.
flowchart TD
A["Input image<br/>G = H·F + N"] --> B["Divide by H(u,v)<br/>Pointwise in frequency"]
B --> C{"H ≈ 0?<br/>High-freq region"}
C -->|No Low-freq| D["Low-freq components<br/>Recovered well"]
C -->|Yes High-freq| E["N/H → ∞<br/>Noise explodes"]
D --> F["Result: clear low-freq<br/>High-freq pure noise"]
E --> F
classDef input fill:#2196F3,color:#fff
classDef proc fill:#9C27B0,color:#fff
classDef check fill:#FF9800,color:#fff
classDef good fill:#4CAF50,color:#fff
classDef bad fill:#f44336,color:#fff
class A input
class B proc
class C check
class D good
class E bad
class F badThe failure of inverse filtering reveals a deep truth: restoration is fundamentally ill-posed. Information lost during degradation cannot be recovered unconditionally. Noise permanently corrupts high-frequency data. Any practical restoration method must trade off between signal recovery and noise suppression.
Wiener Filtering: Balancing Signal and Noise
Goal: Minimum Mean Square Error
Wiener filtering does not follow the extreme path of inverse filtering. Its goal is more measured: find an estimate $\hat{F}$ that minimizes the Mean Square Error (MSE) between the estimate and the true image:
$$\min ; E\left[ |F(u,v) - \hat{F}(u,v)|^2 \right]$$
This is a statistical optimization. It does not chase perfect recovery of a single image. Instead, it minimizes the expected error across the ensemble.
The Formula
After derivation (omitting the statistical signal processing details), the frequency domain Wiener filter is:
$$\hat{F}(u,v) = \left[ \frac{H^*(u,v)}{|H(u,v)|^2 + K} \right] G(u,v)$$
where:
- $H^*(u,v)$ is the complex conjugate of the degradation function
- $|H(u,v)|^2$ is the power spectrum of the degradation function
- $K$ is the ratio of noise power spectrum to signal power spectrum: $K = S_n(u,v) / S_f(u,v)$
You can think of $H^* / (|H|^2 + K)$ as a corrected inverse filter. Compared to $1/H$ from inverse filtering, the denominator has an extra $K$ term. That term is what makes everything work.
Intuition for the K Parameter
Understanding Wiener filtering means understanding how $K$ behaves. $K$ is the inverse of the signal-to-noise ratio. It controls the filter’s behavior at different frequencies:
When signal dominates noise ($|H|^2 \gg K$, high SNR at low frequencies):
$$\frac{H^}{|H|^2 + K} \approx \frac{H^}{|H|^2} = \frac{1}{H}$$
The Wiener filter degenerates into inverse filtering, recovering signal aggressively. Since signal is strong and noise is weak here, inverse filtering’s noise amplification is not a problem.
When noise dominates signal ($|H|^2 \ll K$, low SNR at high frequencies):
$$\frac{H^}{|H|^2 + K} \approx \frac{H^}{K} \to 0$$
The filter output approaches zero, actively abandoning these frequencies. Noise has overwhelmed the signal here, so any recovery attempt would only amplify noise.
flowchart TD
A["Input G(u,v)"] --> B["Wiener filter<br/>H* / (|H|² + K)"]
B --> C{"Frequency region"}
C -->|Low-freq |H|²≫K| D["≈ 1/H<br/>Inverse filter<br/>Recover signal"]
C -->|High-freq |H|²≪K| E["≈ H*/K → 0<br/>Suppress noise<br/>Abandon recovery"]
D --> F["Clear low-freq<br/>Clean high-freq"]
E --> F
classDef input fill:#2196F3,color:#fff
classDef proc fill:#9C27B0,color:#fff
classDef check fill:#FF9800,color:#fff
classDef good fill:#4CAF50,color:#fff
classDef suppress fill:#f44336,color:#fff
class A input
class B proc
class C check
class D good
class E suppress
class F goodThis is the elegance of Wiener filtering: it recovers at low frequencies and suppresses at high frequencies. The switch is decided automatically by the data itself (the SNR), with no need for a manually set cutoff frequency.
Practical K Values
In theory, $K = S_n / S_f$ requires knowledge of both noise and signal power spectra. In practice these are usually unknown, so $K$ is treated as a tunable constant. Larger $K$ makes the filter conservative: weaker recovery but less noise. Smaller $K$ makes it aggressive: stronger recovery but more noise.
A typical approach is to start with a small $K$ and gradually increase it, watching the balance between noise level and sharpness in the result. A common empirical range is $K = 0.01 \sim 0.1$ (depending on image dynamic range), but this varies entirely with the scenario.
Why Wiener Filtering Is the Most Practical
Wiener filtering is the most practical classical method for three reasons:
- Mathematically optimal: It is the optimal linear estimator in the MSE sense, with theoretical guarantees
- Computationally efficient: Pointwise operations in the frequency domain, $O(N \log N)$ complexity (FFT), far faster than spatial domain iterative methods
- Simple parameter: Only one $K$ parameter, easy to implement in engineering
Its limitation is the assumption that noise and signal are stationary random processes (power spectra do not vary with position), which real images often violate. But as a classical method, it balances effectiveness and efficiency very well.
Constrained Least Squares: Incorporating Prior Knowledge
The Idea
Wiener filtering relies on statistical assumptions (stationarity) that may not hold. Constrained Least Squares (CLS) takes a different angle: instead of statistical assumptions, use prior knowledge about image structure.
The specific prior: natural images are smooth, and gradients should not be too large. Add this constraint to the optimization objective:
$$\min_f ; | g - h * f |^2 + \lambda | \nabla f |^2$$
The first term $|g - h * f|^2$ is the data fidelity term, ensuring the recovered image matches the observation after degradation. The second term $\lambda |\nabla f|^2$ is the regularization term, penalizing large gradients to suppress the wild oscillations caused by noise amplification. $\lambda$ is the regularization parameter controlling the tradeoff.
Frequency Domain Solution
The gradient penalty is approximated using the Laplacian operator $\nabla^2$ (Laplacian is the divergence of the gradient, measuring second-order change rate). After solving in the frequency domain:
$$\hat{F}(u,v) = \frac{H^*(u,v)}{|H(u,v)|^2 + \lambda |P(u,v)|^2} G(u,v)$$
where $P(u,v)$ is the frequency domain representation of the Laplacian operator, and $\lambda$ is the regularization parameter.
The formula structure is nearly identical to Wiener filtering. The difference is the regularization term changes from constant $K$ to $\lambda |P(u,v)|^2$. $P(u,v)$ is large at high frequencies and small at low frequencies, meaning regularization strength adapts with frequency: strong constraint at high frequencies (suppressing noise), weak constraint at low frequencies (preserving signal).
Wiener vs Constrained Least Squares
The two methods share a similar formula structure but differ in design philosophy:
- Wiener filtering derives $K$ from SNR statistics, requiring estimates of noise and signal power spectra
- Constrained least squares derives $\lambda$ from image smoothness assumptions, requiring tuning of regularization strength
The choice depends on the scenario. If you have a good estimate of the noise level, Wiener filtering is more appropriate. If you have prior knowledge about image smoothness, constrained least squares is more flexible. In practice, the difference in results is usually small, and Wiener filtering is more commonly used because its parameter has a more intuitive meaning.
Comparing the Three Methods
flowchart TD
A["Blurred image G(u,v)"] --> B["Choose strategy"]
B --> C["Inverse filter<br/>F̂ = G / H"]
B --> D["Wiener filter<br/>H* / (|H|² + K)"]
B --> E["Constrained LS<br/>H* / (|H|² + λ|P|²)"]
C --> F["Optimal without noise<br/>Fails with noise"]
D --> G["MMSE optimal<br/>Needs K estimate"]
E --> H["Smoothness constraint<br/>Needs λ tuning"]
classDef input fill:#2196F3,color:#fff
classDef decision fill:#FF9800,color:#fff
classDef method fill:#9C27B0,color:#fff
classDef result fill:#4CAF50,color:#fff
class A input
class B decision
class C,D,E method
class F,G,H resultThe diagram shows the positioning of each method:
- Inverse filtering is the theoretical baseline, valid only in the noiseless case, unusable in practice, but essential for understanding the nature of degradation
- Wiener filtering is the engineering default, statistically optimal, intuitive parameter, suitable for most scenarios
- Constrained least squares offers an alternative regularization approach, suitable when you have clear priors about image structure
Summary
The fundamental tension in frequency domain restoration is this: information lost during degradation is mixed with noise, and you cannot recover signal without also amplifying noise.
Inverse filtering ignores this tension, does raw division, and explodes at high frequencies. Wiener filtering acknowledges the tension, using the $K$ parameter to adaptively balance recovery and suppression: recover at low frequencies, suppress at high frequencies. Constrained least squares approaches the same goal from a different angle, using a smoothness constraint to achieve a similar effect.
The key to understanding these three methods is not memorizing formulas, but grasping how each one handles the fundamental tradeoff between signal recovery and noise suppression. Modern methods that follow (Tikhonov regularization, variational methods, deep learning restoration) all solve the same problem, just with more sophisticated regularization.