EE264: Lecture 12

Wiener Filtering

In this lecture we will take a different view of filtering. Previously, we have depended on frequency-domain
specifications to make some sort of LP/ BP/ HP/ BS filter, which would extract the desired information from
an input signal.
Now, we wish to filter a signal x[n] to modify it such that it approximates some other signal d[n] in some
statistical sense. That is, the output of the filter y[n] is a good estimate of d[n]. The output error e[n]
represents the mismatch between y[n] and d[n]. This can be considered a time-domain specification of the
filter.
x[n]

y[n]

FILTER

e[n]

+
+
d[n]

Why would we want to do this? Well, consider the problem of predicting the stock market. Assume that
x[n] represents the price of a particular stock on day n and d[n] is the price of that stock one day in the
future. The goal is to find a filter that will predict d[n] given x[n]. Moreover, we would like to find the best
such linear predictor.
Q: Is there an obvious solution for this case?
A: Sure, pick H(z) = z, but this is worthless for this real-time example!
In order to talk about an optimal filter which estimates d[n] from x[n], we must have a method of measuring how good a job the filter does. A cost function is used to judge the performance, and could take on
many different forms.

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A.
131
Flores

132

Wiener Filtering

Lect. 12

For example,
=
E[u[|e[n]| k]]

Mean
Threshold
Error

error

Mean
Absolute
Error

= E[|e[n]|]

Mean
Square
Error

=
E[e2 [n]]

error

error

We most commonly use the mean square error (MSE) as our cost function.
= E[e2 [n]]
where E[] represents statistical expectation, and
2

E[e [n]] =

x2 pe (x) dx

where pe (x) is the probability density function of the error. The filter that is optimum in the MSE sense is
called a Wiener filter.
In all our analyses here we will assume:
1. x[n] is wide-sense stationary, i.e. it has constant (and finite) mean and variance, and a correlation
function which is a function of time shift only:
E[x[m1 ] x[m1 + n]] = E[x[m2 ] x[m2 + n]]

m 1 , m2

Incidentally, a stronger condition would be to assume that x[n] is strong-sense stationary, i.e. all of its
moments are constant and finite. An even stronger condition would be to assume that x[n] is ergodic,
i.e. all of its sample averages converge to the true moments. Fortunately, these latter two conditions
are not necessary for the present development.
2. All of the signals are zero-mean.
3. We use MSE as our error criterion.

Correlation
The definition of correlation between two random processes x[n] and y[n] (when x and y are wide sense
stationary) is:
xy [k] = E[x[n]y[n + k]]
c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

Lect. 12

Wiener Filtering

133

It is a measure of the similarity between two random variables, and is related to the amount of information
that one random variable contains about another.
Example:
xx

xx

Highly
Correlated

Uncorrelated
(white)

Two-Sided (Unconstrained) Wiener Filters

=
=
=
=

E[e2 [n]]
E[(d[n] y[n])2 ]
E[d2 [n]] 2E[y[n]d[n]] + E[y 2 [n]]
dd [0] 2yd [0] + yy [0]

So, the cost function is itself a function of the correlation and cross correlation of the output of the filter and
the desired output. We know that:
y[n] =

m=

h[m]x[n m]

Observe that, in general, h[m] represents a non-causal filter with infinite impulse response (IIR). (Later, we
consider the special case where h[m] is either a causal IIR filter or a causal FIR filter.) Substituting the
definition for y[n] into the previous result, we get:
yd [0] = E[y[n]d[n]]
= E
=

"

h[m]x[n m]d[n]

m=

h[m]xd [m]

m=

yy [0] = E[y[n]y[n]]

= E
=

h[m]x[n m]

m=

X X

m= l=

= dd [0] 2

"

l=

h[l]x[n l]

h[m]h[l]xx [l m]

m=

h[m]xd [m] +

m= l=

h[m]h[l]xx [l m]

Our goal is to find h[n] which minimizes . The solution is to take the derivative

d
d h[ni ]

and set it to zero:

X
d
= 0 2xd [ni ] + 2
hopt [m]xx [ni m] = 0
d h[ni ]
m=

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

134

Wiener Filtering

m=

Lect. 12

hopt [m]xx [ni m] = xd [ni ]

which must be satisfied for all ni .

