Wiener Filter

The Wiener filter is the linear estimator of an image that minimizes the expected squared error, given a degradation model and the second-order statistics of the image and the noise. For a stationary image observed through a known blur and additive noise, it becomes a single gain per spatial frequency that depends on how much signal and how much noise that frequency carries: frequencies dominated by signal are restored, and frequencies dominated by noise are suppressed. With no blur it is a denoiser; with a blur it is a regularized deconvolution, the Wiener deconvolution, which scikit-image’s documentation calls Wiener–Hunt deconvolution (scikit-image 0.25.2 source). In estimation terms it is the linear minimum mean-square-error (LMMSE) estimate, and for Gaussian image and noise models it is also the posterior mean.

The filter is named after Norbert Wiener, and Sage et al. (2017) cite his 1949 monograph on the extrapolation, interpolation and smoothing of stationary time series for it. The Kalman filter paper recounts how Wiener’s formulation of optimal linear estimation led to the Wiener–Hopf equation and how Kalman later recast the problem in state-space form. Image processing gives the filter’s formula among the restoration problems, and image deblurring derives it as a special case of regularized deconvolution and compares it with iterative methods. This article covers the minimum-mean-square-error filter itself: where the formula comes from, the error it achieves, what happens when the spectra it needs are wrong or estimated, and the local and transform-domain filters that reuse the same estimate.

Problem

The observation model is the one used throughout restoration,

$$y = h * x + n ,$$

where $x$ is the unknown image, $h$ a known point-spread function, $n$ noise independent of $x$, and $*$ convolution. The task is to find the filter $w$ for which the estimate $\hat x = w * y$ has the smallest expected squared error $\mathrm{E}\,\lVert \hat x - x \rVert^2$, the expected mean squared error summed over pixels.

Two cases are distinguished. In denoising, $h$ is the identity and the only degradation is noise, as in image denoising. In deconvolution, $h$ is a blur, and dividing by its transfer function amplifies noise wherever the transfer function is small; the Wiener filter decides, frequency by frequency, how far to go towards that division.

The expectation is over both the noise and the image. The filter therefore treats the image as a random field and needs a statistical model of it: its mean, and its covariance or, for a stationary field, its power spectrum.

Inputs and Outputs

• Inputs: the degraded image $y$; the transfer function $H$ of the blur (or $H = 1$ for denoising); the power spectrum $S_x$ of the image; and the power spectrum $S_n$ of the noise, which is the constant $\sigma^2$ for white noise of variance $\sigma^2$.
• Output: the estimate $\hat x$, and, when the spectra are correct, the expected error of that estimate in closed form.

In practice $S_x$ is not known. It is replaced by a model with a few parameters, by an estimate from the data, or by a constant, and the filter is then no longer the minimum-error one. The noise level can come from calibrating the sensor or be estimated from the image, as image noise describes.

Intuition

At a frequency where the image has much more power than the noise, the observation is reliable, and the best estimate inverts the blur there. At a frequency where the noise dominates, the observation says little about the image, and the best estimate is close to zero, the image’s mean at that frequency. In between, the filter takes a fraction of the observation set by the ratio of the two powers. Natural images have most of their power at low frequencies (Torralba, Isola & Freeman, 2024, §27.4.1) while white noise has the same power everywhere, so, in our reading, the Wiener filter of a photograph is usually a low-pass filter whose cut-off moves down as the noise grows.

It follows from the formula below that the filter never amplifies a frequency more than the inverse filter would, and amplifies less wherever the noise matters. In the limit of no noise it becomes the inverse filter.

Algorithm

For a periodic image, where the blur is a circular convolution and the image statistics are stationary, the filter is applied in the discrete Fourier domain:

1. Compute the transfer function $H_k$ of the point-spread function at every frequency $k$, with the kernel centred at the origin.
2. Obtain the image power spectrum $S_k$ and the noise power spectrum $N_k$ (for white noise, $N_k = \sigma^2$).
3. Form the gain $W_k = H_k^{*} S_k / (|H_k|^2 S_k + N_k)$.
4. Transform the observation, $Y_k = \mathcal{F}(y)_k$, multiply, $\hat X_k = W_k Y_k$, and transform back.

The steps cost two FFTs and a pointwise product. Real images are not periodic, so implementations that apply the filter through the FFT impose a circulant blur, with the border effects that image deblurring studies. scikit-image’s wiener notes the circulant assumption in its documentation and applies exactly this filter with a regularization operator in place of $N_k/S_k$, as described under Parameters (scikit-image 0.25.2 source).

Mathematical Formulation

The minimum-error gain

With a unitary discrete Fourier transform, a stationary periodic field has uncorrelated Fourier coefficients, with $\mathrm{E}|X_k|^2 = S_k$; Torralba, Isola and Freeman (2024, §27.5) state the same fact as the circulant covariance matrix of a stationary image being diagonalized by the Fourier basis, with the power spectrum on the diagonal. The convolution theorem of the Fourier transform turns the model into $Y_k = H_k X_k + N_k'$, one equation per frequency, with noise coefficients of variance $N_k$. For any filter with gains $W_k$, the error at frequency $k$ is $X_k - W_k Y_k = (1 - W_k H_k) X_k - W_k N_k'$, and since image and noise are independent its expected power is

