Time domain Wiener filter - Unknown signal in white Gaussian noise

1) IntroductionPermalink

The objective for the Wiener filter is to retrieve a desired signal $s(n)$ in noise $v(n)$ given the observed signal $x(n)$ i.e.:

$x(n) = s(n) + v(n).$

The time domain Wiener filter is a linear minimum mean square error estimator (LMMSEE) with optimal filter coefficients $\textbf{w}^\ast$ that satisfy following optimization problem:

$\textbf{w}^\ast = \underset{\textbf{w}}{\text{arg min}} \quad \mathbb{E} \left[ \left( s(n) - \textbf{w}^T \textbf{x}(n) \right)^2 \right]$,

where $\textbf{w} = [w_0,…,w_{M-1}]^T \in \mathbb{R}^{M}$ and $\textbf{x}(n) = [x(n),x(n-1)…,x(n-M+1)]^T \in \mathbb{R}^{M}$ such that the LMMSE estimate of $s(n)$ given $\textbf{x}(n)$ is:

$\hat{s}(n) = \textbf{w}^{\ast T}\textbf{x}(n)$.

2) AssumptionsPermalink

3) SolutionPermalink

The Wiener filter is given as

$\textbf{w}^\ast = \hat{\textbf{R}}_{xx}^{-1} \hat{\textbf{r}}_{xs}$,

where $\hat{\textbf{R}}_{xx} \in \mathbb{R}^{M\times M}$ is the sample estimate of the autocorrelation matrix of $x(n)$ and $\hat{\textbf{r}}_{xs} \in \mathbb{R}^M$ is the sample estimate of the cross-correlation vector between $s(n)$ and $x(n)$.

4) Step-by-step guidePermalink

1. Set filter length e.g. $M=20$.

2. Set autocorrelation of the noise to

   $r_{vv}(k) = \sigma_v^2 \delta (k), \quad \text{for }k = 0,…,M-1$

   where $\delta (k)$ is the Kronecker delta function and $k$ is the lag-index.

3. Compute the unbiased sample estimate of the autocorrelation of $x(n)$

   $$\hat{r}_{xx}(k) = \frac{1}{N-k} \sum\limits_{n=0}^{N-k-1} x(n)x(n+k), \quad \text{for }k = 0,...,M-1$$

   or equivalently

   $$\hat{r}_{xx}(k) = \frac{1}{N-k} \textbf{x}^T(n)\textbf{x}(n+k), \quad \text{for }k = 0,...,M-1$$

   where $\textbf{x}(n) = [\textbf{x}(0),...,\textbf{x}(N-k-1)]^T$ and $\textbf{x}(n+k) = [\textbf{x}(k),...,\textbf{x}(N-1)]^T$.

4. Estimate the autocorrelation of $s(n)$ as

   $$\hat{r}_{ss}(k) = \hat{r}_{xx}(k) - r_{vv}(k)$$
5. Form the autocorrelation matrix $\hat{\textbf{R}}_{xx}$ of $x(n)$

6. Form the cross-correlation vector between $x(n)$ and $s(n)$

   $$\hat{r}_{xs}(k) = \hat{r}_{ss}(k)$$

7. Compute the Wiener filter coefficients

   $$\text{w}^{\ast} = \hat{\textbf{R}}_{xx}^{-1} \hat{\textbf{r}}_{xs}$$

8. Perform filtering

   $\hat{s}(n) = \text{w}^{\ast T} \textbf{x}(n)$

5) MATLAB codePermalink

A MATLAB example is provided where the signals are artificially generated.

%   Copyright 2020: DSPCookbook
clc, clear, close all

% Constants in simuation
N       = 20000;
varu    = 1;
varv    = 100;
a       = 0.999;

% Generate the noisy observed signal
u       = sqrt(varu)*randn(N,1);
s       = 0; % Initialize s to be the mean
for n = 2:N
    s(n,1) = a*s(n-1,1) + u(n);
end
v       = sqrt(varv)*randn(N,1);
x       = s + v;

% Step 1
M       = 50;
tau     = (0:(M-1))';

% Step 2
rvv     = zeros(M,1);
rvv(1)  = varv;            

% Step 3
kk = 1;
for k = 0:M-1
    rxx(kk,1) = 1/(N-k)*x(1:(N-k))'*x(kk:N);
    kk = kk + 1;
end

% Step 4
rss = rxx - rvv;

% Step 5
Rxx = toeplitz(rxx);

% Step 6
rxs     = rss;

% Step 7
w       = inv(Rxx)*rxs;

% Step 8
for n = M:N
    y(n,1) = w'*flipud(x((n-M+1):n));
end

%% Plot result
plot(x)
grid, xlabel('time index n'), ylabel('Amplitude'), hold on
plot(y)
plot(s)
legend('x(n)','y(n)','s(n)')
title('Wiener filtering')

Derivation of the Wiener filterPermalink

For the derivation of the time domain Wiener filter, check out the extra material.