This equation specifies the impulse response of the optimal filter. Observe that it depends on the cross
correlation between x[n] and d[n] and the autocorrelation function of x[n].
The optimal filter is known as the Wiener solution. Observing that the left hand side of the equation is equal
to the convolution of hopt [n] and xx [n],
Hopt (z)xx (z) = xd (z)
Hopt (z) = xd (z)1
xx (z)
Note:
The filter obtained may be non-causal!
min = dd [0]

n=

hopt [n]xd [n] 6= 0

in general.

Example: Noise Canceling:

x[n]= s[n] + v[n]

y[n]

Noise
Canceler
e[n]

+
+
d[n] = s[n]

Assuming that s[n] and v[n] are uncorrelated, zero mean:

E[x[n]x[n + k]] =
=
xx (z) =
E[x[n]d[n + k]] =
=
xd (z) =
Thus,

E[(s[n] + v[n])(s[n + k] + v[n + k])]

ss [k] + vv [k]
ss (z) + vv (z)
E[(s[n] + v[k])s[n + k]]
ss [k]
ss (z)
ss (z)
Hopt (z) =
ss (z) + vv (z)

We notice two properties from this solution:

1. At frequencies where vv (z) 0 (i.e. the noise is zero) there is no need to filter the signal and
Hopt (z) 1.
2. At frequencies where vv (z) (i.e. where the noise dominates) the filter is turned off and
Hopt (z) 0.
c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

Lect. 12

Wiener Filtering

135

Numerical Example: Assume that the autocorrelation function of the input signal is:
 |k|

1
2

10
ss [k] =
27

and that the autocorrelation of the (white) noise source is:

2
vv [k] = [k]
3
The signal and noise are assumed to be uncorrelated with each other.
ss

vv

To find the Wiener filter, we need to obtain xd (z) and xx (z). We can first find xd [k] and xx [k], and
then z-transform them. Since the signal and noise are uncorrelated with each other,
xx [k] = ss [k] + vv [k]
 
10 1 |k| 2
=
+ [k],
27 2
3
xx (z) = ss (z) + vv (z).
The z-transform of the noise autocorrelation is
vv (z) =

2
.
3

The z-transform of the signal autocorrelation (after some math) is

ss (z) =

5/18
(1

1 1
2 z )(1

12 z)

Combining these:
xx (z) =

5/18
1 1
2 z )(1

12 z)
6z 1