$$e_k(W_k) = |1 - W_k H_k|^2 S_k + |W_k|^2 N_k .$$

The total expected squared error is the sum of the $e_k$ (Parseval’s theorem for the unitary transform), and each $e_k$ depends only on its own $W_k$. Setting the derivative with respect to $W_k^{*}$ to zero gives, in our derivation,

$$W_k = \frac{H_k^{*} S_k}{|H_k|^2 S_k + N_k} = \frac{H_k^{*}}{|H_k|^2 + N_k / S_k} ,$$

the second form being the one image processing gives. The ratio $N_k/S_k$ is the noise-to-signal ratio at frequency $k$. Substituting back gives the minimum error per frequency and the minimum mean squared error per pixel of an $M$-pixel image,

$$e_k^{\min} = \frac{S_k N_k}{|H_k|^2 S_k + N_k}, \qquad \mathrm{MMSE} = \frac{1}{M} \sum_k \frac{S_k N_k}{|H_k|^2 S_k + N_k} .$$

Both formulas follow from the quadratic form above; because $e_k$ is a convex quadratic in $W_k$, every other choice of gains has a larger expected error, by $|H_k|^2 S_k + N_k$ times $|W_k - W_k^{\text{Wiener}}|^2$ at each frequency (our derivation). That last fact explains why moderate errors in the spectra cost little: the excess error is second order in the error of the gain.

The formula covers the inverse filter as the limit $N_k \to 0$, where $W_k = 1/H_k$ and the expected error is $\sum_k N_k / |H_k|^2$, unbounded as $|H_k|$ approaches zero. It covers denoising with $H_k = 1$, where

$$W_k = \frac{S_k}{S_k + \sigma^2} ,$$

the form that Torralba, Isola and Freeman (2024, §27.5.2) write for a power-law spectrum $S_k = A / |k|^{2\alpha}$.

The Bayesian view

If the image is a Gaussian random field with covariance $C$ and the noise is white Gaussian, the posterior of $x$ given $y$ is Gaussian, and its mode and its mean coincide with the linear estimate. Torralba, Isola and Freeman (2024, §27.5.2) derive the denoiser this way, as the maximum a posteriori estimate $\hat x = C (C + \sigma^2 I)^{-1} y$, and present the Gaussian model as the maximum-entropy distribution when only the correlations between pixels are constrained (§27.5). Their likelihood is printed as $\exp(-\lVert y - x \rVert^2 / \sigma^2)$, without the factor $\tfrac{1}{2}$ of Gaussian noise of variance $\sigma^2$; the closed form they give is the one for that variance (our derivation). Without the Gaussian assumption the same filter is still the best linear estimator, since $e_k$ involves only second moments; the restriction to linear estimators is what nonlinear denoisers give up.

The matrix form does not need stationarity. NL-Bayes, described under Variants, applies $\bar P + C (C + \sigma^2 I)^{-1}(\tilde P - \bar P)$ to small groups of similar patches with a covariance $C$ estimated from the group, which is the same estimate with a local, data-driven Gaussian model (Lebrun, Buades & Morel, 2013, Eq. 9).

Relation to regularized deconvolution

The regularized inverse filter minimizes $\lVert y - Hx \rVert^2 + \lambda \lVert Lx \rVert^2$ for a penalty operator $L$, and with circular boundaries its solution is the gain $H_k^{*} / (|H_k|^2 + \lambda |L_k|^2)$, as image deblurring derives. Comparing with the Wiener gain, the two coincide when $\lambda |L_k|^2 = N_k / S_k$: in our reading, a quadratic penalty is a Wiener filter for the image spectrum $S_k \propto 1 / |L_k|^2$, with $\lambda$ playing the role of the noise level. Sage et al. (2017) make the same identification: when the filtering by $L^{\mathsf{T}} L$ has a whitening effect on the image and the regularization is set by the noise, the regularized inverse filter amounts to Wiener filtering. A constant penalty ($L = I$, Tikhonov regularization) corresponds to a white image spectrum, and a Laplacian penalty, whose $|L_k|^2$ grows like $|k|^4$ at low frequencies, to a spectrum falling like $|k|^{-4}$.

Parameters

The filter has no free parameters when the spectra are known. In practice its parameters are the choices made in place of them.

