Wiener filter

In signal processing, the Wiener filter is a filter proposed by Norbert Wiener during the 1940s and published in 1949.^[1] Its purpose is to reduce the amount of noise present in a signal by comparison with an estimation of the desired noiseless signal.

Description

The goal of the Wiener filter is to filter out noise that has corrupted a signal. It is based on a statistical approach.

Typical filters are designed for a desired frequency response. The Wiener filter approaches filtering from a different angle. One is assumed to have knowledge of the spectral properties of the original signal and the noise, and one seeks the LTI filter whose output would come as close to the original signal as possible. Wiener filters are characterized by the following:^[2]

1. Assumption: signal and (additive) noise are stationary linear stochastic processes with known spectral characteristics or known autocorrelation and cross-correlation
2. Requirement: the filter must be physically realizable, i.e. causal (this requirement can be dropped, resulting in a non-causal solution)
3. Performance criteria: minimum mean-square error

This filter is frequently used in the process of deconvolution; for this application, see Wiener deconvolution.

Model/problem setup

The input to the Wiener filter is assumed to be a signal, Failed to parse (Missing texvc executable; please see math/README to configure.): s(t) , corrupted by additive noise, Failed to parse (Missing texvc executable; please see math/README to configure.): n(t) . The output, Failed to parse (Missing texvc executable; please see math/README to configure.): \hat{s}(t) , is calculated by means of a filter, Failed to parse (Missing texvc executable; please see math/README to configure.): g(t) , using the following convolution:

Failed to parse (Missing texvc executable; please see math/README to configure.): \hat{s}(t) = g(t) * (s(t) + n(t))

where

• Failed to parse (Missing texvc executable; please see math/README to configure.): s(t)

is the original signal (to be estimated)

• Failed to parse (Missing texvc executable; please see math/README to configure.): n(t)

is the noise

• Failed to parse (Missing texvc executable; please see math/README to configure.): \hat{s}(t)

is the estimated signal (which we hope will equal Failed to parse (Missing texvc executable; please see math/README to configure.): s(t)

)

• Failed to parse (Missing texvc executable; please see math/README to configure.): g(t)

is the Wiener filter

The error is Failed to parse (Missing texvc executable; please see math/README to configure.): e(t) = s(t + \alpha) - \hat{s}(t)

and the squared error is Failed to parse (Missing texvc executable; please see math/README to configure.): e^2(t) = s^2(t + \alpha) - 2s(t + \alpha)\hat{s}(t) + \hat{s}^2(t)
where

• Failed to parse (Missing texvc executable; please see math/README to configure.): s(t + \alpha)

is the desired output of the filter

• Failed to parse (Missing texvc executable; please see math/README to configure.): e(t)

is the error

Depending on the value of α the problem name can be changed:

• If Failed to parse (Missing texvc executable; please see math/README to configure.): \alpha > 0

then the problem is that of prediction

• If Failed to parse (Missing texvc executable; please see math/README to configure.): \alpha = 0

then the problem is that of filtering

• If Failed to parse (Missing texvc executable; please see math/README to configure.): \alpha < 0

then the problem is that of smoothing

Writing Failed to parse (Missing texvc executable; please see math/README to configure.): \hat{s}(t)

as a convolution integral: Failed to parse (Missing texvc executable; please see math/README to configure.): \hat{s}(t) = \int_{-\infty}^{\infty}{g(\tau)\left[s(t - \tau) + n(t - \tau)\right]\,d\tau}

.

Taking the expected value of the squared error results in

Failed to parse (Missing texvc executable; please see math/README to configure.): E(e^2) = R_s(0) - 2\int_{-\infty}^{\infty}{g(\tau)R_{x\,s}(\tau + \alpha)\,d\tau} + \int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{g(\tau)g(\theta)R_x(\tau - \theta)\,d\tau}\,d\theta}

where

• Failed to parse (Missing texvc executable; please see math/README to configure.): x(t) = s(t) + n(t)

is the observed signal

• Failed to parse (Missing texvc executable; please see math/README to configure.): \,\!R_s

is the autocorrelation function of Failed to parse (Missing texvc executable; please see math/README to configure.): s(t)

• Failed to parse (Missing texvc executable; please see math/README to configure.): \,\!R_x

is the autocorrelation function of Failed to parse (Missing texvc executable; please see math/README to configure.): x(t)

• Failed to parse (Missing texvc executable; please see math/README to configure.): R_{x\,s}

is the cross-correlation function of Failed to parse (Missing texvc executable; please see math/README to configure.): x(t)
and Failed to parse (Missing texvc executable; please see math/README to configure.): s(t)

If the signal Failed to parse (Missing texvc executable; please see math/README to configure.): s(t)

and the noise Failed to parse (Missing texvc executable; please see math/README to configure.): n(t)
are uncorrelated (i.e., the cross-correlation Failed to parse (Missing texvc executable; please see math/README to configure.):  \,R_{sn}\, 

is zero) then note the following

