Wiener Filter: Optimal Linear Filtering Theory and Python Implementation

From Wiener filter derivation (MSE minimization, Wiener-Hopf equation) to frequency-domain implementation for signal and image denoising.

Introduction

The Wiener filter is the optimal linear filter that minimizes mean squared error (MSE). For estimating a signal from noisy observations, it provides the theoretically optimal solution within the framework of stationary stochastic processes.

Problem Formulation

Noisy observation:

\[y(t) = s(t) + n(t) \tag{1}\]

where \(s(t)\) is the desired signal and \(n(t)\) is noise. We estimate \(s(t)\) by filtering:

\[\hat{s}(t) = \int_{-\infty}^{\infty} h(\tau) y(t - \tau) d\tau \tag{2}\]

MSE Minimization

Minimize the mean squared estimation error:

\[\min_h E\left[|s(t) - \hat{s}(t)|^2\right] \tag{3}\]

Wiener-Hopf Equation

The optimality condition (orthogonality principle) yields:

\[R_{yy}(\tau) * h(\tau) = R_{sy}(\tau) \tag{4}\]

where \(R_{yy}\) is the autocorrelation of the observed signal and \(R_{sy}\) is the cross-correlation between desired and observed signals.

Frequency Domain Solution

In the frequency domain, convolution becomes multiplication:

\[H(f) = \frac{S_{sy}(f)}{S_{yy}(f)} \tag{5}\]

When signal and noise are uncorrelated:

\[H(f) = \frac{S_{ss}(f)}{S_{ss}(f) + S_{nn}(f)} = \frac{\text{SNR}(f)}{\text{SNR}(f) + 1} \tag{6}\]

where \(\text{SNR}(f) = S_{ss}(f) / S_{nn}(f)\).

Intuition: Frequencies with high SNR pass through; frequencies with low SNR are attenuated.

Python Implementation

1D Signal Denoising

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

np.random.seed(42)

# --- Signal generation ---
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise

# --- Wiener filter (frequency domain) ---
N = len(observed)
Y = fft(observed)
freqs = fftfreq(N, 1/fs)

# Power spectrum estimation
S_yy = np.abs(Y)**2 / N
noise_power = np.var(noise) * N
S_nn = noise_power / N * np.ones(N)
S_ss = np.maximum(S_yy - S_nn, 0)

# Wiener filter
H = S_ss / (S_ss + S_nn + 1e-10)
S_filtered = ifft(Y * H).real

# --- Visualization ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(t, signal, 'b-', alpha=0.7, label='Original')
axes[0].set_title('Original Signal')
axes[0].set_ylabel('Amplitude')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, observed, 'gray', alpha=0.5, label='Noisy')
axes[1].set_title('Observed (Signal + Noise)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, signal, 'b-', alpha=0.3, label='Original')
axes[2].plot(t, S_filtered, 'r-', alpha=0.8, label='Wiener filtered')
axes[2].set_title('Wiener Filter Output')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

mse_noisy = np.mean((signal - observed)**2)
mse_filtered = np.mean((signal - S_filtered)**2)
print(f"MSE (noisy):    {mse_noisy:.4f}")
print(f"MSE (filtered): {mse_filtered:.4f}")
print(f"Improvement:    {(1 - mse_filtered/mse_noisy)*100:.1f}%")

Relationship with Kalman Filter

FeatureWiener FilterKalman Filter
AssumptionStationary processHandles non-stationary
ProcessingBatch (frequency domain)Recursive (time domain)
OptimalityMSE-optimal (stationary)MSE-optimal (linear Gaussian)
ApplicationSignal processingState estimation, control

The Wiener filter is optimal for stationary processes; the Kalman filter extends this to non-stationary, recursive settings.

References

  • Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 2-3.
  • Wiener, N. (1949). Extrapolation, Interpolation, and Smoothing of Stationary Time Series. MIT Press.
  • Vaseghi, S. V. (2008). Advanced Digital Signal Processing and Noise Reduction (4th ed.). Wiley.