• The image spectrum $S_k$. Torralba, Isola and Freeman (2024, §27.4.1) describe the Fourier amplitude of natural images as falling off roughly as a power of the radial frequency, $1/w^{\alpha}$ with $\alpha \approx 1$, and their three 512 × 512 examples as decaying roughly as $1/w^{1.5}$. Anger, Facciolo and Delbracio (2018, Eqs. 20–21) use the power-spectrum model $|\hat u(\xi)|^2 \approx (\epsilon + c^{-1} \lVert \xi \rVert^{\beta})^{-1}$, valid at all frequencies, and note that $\beta = 2$ is expected if the relative contrast energy of an image does not depend on the viewing distance. A parametric spectrum of this kind can be fitted to the observation, after subtracting the noise power and dividing by $|H|^2$, as the example below does.
• A constant noise-to-signal ratio $K$. Replacing $N_k/S_k$ by one constant $K$ gives $W_k = H_k^{*} / (|H_k|^2 + K)$. This is the Wiener filter for a white image spectrum, or Tikhonov regularization with $L = I$, and $K$ becomes a tuning parameter.
• A regularization operator and balance. scikit-image’s wiener(image, psf, balance, reg=None, is_real=True, clip=True) computes $H_k^{*} / (|H_k|^2 + \lambda |D_k|^2)$, with the balance $\lambda$ and a regularization operator $D$ that is the discrete Laplacian by default (the 3 × 3 kernel with 4 at the centre and −1 at its four neighbours) and can be passed as an impulse response or a transfer function; a transfer function passed for the point-spread function or for $D$ is reduced to its real part before use (scikit-image 0.25.2 source), which in our reading is exact only for kernels symmetric about the origin, such as those used here. Passing $D$ as the transfer function $\sqrt{N_k / S_k}$ with $\lambda = 1$ gives the Wiener filter with a known spectrum, which is how the example checks its hand-written filter. With the default clip=True the result is clipped to $[-1, 1]$, so data on other scales need clip=False.
• The noise level. For white noise the only noise parameter is $\sigma^2$. When it is unknown, scikit-image’s unsupervised_wiener estimates the noise precision and the prior precision of the Laplacian model together with the image, by a Gibbs sampler that returns the posterior mean and the two sampled chains; its defaults are a burn-in of 15 samples, at least 30 and at most 200 iterations, and a stopping threshold of $10^{-4}$ on the change of the running mean (scikit-image 0.25.2 source). The source cites Orieux, Giovannelli and Rodet (2010) for the algorithm.

Complexity

For an image of $M$ pixels, the global filter costs two FFTs, $O(M \log M)$, plus $O(M)$ to form the gains; estimating a radially averaged spectrum adds $O(M)$. unsupervised_wiener repeats an FFT-domain sampling step at every iteration, up to 200 by default. The matrix form $C(C + \sigma^2 I)^{-1}$ costs $O(d^3)$ for a $d$-dimensional covariance, which is affordable for patches but not for whole images; that is why the image-level filter relies on the Fourier diagonalization of a stationary covariance, and why patch-based methods such as NL-Bayes work with small patches.

Implementation

The example computes the Wiener filter where its assumptions hold exactly, so that the closed-form minimum error can be compared with the error actually achieved, and then replaces the true spectrum with the substitutes used in practice. The image is a 256 × 256 periodic Gaussian random field whose power spectrum is $S(f) = A / (f_0^2 + |f|^2)$ with $f_0 = 1/64$ cycles per pixel, the form of Anger et al.’s Eq. 21 with $\beta = 2$, scaled to a pixel variance of 1. It is generated by filtering white Gaussian noise with $\sqrt{S}$ in the unitary Fourier domain, so that its Fourier coefficients have exactly the expected power $S$. For deconvolution it is blurred by circular convolution with a sampled Gaussian point-spread function ($\sigma_b = 1.5$ pixels, 13 × 13, summing to one), whose transfer function is close to zero at the highest frequencies; for denoising $H = 1$. White Gaussian noise of standard deviation 0.05, 0.5 or 2.0 is added, giving six conditions, each run with seeds 0 to 4; a seed fixes the field and a unit noise pattern that is scaled to each noise level.

Seven methods are compared:

• the inverse filter $Y/H$, which in denoising is the identity;
• the Wiener filter with the true spectrum;
• the Wiener filter with an estimated spectrum: a power law $e^{a} |f|^{-\beta}$ fitted by least squares on a log–log scale to the radially averaged periodogram of the observation, minus $\sigma^2$, divided by $|H|^2$, using only the radial bins whose power exceeds $2\sigma^2$, with the noise level known;
• a constant noise-to-signal ratio $K = \sigma^2 / \sigma_x^2$, the right ratio of total variances with a white spectrum;
• a constant $K$ at its best value, found by minimizing the closed-form expected error over $K$, an oracle no real data would allow;
• scikit-image’s default Laplacian penalty at its best balance, found the same way;
• scikit-image’s unsupervised_wiener, which estimates its own noise and prior levels, run with rng set to the seed.

For each method the code prints the mean squared error divided by the closed-form MMSE of the condition, as the mean and range over the five seeds. For the methods that are fixed linear filters it also prints the closed-form ratio and the ratio of the empirical error to that method’s own closed form, and the number of seeds on which the true-spectrum Wiener filter had no larger error. For the estimated spectrum it prints the fitted exponent, and for unsupervised_wiener its noise estimate divided by the true $\sigma$, computed as one over the square root of the mean noise precision sampled in the iterations whose image samples enter its posterior mean. The first lines check the data (the periodogram of each field against $S$, and the standard deviation of each noise pattern), the inverse filter on noise-free blurred data, and finally the agreement between the hand-written filters and scikit-image’s wiener.

The checks and their tolerances were fixed before the example was first run, from closed forms computed on the same set-up. In our derivation, the empirical squared error of a fixed linear filter is a sum of independent terms with relative standard deviation $\sqrt{2 \sum_k e_k^2} / \sum_k e_k$, between 0.0055 and 0.065 here, and each ratio to its closed form was required to lie within four of these standard deviations of 1. The data checks used the same rule (tolerances 0.022 for the periodogram ratio and 0.011 for the noise standard deviation), the noise-free inverse filter had to reach an RMSE below $10^{-6}$, and the library comparison a maximum difference below $10^{-9}$.