• Failed to parse (Missing texvc executable; please see math/README to configure.): R_{x\,s} = R_s

• Failed to parse (Missing texvc executable; please see math/README to configure.): \,\!R_x = R_s + R_n

For most applications, the assumption of uncorrelated signal and noise is reasonable because the source of the noise (e.g. sensor noise or quantization noise) do not depend upon the signal itself.

The goal is to then minimize Failed to parse (Missing texvc executable; please see math/README to configure.): E(e^2)

by finding the optimal Failed to parse (Missing texvc executable; please see math/README to configure.): g(t)

.

Stationary solution

The Wiener filter has solutions for two possible cases: the case where a causal filter is desired, and the one where a non-causal filter is acceptable. The latter is simpler but is not suited for real-time applications. Wiener's main accomplishment was solving the case where the causality requirement is in effect.

Noncausal solution

Failed to parse (Missing texvc executable; please see math/README to configure.): G(s) = \frac{S_{x,s}(s)e^{\alpha s}}{S_x(s)}

Provided that Failed to parse (Missing texvc executable; please see math/README to configure.): g(t)

is optimal then the MMSE equation reduces to Failed to parse (Missing texvc executable; please see math/README to configure.): E(e^2) = R_s(0) - \int_{-\infty}^{\infty}{g(\tau)R_{x,s}(\tau + \alpha)\,d\tau}

And the solution, Failed to parse (Missing texvc executable; please see math/README to configure.): g(t)

is the inverse two-sided Laplace transform of Failed to parse (Missing texvc executable; please see math/README to configure.): G(s)

.

Causal solution

Failed to parse (Missing texvc executable; please see math/README to configure.): G(s) = \frac{H(s)}{S_x^{+}(s)}

Where

• Failed to parse (Missing texvc executable; please see math/README to configure.): H(s)

consists of the causal part of Failed to parse (Missing texvc executable; please see math/README to configure.): \frac{S_{x,s}(s)e^{\alpha s}}{S_x^{-}(s)}
(that is, that part of this fraction having a positive time solution under the inverse Laplace transform)

• Failed to parse (Missing texvc executable; please see math/README to configure.): S_x^{+}(s)

is the causal component of Failed to parse (Missing texvc executable; please see math/README to configure.): S_x(s)
(i.e. the inverse Laplace transform of Failed to parse (Missing texvc executable; please see math/README to configure.): S_x^{+}(s)
is non-zero only for Failed to parse (Missing texvc executable; please see math/README to configure.): t\ge 0

)

• Failed to parse (Missing texvc executable; please see math/README to configure.): S_x^{-}(s)

is the anti-causal component of Failed to parse (Missing texvc executable; please see math/README to configure.): S_x(s)
(i.e. the inverse Laplace transform of Failed to parse (Missing texvc executable; please see math/README to configure.): S_x^{-}(s)
is non-zero only for negative t)

This general formula is complicated and deserves a more detailed explanation. To write down the solution Failed to parse (Missing texvc executable; please see math/README to configure.): G(s)

in a specific case, one should follow these steps (see LLoyd R. Welch: Wiener Hopf Theory):

1. Start with the spectrum Failed to parse (Missing texvc executable; please see math/README to configure.): S_x(s)

in rational form and factor it into causal and anti-causal components:

Failed to parse (Missing texvc executable; please see math/README to configure.): S_x(s) = S_x^{+}(s) S_x^{-}(s)

where Failed to parse (Missing texvc executable; please see math/README to configure.): S^{+}

contains all the zeros and poles in the left hand plane (LHP) and Failed to parse (Missing texvc executable; please see math/README to configure.): S^{-}
contains the zeroes and poles in the RHP.

2. Divide Failed to parse (Missing texvc executable; please see math/README to configure.): S_{x,s}(s)e^{\alpha s}

by Failed to parse (Missing texvc executable; please see math/README to configure.): {S_x^{-}(s)}
and write out the result as a partial fraction expansion.

3. Select only those terms in this expansion having poles in the LHP. Call these terms Failed to parse (Missing texvc executable; please see math/README to configure.): H(s) .

4. Divide Failed to parse (Missing texvc executable; please see math/README to configure.): H(s)

by Failed to parse (Missing texvc executable; please see math/README to configure.): S_x^{+}(s)

. The result is the desired filter transfer function Failed to parse (Missing texvc executable; please see math/README to configure.): G(s)

FIR Wiener filter for discrete series

In order to derive the coefficients of the Wiener filter, we consider a signal w[n] being fed to a Wiener filter of order N and with coefficients Failed to parse (Missing texvc executable; please see math/README to configure.): \{a_i\} , Failed to parse (Missing texvc executable; please see math/README to configure.): i=0,\ldots, N . The output of the filter is denoted x[n] which is given by the expression

Failed to parse (Missing texvc executable; please see math/README to configure.): x[n] = \sum_{i=0}^N a_i w[n-i]