(1
20 6z
=
.
18(1 12 z 1 )(1 12 z)

2
3

We saw earlier,
xd [k] = ss [k],
xd (z) = ss (z).
Combining these results:
Hopt (z) =

5/18
(1 13 z 1 )(1 31 z)

5
hopt [n] =
16

 |n|

1
3

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

136

Wiener Filtering

Lect. 12

Q: How did we pick this hopt [n] without knowing the ROC of Hopt (z)?
A: We need a stable noise canceller.
Clearly this resulting filter is non-causal. We have two recourses:
1. Limit the filter to be a causal IIR filter.
2. Compute the optimal FIR filter instead.

Causal (Shannon-Bode) Wiener Filters

Our interest now focuses on the realization of causal Wiener filters, whose impulse responses are constrained
to be zero for negative time. The optimal causal impulse response has zero response for negative time and
has zero derivatives of with respect to impulse response for all times equal to and greater than zero. The
causal Wiener equation becomes:

m=0

hopt causal [m]xx [n m] = xd [n],

hopt causal [n] = 0,

n0

n < 0.

This is not a simple convolution like the unconstrained Wiener equation, and special methods will be needed
to find useful solutions. The approach developed by Shannon and Bode will be used. We begin with a simple
case. Let the filter input be white with zero-mean and unit variance, so that
xx [n] = [n].
For this input, the causal Wiener equations become

m=0

hopt causal [m]xx [n m] = hopt causal [n] = xd [n],

hopt causal [n] = 0,

n0

n<0

With the same white input, but without the causality constraint, the Wiener equation would be

m=

hopt [m]xx [n m] = hopt [n] = xd [n],

n.

From this we may conclude that when the input to the Wiener filter is white, the optimal solution with a
causality constraint is the same as the optimal solution without the constraint, except that with the causality
constraint the impulse response is set to zero for negative time. With a white input, the causal solution is
easy to obtain, and it is key to Shannon-Bode. You find the unconstrained two-sided Wiener solution and
remove the non-causal part in the time domain. You dont do this if the input is not white.
Usually, the input to the Wiener filter is not white. Therefore the first step in the Shannon-Bode realization
is to whiten the input signal. A whitening filter can always be designed to do this, using a priori knowledge
of the input autocorrelation function or its z-transform.
Assume that the z-transform of the autocorrelation function is a rational function of z that can be written as
the ratio of a numerator polynomial in z to a denominator polynomial in z. Factor both the numerator and
denominator polynomials. The autocorrelation function is symmetric, and its z-transform has the following
symmetry:
xx [n] = xx [n],
xx (z) = xx (z 1 ).
c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

Lect. 12

Wiener Filtering

137

Accordingly, there must be symmetry in the numerator and denominator factors:

xx (z) = A

(1 az 1 )(1 az)(1 bz 1 )(1 bz) . . .

.
(1 z 1 )(1 z)(1 z 1 )(1 z) . . .

With no loss in generality,1 assume that all the parameters a, b, c, . . . , , , . . . have magnitudes less
than one. We can then factor xx (z) as

xx (z) = +
xx (z) xx (z),

where
(1 az 1 )(1 bz 1 ) . . .
A
(1 z 1 )(1 z 1 ) . . .
(1 az)(1 bz) . . .

A
.
xx (z) =
(1 z)(1 z) . . .

+
xx (z) =

All poles and zeros of +

xx (z) will be inside the unit circle in the z-plane. All poles and zeros of xx (z)
will be outside the unit circle in the z-plane. Furthermore,

1
+
1

+
xx (z) = xx (z ), xx (z ) = xx (z).

The whitening filter can now be designed. Let it have the transfer function
Hwhitening (z) =

.
+
xx (z)

To verify its whitening properties, let its input be x[n] with autocorrelation function xx [n], and let its output
be x
[n] with autocorrelation function xx [n]. The transform of the output autocorrelation function is
xx (z) = Hwhitening (z 1 ) Hwhitening (z) xx (z)
1
1
xx (z)
= + 1 +
xx (z ) xx (z)
1
1
=
+
xx (z)
xx (z) xx (z)
= 1.
The autocorrelation function of the output is
xx [n] = [n].
Therefore the output is white with unit mean square.
It is important to note that the whitening filter is both causal and stable, since the zeros of +
xx (z) are all
inside the unit circle in the z-plane. Furthermore, using this whitening filter does nothing irreversible to the
signal. The inverse of the whitening filter has a transfer function equal to +
xx (z), which is also stable and
causal because +
(z)
has
all
of
its
poles
inside
the
unit
circle.
The
whitening
filter can be readily utilized
xx
since it is causal, stable, and invertible with a causal, stable filter.
1

The only case not included involves zeros of xx (z) exactly on the unit circle. This is a special case and it requires special
treatment. These zeros are assumed to be somewhere off the unit circle, and they are placed back on it by a limiting process.

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

138

Wiener Filtering

x
[n]

x[n]

+
xx (z)

Yopt (z)

Lect. 12

y[n]
+

e[n]

+
d[n]

The subsequent filter is easily designed, since its input is white. We represent its transfer function by
Yopt causal (z) and can design it by first obtaining Y opt (z) disregarding causality, and then removing the
non-causal part of its impulse response.
We wish to find Yopt (z) disregarding causality in terms of the autocorrelation of x[n] and crosscorelation
with d[n]. To do so, we could work out expressions for xx (z) and xd (z) in terms of xx (z) and xd (z)
and then use the formula for the uncostrained Wiener filter. However, there is an easier way; notice that the
xd (z)
cascading of +1(z) and Yopt (z) should be equal to the unconstrained Wiener filter given by
, hence
xx (z)
xx

xd (z)
+
xx (z)
xx (z)
xd (z)
=
xx (z)

Yopt (z) =

To obtain Yopt causal (z), we first inverse transform Yopt (z) into the time domain, remove the non-causal
part, and then z-transform the causal part obtained. This operation cannot be done with z-transforms alone.
A special notation has been devised to represent taking the causal part:
"

xd (z)
Yopt causal (z) = [Yopt (z)]+ =

xx (z)

.
+

The Shannon-Bode realization of the causal Wiener filter can now be formulated:
"

xd (z)
1
Hopt causal (z) = +
xx (z)
xx (z)

.
+

Example:
An example helps to illustrate how this formula is used. We will rework the above noise filtering example,
only in this case the Wiener filter will be designed to be causal.
The first step is to factor xx :
xx (z) =

(1 13 z)(1 31 z 1 )
.
(1 12 z)(1 21 z 1 )

Therefore,
+
xx =

(1 13 z 1 )
(1 31 z)

=
.
xx
(1 12 z 1 )
(1 12 z)

xd (z) was found before. Thus:

xd (z) =

5/18
(1

1 1
2 z )(1

21 z)

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

Lect. 12

Wiener Filtering

139

xd (z)
5/18
=

1 1
xx (z)
(1 2 z )(1 31 z)

Yopt (z) =

"

1 1
6z
12 z 1

1
3

1 13 z

The filter Yopt (z) is generally two-sided in the time domain. In order for it to be stable, its two terms must
correspond to time domain components as follows:
1

1 1
6z
12 z 1
1
3
1 13 z

1
6

(right-handed time function).

(left-handed time function).

1
1
12

1
24

1
3

1
3

1
27

1 13 z

1
9

"

1
3

1
27

1
6

1
9

1 1
z
6
12 z 1

1
12

"
1
12

1 1
z
6
12 z 1

1
3

1 13 z

= Yopt (z)

1
24

1
3
1
6

1 1
z
6
12 z 1

1
3

1 13 z

= Yopt causal (z)

+

1
24

Yopt causal (z) =

=

1
1 1 1
+ z + z 2 + . . .
3 6
!12
1

1
3
1 1
2z

Including the whitening filter, the transfer function of the causal Wiener filter is
Hopt causal (z) =
=

1 12 z 1
1 13 z 1
1
3
1 1 .
z
3
n

1

1 1
hopt causal [n] =
3 3

1
3
1 1
2z

n 0.

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

140

Wiener Filtering

1
3
1
9

1
1
27

1 1
z
3
13 z 1

Lect. 12

= Hopt causal (z)

1
81

It is expected that the causal Wiener filter would not perform as well as the unconstrained non-causal Wiener
filter. This is indeed the case, and can be verified with a little math.
min

causal

= dd [0]

hopt causal [n]xd [n].

n=0

Optimal FIR Wiener Filters

For an FIR realization, we make the following changes in the derivation:
y[n] =

N
1
X
k=0

h[k]x[n k]

This time, we will solve the system using matrix notation:

Filter
H = [h0 , h1 , . . . , hN 1 ]T
Input
X = [xn , xn1 , . . . , xnN +1 ]T
Output y[n] = H T X = X T H
(scalar)
Error
e[n] = d[n] y[n] = d[n] H T X
= E[e2 [n]]
= E[d2 [n]] 2E[H T Xd[n]] + E[H T XX T H]
= dd [0] 2H T E[Xd[n]] +H T E[XX T ] H
| {z }

d2

{z
P

{z
R

Note: is a quadratic function of the filter H. This implies (among other things) that there is a unique global
optimum.
P = E[Xd[n]]

E[x[n]d[n]]
..
= .

E[x[n N + 1]d[n]

p[0]

..
= .

p[N 1]

If X and d[n] are uncorrelated, then P = 0.

R = E[X(X)T ]

E[x[n]x[n]]
..
= .

..
.
E[x[n]x[n N + 1]

E[x[n N + 1]x[n]]
..
.
E[x[n N + 1]x[n N + 1]]

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores

Lect. 12

Wiener Filtering

141

r[0]

..
..
= .
.
r[N 1]

r[N 1]

..

r[0]

If x[n] is a white noise process with unit variance, then R = I, the identity matrix.
Solution: The minimum error location will be where:

T

where
(H) =
,...,
(Hopt ) = 0
h[0]
h[N 1]

As a result,

2P + 2RHopt = 0
RHopt = P
A correlation matrix R is always positive semi-definite, so unless R = 0, its inverse will exist:
Hopt = R1 P
The minimum error can be found too:
min =
=
=
=

(Hopt )
d2 2P T RT P + P T RT RR1 P
d2 P T R1 P
d2 P T Hopt

R T =
R1 since
R is a cov.
matrix

and since R is positive semi-definite,

min d2

Without filtering, the error is d2 , so we do better!

c
c
c
Copyright 19951996
by G. Plett. Copyright 19982004
by M. Kamenetsky. Copyright 2005
by A. Flores