The plan also recorded five predictions. Three follow from the formulas: every fixed linear filter has a closed-form error above the MMSE, and the true-spectrum filter was predicted to be no worse than each of them on every realization; in denoising, the best constant $K$ equals $\sigma^2/\sigma_x^2$ (with $H = 1$ the error of the constant gain $1/(1+K)$ is $(K^2 \sigma_x^2 + \sigma^2)/(1+K)^2$, minimized at that value), while in deblurring the best constant loses more as the noise grows; and the inverse filter is catastrophic in deblurring and equals doing nothing in denoising. Two were predictions without a derivation: that the estimated spectrum stays within 2% of the MMSE in every condition and seed, and that unsupervised_wiener stays above the best Laplacian filter’s closed-form ratio (less four standard deviations) and estimates the noise level within 10%. The example uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 on one thread, runs in about two seconds, and printed identical output on two runs. All settings were fixed before the first run, and nothing was changed after it.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"                 # one BLAS thread, before NumPy
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from scipy.optimize import minimize_scalar
from skimage.restoration import wiener, unsupervised_wiener

N, f0, sb, R = 256, 1 / 64, 1.5, 6
fy, fx = np.meshgrid(np.fft.fftfreq(N), np.fft.fftfreq(N), indexing="ij")
fr = np.sqrt(fx**2 + fy**2)                               # radial frequency, cycles/px
S = 1.0 / (f0**2 + fr**2); S /= S.mean()                  # power spectrum, pixel variance 1

r = np.arange(-R, R + 1)                                  # Gaussian PSF, sigma 1.5 px, 13 x 13
psf = np.exp(-(r[:, None]**2 + r[None, :]**2) / (2 * sb**2)); psf /= psf.sum()
def tf(kernel):                                           # transfer function, kernel centred at 0
    pad = np.zeros((N, N)); c = kernel.shape[0] // 2
    pad[:kernel.shape[0], :kernel.shape[1]] = kernel
    return np.fft.fft2(np.roll(pad, (-c, -c), axis=(0, 1)))
