Generalized Wiener filter

The Wiener filter as originally proposed by Norbert Wiener is a signal processing filter which uses knowledge of the statistical properties of both the signal and the noise to reconstruct an optimal estimate of the signal from a noisy one-dimensional time-ordered data stream. The generalized Wiener filter generalizes the same idea beyond the domain of one-dimensional time-ordered signal processing, with two-dimensional image processing being the most common application.

Description
Consider a data vector $$d$$ which is the sum of independent signal and noise vectors $$d = s+n$$ with zero mean and covariances $$\langle ss^T\rangle=S$$ and $$\langle nn^T\rangle=N$$. The generalized Wiener Filter is the linear operator $$G$$ which minimizes the expected residual between the estimated signal and the true signal, $$e = \langle(Gd-s)^T(Gd-s)\rangle$$. The $$G$$ that minimizes this is $$G = S(S+N)^{-1}$$, resulting in the Wiener estimator $$\hat s = S(S+N)^{-1}d$$. In the case of Gaussian distributed signal and noise, this estimator is also the maximum a posteriori estimator.

The generalized Wiener filter approaches 1 for signal-dominated parts of the data, and S/N for noise-dominated parts.

An often-seen variant expresses the filter in terms of inverse covariances. This is mathematically equivalent, but avoids excessive loss of numerical precision in the presence of high-variance modes. In this formulation, the generalized Wiener filter becomes $$G = (S^{-1}+N^{-1})^{-1}N^{-1}$$ using the identity $$A^{-1}+B^{-1}=A^{-1}(A+B)B^{-1}$$.

An example
The cosmic microwave background (CMB) is a homogeneous and isotropic random field, and its covariance is therefore diagonal in a spherical harmonics basis. Any given observation of the CMB will be noisy, with the noise typically having different statistical properties than the CMB. It could for example be uncorrelated in pixel space. The generalized Wiener filter exploits this difference in behavior to isolate as much as possible of the signal from the noise.

The Wiener-filtered estimate of the signal (the CMB in this case) $$\hat s = S(S+N)^{-1}d$$ requires the inversion of the usually huge matrix $$S+N$$. If S and N were diagonal in the same basis this would be trivial, but often, as here, that isn't the case. The solution must in these cases be found by solving the equivalent equation $$(S+N)S^{-1}\hat s = d$$, for example via conjugate gradients iteration. In this case all the multiplications can be performed in the appropriate basis for each matrix, avoiding the need to store or invert more than their diagonal. The result can be seen in the figure.