The residual error is denoted e[n] and is defined as e[n] = x[n] − s[n] (See the corresponding block diagram). The Wiener filter is designed so as to minimize the mean square error (MMSE criteria) which can be stated concisely as follows:

Failed to parse (Missing texvc executable; please see math/README to configure.): a_i = \arg \min ~E\{e^2[n]\}

where Failed to parse (Missing texvc executable; please see math/README to configure.): E\{\cdot\}

denote the expectation operator. In the general case, the coefficients Failed to parse (Missing texvc executable; please see math/README to configure.): a_i
may be complex and may be derived for the case where w[n] and s[n] are complex as well. For simplicity, we will only consider the case where all these quantities are real. The mean square error may be rewritten as:

Failed to parse (Missing texvc executable; please see math/README to configure.): \begin{array}{rcl} E\{e^2[n]\} &=& E\{(x[n]-s[n])^2\} \\ &=& E\{x^2[n]\} + E\{s^2[n]\} - 2E\{x[n]s[n]\}\\ &=& E\{\big( \sum_{i=0}^N a_i w[n-i] \big)^2\} + E\{s^2[n]\} -2 E\{ \sum_{i=0}^N a_i w[n-i]s[n]\} \end{array}

To find the vector Failed to parse (Missing texvc executable; please see math/README to configure.): [a_0,\ldots,a_N]

which minimizes the expression above, let us now calculate its derivative w.r.t Failed to parse (Missing texvc executable; please see math/README to configure.): a_i

Failed to parse (Missing texvc executable; please see math/README to configure.): \begin{array}{rcl} \frac{\partial}{\partial a_i} E\{e^2[n]\} &=& 2E\{ \big( \sum_{j=0}^N a_j w[n-j] \big) w[n-i] \} - 2E\{s[n]w[n-i]\} \quad i=0, \ldots ,N \\ &=& 2 \sum_{j=0}^N E\{w[n-j]w[n-i]\} a_j - 2 E\{ w[n-i]s[n] \} \end{array}

If we suppose that w[n] and s[n] are stationary, we can introduce the following sequences Failed to parse (Missing texvc executable; please see math/README to configure.): R_w[m] ~\textit{ and }~ R_{ws}[m]

known respectively as the autocorrelation of w[n] and the cross-correlation between w[n] and s[n] defined as follows

Failed to parse (Missing texvc executable; please see math/README to configure.): \begin{align} R_w[m] =& E\{w[n]w[n+m]\} \\ R_{ws}[m] =& E\{w[n]s[n+m]\} \end{align}

The derivative of the MSE may therefore be rewritten as (notice that Failed to parse (Missing texvc executable; please see math/README to configure.): R_{ws}[-i] = R_{sw}[i] )

Failed to parse (Missing texvc executable; please see math/README to configure.): \frac{\partial}{\partial a_i} E\{e^2[n]\} = 2 \sum_{j=0}^{N} R_w[j-i] a_j - 2 R_{sw}[i] \quad i=0, \ldots ,N

Letting the derivative be equal to zero, we obtain

Failed to parse (Missing texvc executable; please see math/README to configure.): \sum_{j=0}^N R_w[j-i] a_j = R_{sw}[i] \quad i=0, \ldots ,N

which can be rewritten in matrix form

Failed to parse (Missing texvc executable; please see math/README to configure.): \begin{bmatrix} R_w[0] & R_w[1] & \cdots & R_w[N] \\ R_w[1] & R_w[0] & \cdots & R_w[N-1] \\ \vdots & \vdots & \ddots & \vdots \\ R_w[N] & R_w[N-1] & \cdots & R_w[0] \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_N \end{bmatrix} = \begin{bmatrix} R_{sw}[0] \\R_{sw}[1] \\ \vdots \\ R_{sw}[N] \end{bmatrix}

These equations are known as the Wiener-Hopf equations. The matrix appearing in the equation is a symmetric Toeplitz matrix. These matrices are known to be positive definite and therefore non-singular yielding a single solution to the determination of the Wiener filter coefficients. Furthermore, there exists an efficient algorithm to solve the Wiener-Hopf equations known as the Levinson-Durbin algorithm.

The FIR Wiener filter is related to the Least mean squares filter, but minimizing its error criterion does not rely on cross-correlations or auto-correlations. Its solution converges to the Wiener filter solution.

See also

• Norbert Wiener
• Kalman filter
• Wiener deconvolution
• Eberhard Hopf
• Least mean squares filter
• Similarities between Wiener and LMS

References

• ^ [1]: Wiener, Norbert (1949), Extrapolation, Interpolation, and Smoothing of Stationary Time Series. New York: Wiley. ISBN 0-262-73005-7
• ^ [2]: Brown, Robert Grover and Patrick Y.C. Hwang (1996) Introduction to Random Signals and Applied Kalman Filtering. 3 ed. New York: John Wiley & Sons. ISBN 0-471-12839-2de:Wiener-Filter

it:Filtro di Wiener