lap = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], float)
L2 = np.abs(tf(lap))**2                                   # skimage's default penalty
F = lambda a: np.fft.fft2(a, norm="ortho")
iF = lambda A: np.fft.ifft2(A, norm="ortho").real
half = lambda A: A[:, :N // 2 + 1]                        # rfft layout for skimage

def closed(W, H, s2):                                     # expected MSE of a fixed linear filter
    e = np.abs(1 - W * H)**2 * S + np.abs(W)**2 * s2
    return e.mean()

def oracle(make, H, s2):                                  # minimize the closed form over log parameter
    t = minimize_scalar(lambda t: closed(make(np.exp(t)), H, s2), bounds=(-30, 5),
                        method="bounded", options={"xatol": 1e-6}).x
    return np.exp(t)

rbin = np.rint(fr * N).astype(int)
def estimated_wiener(Y, H, s2):                           # power law fitted to the observation
    P = np.abs(Y)**2
    rs, ss = [], []
    for b in range(1, N // 2):
        m = rbin == b
        Pb = P[m].mean()
        if Pb > 2 * s2:
            rs.append(b / N); ss.append((Pb - s2) / np.mean(np.abs(H[m])**2))
    slope, a = np.polyfit(np.log(rs), np.log(ss), 1)
    Sh = np.exp(a) * np.maximum(fr, 1 / N)**slope
    return np.conj(H) * Sh / (np.abs(H)**2 * Sh + s2), -slope

seeds, sigmas = range(5), (0.05, 0.5, 2.0)
data = {}
print("C1 data check  seed: periodogram/S   noise std   (both within 1 +/- 0.022, 0.011)")
for s in seeds:
    rng = np.random.default_rng(s)
    X = np.sqrt(S) * F(rng.standard_normal((N, N)))
    n1 = rng.standard_normal((N, N))
    data[s] = (X, n1)
    m = np.ones((N, N), bool); m[0, 0] = False
    print(f"  seed {s}: {np.mean(np.abs(X[m])**2 / S[m]):.4f}   {n1.std():.4f}")

Hb = tf(psf)
x0 = iF(data[0][0])
xinv = iF(F(iF(Hb * data[0][0])) / Hb)
print(f"C2 noise-free deblurring, inverse filter RMSE {np.sqrt(np.mean((xinv - x0)**2)):.2e} (< 1e-6)")

names = ["inverse", "Wiener, true S", "Wiener, estimated S", "K = s^2/var", "K oracle",
         "Laplacian oracle", "unsupervised_wiener"]
libdiff = 0.0
for case, H, kern in [("deblurring", Hb, psf), ("denoising", np.ones((N, N)), np.ones((1, 1)))]:
    for sig in sigmas:
        s2 = sig**2
        Wt = np.conj(H) * S / (np.abs(H)**2 * S + s2)
        mmse = closed(Wt, H, s2)
        K = oracle(lambda k: np.conj(H) / (np.abs(H)**2 + k), H, s2)
        lam = oracle(lambda l: np.conj(H) / (np.abs(H)**2 + l * L2), H, s2)
        fixed = {"inverse": 1 / H, "Wiener, true S": Wt,
                 "K = s^2/var": np.conj(H) / (np.abs(H)**2 + s2),
                 "K oracle": np.conj(H) / (np.abs(H)**2 + K),
                 "Laplacian oracle": np.conj(H) / (np.abs(H)**2 + lam * L2)}
        res = {n: [] for n in names}; emp_cf = {n: [] for n in fixed}; betas, shat = [], []
        for s in seeds:
            X, n1 = data[s]
            x = iF(X)
            y = iF(H * X) + sig * n1
            Y = F(y)
            mse = lambda xh: np.mean((xh - x)**2)
            for n, W in fixed.items():
                v = mse(iF(W * Y)); res[n].append(v / mmse); emp_cf[n].append(v / closed(W, H, s2))
            lib = {"Wiener, true S": wiener(y, kern, 1.0, reg=half(np.sqrt(s2 / S)).astype(complex), clip=False),
                   "K = s^2/var": wiener(y, kern, s2, reg=half(np.ones((N, N))).astype(complex), clip=False),
                   "K oracle": wiener(y, kern, K, reg=half(np.ones((N, N))).astype(complex), clip=False),
                   "Laplacian oracle": wiener(y, kern, lam, clip=False)}
            for n, xl in lib.items():
                libdiff = max(libdiff, np.abs(xl - iF(fixed[n] * Y)).max())
            We, beta = estimated_wiener(Y, H, s2)
            res["Wiener, estimated S"].append(mse(iF(We * Y)) / mmse); betas.append(beta)
            xu, chains = unsupervised_wiener(y, kern, clip=False, rng=s)
            res["unsupervised_wiener"].append(mse(xu) / mmse)
            shat.append(1 / np.sqrt(np.mean(chains["noise"][17:])) / sig)
        print(f"\n{case}, noise sigma {sig}: MMSE {mmse:.4e}   (K oracle {K:.3e}, Laplacian lambda {lam:.3e})")
        print(f"  {'method':<21} MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=")
        for n in names:
            v = np.array(res[n])
            fmt = (lambda z: f"{z:.4e}") if v.max() > 100 else (lambda z: f"{z:.4f}")
            line = f"  {n:<21} {fmt(v.mean())} [{fmt(v.min())}, {fmt(v.max())}]"
            if n in fixed:
                c = closed(fixed[n], H, s2) / mmse
                le = sum(a <= b for a, b in zip(res["Wiener, true S"], res[n]))
                line = f"{line:<58} {fmt(c):<19} [{min(emp_cf[n]):.4f}, {max(emp_cf[n]):.4f}]   {le}/5"
            print(line)
        print(f"  fitted spectral exponent beta {min(betas):.3f}-{max(betas):.3f} (true spectrum: 2 at high f)")
        print(f"  unsupervised noise estimate / sigma {min(shat):.4f}-{max(shat):.4f}")
print(f"\nC3 max |hand - skimage wiener| over fixed filters: {libdiff:.1e} (< 1e-9)")

C1 data check  seed: periodogram/S   noise std   (both within 1 +/- 0.022, 0.011)
  seed 0: 0.9989   1.0031
  seed 1: 0.9919   1.0008
  seed 2: 0.9940   1.0028
  seed 3: 0.9952   1.0007
  seed 4: 0.9988   0.9979
C2 noise-free deblurring, inverse filter RMSE 5.71e-09 (< 1e-6)

deblurring, noise sigma 0.05: MMSE 2.2697e-01   (K oracle 3.921e-03, Laplacian lambda 7.471e-04)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               1.5174e+13 [1.4470e+13, 1.5869e+13] 1.5288e+13          [0.9465, 1.0380]   5/5
  Wiener, true S        0.9951 [0.9892, 0.9980]            1.0000              [0.9892, 0.9980]   5/5
  Wiener, estimated S   0.9986 [0.9933, 1.0018]
  K = s^2/var           1.0118 [1.0060, 1.0151]            1.0159              [0.9903, 0.9992]   5/5
  K oracle              0.9978 [0.9923, 1.0004]            1.0025              [0.9899, 0.9979]   5/5
  Laplacian oracle      0.9959 [0.9896, 0.9991]            1.0009              [0.9888, 0.9983]   5/5
  unsupervised_wiener   1.0359 [1.0281, 1.0400]
  fitted spectral exponent beta 1.543-1.622 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 1.0174-1.0224

deblurring, noise sigma 0.5: MMSE 3.8060e-01   (K oracle 1.291e-01, Laplacian lambda 1.605e-01)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               9.0493e+14 [8.6292e+14, 9.4637e+14] 9.1170e+14          [0.9465, 1.0380]   5/5
  Wiener, true S        0.9976 [0.9868, 1.0086]            1.0000              [0.9868, 1.0086]   5/5
  Wiener, estimated S   1.0061 [0.9966, 1.0148]
  K = s^2/var           1.0872 [1.0759, 1.0952]            1.0926              [0.9847, 1.0024]   5/5
  K oracle              1.0319 [1.0202, 1.0435]            1.0349              [0.9858, 1.0084]   5/5
  Laplacian oracle      1.0024 [0.9921, 1.0128]            1.0049              [0.9873, 1.0079]   5/5
  unsupervised_wiener   1.0325 [1.0238, 1.0413]
  fitted spectral exponent beta 1.361-1.458 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 1.0076-1.0139

deblurring, noise sigma 2.0: MMSE 5.6774e-01   (K oracle 7.082e-01, Laplacian lambda 8.093e+00)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               9.7063e+15 [9.2557e+15, 1.0151e+16] 9.7789e+15          [0.9465, 1.0380]   5/5
  Wiener, true S        0.9973 [0.9875, 1.0113]            1.0000              [0.9875, 1.0113]   5/5
  Wiener, estimated S   1.0249 [1.0139, 1.0442]
  K = s^2/var           1.4297 [1.4161, 1.4732]            1.4588              [0.9707, 1.0099]   5/5
  K oracle              1.1236 [1.1116, 1.1376]            1.1319              [0.9821, 1.0050]   5/5
  Laplacian oracle      1.0162 [1.0044, 1.0299]            1.0174              [0.9873, 1.0123]   5/5
  unsupervised_wiener   1.0402 [1.0292, 1.0517]
  fitted spectral exponent beta 0.997-1.223 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 0.9997-1.0055

denoising, noise sigma 0.05: MMSE 2.4769e-03   (K oracle 2.500e-03, Laplacian lambda 4.392e-04)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               1.0115 [1.0051, 1.0157]            1.0093              [0.9958, 1.0063]   5/5
  Wiener, true S        1.0022 [0.9964, 1.0067]            1.0000              [0.9964, 1.0067]   5/5
  Wiener, estimated S   1.0024 [0.9965, 1.0068]
  K = s^2/var           1.0091 [1.0030, 1.0134]            1.0068              [0.9962, 1.0066]   5/5
  K oracle              1.0091 [1.0030, 1.0134]            1.0068              [0.9962, 1.0066]   5/5
  Laplacian oracle      1.0028 [0.9968, 1.0073]            1.0006              [0.9962, 1.0067]   5/5
  unsupervised_wiener   38.1363 [37.2221, 38.9056]
  fitted spectral exponent beta 1.722-1.774 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 8.0108-8.1406

denoising, noise sigma 0.5: MMSE 1.4197e-01   (K oracle 2.500e-01, Laplacian lambda 5.613e-02)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               1.7647 [1.7535, 1.7720]            1.7609              [0.9958, 1.0063]   5/5
  Wiener, true S        1.0010 [0.9957, 1.0043]            1.0000              [0.9957, 1.0043]   5/5
  Wiener, estimated S   1.0075 [1.0026, 1.0110]
  K = s^2/var           1.4057 [1.3996, 1.4187]            1.4087              [0.9935, 1.0071]   5/5
  K oracle              1.4057 [1.3996, 1.4187]            1.4087              [0.9935, 1.0071]   5/5
  Laplacian oracle      1.0238 [1.0207, 1.0263]            1.0227              [0.9980, 1.0035]   5/5
  unsupervised_wiener   1.5825 [1.5642, 1.6154]
  fitted spectral exponent beta 1.680-1.734 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 1.3528-1.3660

denoising, noise sigma 2.0: MMSE 4.6796e-01   (K oracle 4.000e+00, Laplacian lambda 3.433e+00)
  method                MSE/MMSE mean [min, max]      closed-form ratio   empirical/closed [min, max]  seeds M2 <=
  inverse               8.5660 [8.5116, 8.6016]            8.5477              [0.9958, 1.0063]   5/5
  Wiener, true S        0.9981 [0.9877, 1.0084]            1.0000              [0.9877, 1.0084]   5/5
  Wiener, estimated S   1.1834 [1.1332, 1.2713]
  K = s^2/var           1.6782 [1.6609, 1.7290]            1.7095              [0.9715, 1.0114]   5/5
  K oracle              1.6782 [1.6609, 1.7290]            1.7095              [0.9715, 1.0114]   5/5
  Laplacian oracle      1.0577 [1.0473, 1.0651]            1.0586              [0.9894, 1.0062]   5/5
  unsupervised_wiener   1.0844 [1.0738, 1.0910]
  fitted spectral exponent beta 1.163-1.315 (true spectrum: 2 at high f)
  unsupervised noise estimate / sigma 1.0325-1.0381

C3 max |hand - skimage wiener| over fixed filters: 4.9e-15 (< 1e-9)

All five checks passed. The periodogram ratios (0.9919 to 0.9989) and noise standard deviations (0.9979 to 1.0031) are within their tolerances, the inverse filter recovers the noise-free field to an RMSE of $5.7 \times 10^{-9}$, in our reading round-off amplified by the near-zero transfer-function values, and the hand-written filters agree with scikit-image’s wiener to $5 \times 10^{-15}$. Every fixed filter’s empirical error lies within its tolerance of its closed form, the true-spectrum Wiener filter between 0.9868 and 1.0113 of the MMSE.

The closed form predicts the error. The MMSE is an expectation, and a single realization scatters around it by about one percent here; the ratio of the true-spectrum filter is below 1 for some seeds, because those realizations happened to be easier, not because the filter beat the bound. Across conditions, the MMSE expresses how much information survives. In denoising at $\sigma = 0.05$ it is 99% of the noise variance; with $H = 1$ the MMSE is $\sigma^2$ times the average of $S_k/(S_k + \sigma^2)$, so a value this close to $\sigma^2$ means that the field dominates the noise at almost every frequency and the best filter barely changes the image; at $\sigma = 2$ it is 12% of the noise variance. Deblurring at $\sigma = 0.05$ leaves an error of 0.23, almost a quarter of the field’s variance, because the frequencies that the blur attenuates below the noise level cannot be recovered.

The true spectrum wins on every realization. The first prediction held in all six conditions: every fixed filter’s closed-form ratio is above 1, and the true-spectrum filter had no larger error than each fixed alternative on all five seeds, including the conditions where the Laplacian filter’s expected excess is below 0.1%. The prediction for the inverse filter held as well. Its error in deblurring is $10^{13}$ to $10^{16}$ times the MMSE; its closed form is the average of $\sigma^2/|H_k|^2$, and multiplying the printed closed-form ratio by the MMSE and dividing by $\sigma^2$ gives an average of $1/|H_k|^2$ of about $1.4 \times 10^{15}$ at every noise level, so the transfer function falls to $3 \times 10^{-8}$ or below. In denoising its error is the noise variance itself, from 1.01 times the MMSE at $\sigma = 0.05$ to 8.6 times at $\sigma = 2$.

What a wrong spectrum costs. A constant noise-to-signal ratio treats the image as white. In denoising a constant ratio can only scale the whole image, and the rows for $K = \sigma^2/\sigma_x^2$ and the oracle $K$ are identical, as derived; the oracle values printed in the headers equal $\sigma^2$. The prediction for the constant ratio held in all six conditions: the two rows coincide in denoising, and in deblurring the oracle constant’s cost grows with the noise. The cost of that scaling grows with the noise, to 41% and 68% above the MMSE at $\sigma = 0.5$ and 2. In deblurring the best constant differs from $\sigma^2$, and even at its best its average error is 3% and 12% above the MMSE at $\sigma = 0.5$ and 2, close to its closed forms; at $\sigma = 0.05$ its expected excess is only 0.25%, but it still lost to the Wiener filter on every seed. The Laplacian penalty, whose implied spectrum falls like $|f|^{-4}$ instead of $|f|^{-2}$, does much better at its best balance: its average error is within 2.4% of the MMSE except in denoising at $\sigma = 2$, where it is 5.8% above. In our reading, a decaying penalty gets the low-pass shape right and only misplaces the transition, and the second-order sensitivity of the error to the gain keeps that cheap; the oracle choices of $K$ and $\lambda$ are not available in practice, so these rows are lower bounds for their families.

The estimated spectrum failed its prediction at the highest noise. At $\sigma = 0.05$ and 0.5 the fitted power law gives errors within 1.5% of the MMSE, and the prediction held in those four conditions. At $\sigma = 2$ the prediction of at most 2% failed in both conditions: in deblurring the ratios range from 1.0139 to 1.0442, above 2% on some seeds, and in denoising from 1.1332 to 1.2713, above 2% on every seed. The printed exponents show why the fit goes wrong: they range from 1.0 to 1.8, never 2. In our reasoning from the construction, the local log-slope of $A/(f_0^2 + f^2)$ is $-2f^2/(f_0^2 + f^2)$, which is $-1$ at $f = f_0$ and $-1.6$ at $2 f_0$, so the spectrum looks like a power law with exponent 2 only well above $f_0 = 1/64$. As the noise grows, only low-frequency bins clear the threshold of $2\sigma^2$, the fitted exponent becomes shallower, and the extrapolated spectrum overestimates the signal at high frequencies, where the filter then lets noise through. The planning argument, that the fitted law would be wrong only at low frequencies where the signal dominates anyway, was wrong at $\sigma = 2$. A spectrum model must match the image where the filter’s gain is decided, near the frequency where signal and noise powers cross, not just where it is easy to fit.

The unsupervised filter. unsupervised_wiener never went below the bound set by the best Laplacian filter: that prediction held in all six conditions, since its smallest ratio is above the Laplacian filter’s closed-form ratio in every one. When deblurring it came within 2.4 to 5.2% of the MMSE and estimated the noise level within 2.3%, so the noise-estimate prediction held in the three deblurring conditions. The prediction that its noise estimate would be within 10% failed in two of the three denoising conditions: at $\sigma = 0.05$ it estimated about eight times the true noise level and its error was 38 times the MMSE, and at $\sigma = 0.5$ it overestimated the noise by 35 to 37% with an error 1.58 times the MMSE; at $\sigma = 2$ it was within 4%, and the prediction held. We have not tested why; one untested conjecture is that a Laplacian prior, which expects little high-frequency power, attributes the field’s high-frequency power to noise when there is no blur to remove it.

Properties and Behavior

• Optimal among linear filters, for the assumed statistics. The filter minimizes the expected squared error over all linear estimators given $S_x$ and $S_n$; for Gaussian image and noise it is also the posterior mean, and it follows that it then minimizes the expected squared error over all estimators. The guarantee is about the average over the model’s images, not about any particular image.
• Graceful degradation with wrong spectra. Because the excess error is quadratic in the error of the gain, moderate misspecification costs little, as the Laplacian rows of the example show; misspecification at the frequencies where signal and noise powers cross costs most, as the estimated spectrum at $\sigma = 2$ shows.
• Never worse than doing nothing, in expectation. In denoising, the identity is one of the linear filters, so the expected error of the Wiener filter is at most $\sigma^2$, and the MMSE in units of $\sigma^2$ measures the achievable gain. In deconvolution the zero filter, with error equal to the image variance, bounds it from above.
• Smoothing of edges. A stationary Gaussian model knows only how correlation falls with distance. Torralba, Isola and Freeman (2024, §27.5.2) find the denoising result of the Gaussian model far from satisfactory, with residual image structure in the estimated noise, and note that the model does not capture flat regions, edges and lines (§27.6). Image deblurring reports Szeliski’s view that the filter is almost never used in practice any more, having been replaced by better-performing nonlinear filters.
• Ringing and borders. As a linear frequency-domain filter applied through the FFT, Wiener deconvolution assumes periodic images, and its treatment of borders and the ringing near edges are those described in image deblurring.

Limitations

• Second-order statistics only. The filter cannot use the non-Gaussian, sharply peaked distribution of band-pass filter outputs that Torralba, Isola and Freeman (2024, §27.6) show for natural images and turn into a coring denoiser, nor the patch similarity that non-local means uses (§27.7.1).
• Stationarity. One spectrum for the whole image ignores that images change from textured to flat regions. The local and transform-domain variants below estimate the statistics locally to address this.
• Known spectra. The optimality depends on spectra that are rarely known. Estimating them from the same data introduces errors that the closed form does not account for, and the example shows such an estimate failing when a simple model is extrapolated from the frequencies where it can be fitted.
• Linear, shift-invariant degradation and additive noise. Signal-dependent Poisson noise, spatially varying blur and nonlinear camera responses all violate the model; image deblurring discusses them.

Variants

• Parametric (constant-ratio) Wiener filter. The gain $H^{*}/(|H|^2 + K)$ with a hand-tuned $K$ needs no spectrum; it is Tikhonov-regularized deconvolution, and a frequency-dependent penalty such as the Laplacian gives the regularized inverse filter that image deblurring studies with scikit-image’s wiener.
• Unsupervised Wiener–Hunt deconvolution. scikit-image’s unsupervised_wiener samples the noise and prior precisions together with the image by Gibbs sampling and returns the posterior mean (scikit-image 0.25.2 source), following Orieux, Giovannelli and Rodet (2010).
• Local adaptive Wiener filter. SciPy’s scipy.signal.wiener estimates the local mean $\mu$ and variance $v$ in a small window (3 pixels per dimension by default) and returns $\mu + (1 - \sigma^2/v)(y - \mu)$, or $\mu$ where $v < \sigma^2$; when the noise power is not given, it takes the mean of the local variances as $\sigma^2$ (SciPy 1.16.3 source). Its documentation describes it as similar to MATLAB’s wiener2 and cites Lim’s 1990 textbook. This is the pointwise Wiener gain with a locally estimated signal variance. For the multiplicative speckle of radar images, Kuan’s filter implements the local LMMSE estimate and Lee’s filter the same estimate up to a factor, as synthetic aperture radar describes.
• Empirical Wiener filtering in a transform domain. Lebrun, Buades and Morel (2013) recall that early denoisers applied linear Wiener filters to the Fourier transform, and that these filters are more efficient when applied locally to the discrete cosine transform, after Yaroslavsky and colleagues. The second step of BM3D uses its first estimate as an oracle to compute Wiener coefficients for groups of similar patches in a 3D transform domain (Lebrun, Buades & Morel, 2013, Section 3.1); non-local means and image denoising outline it.
• NL-Bayes. Lebrun, Buades and Morel (2013) model each group of similar patches as a Gaussian vector with an empirical mean and covariance. The first step estimates a patch as $\bar{\tilde P} + (C_{\tilde P} - \sigma^2 I) C_{\tilde P}^{-1} (\tilde P - \bar{\tilde P})$ from the noisy group’s covariance (Eq. 8), and the second as $\bar{\tilde P}^{\text{basic}} + C^{\text{basic}}_{\tilde P} (C^{\text{basic}}_{\tilde P} + \sigma^2 I)^{-1} (\tilde P - \bar{\tilde P}^{\text{basic}})$ with statistics from the first estimate (Eq. 9). Replacing the first estimate by the noise-free image gives an ideal Wiener filter, which in their Table 11 is 1.4 to 2.7 dB better in PSNR than the algorithm, depending on the noise level; the text describes the room for improvement as 2 to 3 dB.
• Burst merging. The HDR+ burst pipeline merges each frame into a reference with a pairwise Wiener-like filter in the frequency domain of each tile, as computational photography describes.
• Recursive estimation. For signals that evolve in time with a state-space model, the Kalman filter computes the minimum-mean-square-error linear estimate recursively, the reformulation the Kalman filter paper introduced.