The Wiener Filter 3.1 The Wiener-Hopf Equation
The Wiener Filter 3.1 The Wiener-Hopf Equation
The Wiener Filter 3.1 The Wiener-Hopf Equation
1. The assumption that both signal and noise are random processed with known spectral
characteristics or, equivalently, known auto- and cross-correlation functions.
2. The criterion for best performance is minimum mean-square error. (This is partially to make
the problem mathematically tractable, but it is also a good physical criterion in many
applications.)
3. A solution based on scalar methods that leads to the optimal filter weighting function (or
transfer function in the stationary case).
{
J = E [ s (n) − sˆ(n) ]
2
}
Fig. 3.1-1 Wiener Filter Problem
1
We now consider the filter optimization problem that Wiener first solved in the 1940s. Referring
to Fig. 3.1-1, we assume the following:
1. The filter input is an additive combination of signal and noise, both of which are jointly wide-
sense stationary (JWSS) with known auto- and cross-correlation functions (or corresponding
spectral functions).
2. The filter is linear and time-invariant. No further assumption is made as to its form.
3. The output is also wide-sense stationary.
4. The performance criterion is minimum mean-square error.
∞ ∞
ŝ(n) = h(n) ∗ z(n) = ∑ h(n − i)z(i) = ∑ h(i)z(n − i),
i =−∞ i =−∞
(3.1-1)
where z(i) is the measurement and h(n) is the impulse response of the estimator. Let H denote
the region of support of h(n) , defined by
H = {n : h(n) ≠ 0}.
2
ŝ(n) = ∑ h(i )z(n − i ) . (3.1-2)
i∈H
Let the mean-square estimation error (MSE) J be
{
J = E [s(n) − s(
ˆ n) ]
2
}
⎧⎪ ⎡ ⎤⎡ ⎤ ⎫⎪
= E ⎨ ⎢s(n) − ∑ h(i )z(n − i ) ⎥ ⎢s(n) − ∑ h( j )z(n − j ) ⎥ ⎬ (3.1-3)
⎪⎩ ⎣ i∈H ⎦⎣ j∈H ⎦ ⎪⎭
= E {s 2 (n)} − 2∑ h(i ) E {s(n)z(n − i )} + ∑ ∑ h(i )h( j ) E {z(n − i )z(n − j )}.
i∈H i∈H j∈H
To minimize the MSE, take the partial derivatives of J with respect to h(i ) , for each h(i ) ≠ 0.
Then, set the result equal to zero
∂J
= −2 E {s(n)z(n − i )} + 2 ∑ h( j ) E {z(n − i )z(n − j )}
∂h(i ) j∈H (3.1-4)
= 0.
3
which we may express in the form of
∑ h( j ) R ( j − i ) = R
j∈H
z sz (i ), i ∈ H , (3.1-6)
where Rz (k ) is the autocorrelation function of z(n) and Rsz (k ) is the crosscorrelation function
of s(n) and z(n).
Eq. (3.1-6) is the discrete-time Wiener-Hopf equation. It is the basis for the derivation of the
Wiener filter.
Suppose the filter has an impulse response h(n) with support H, where H is a finite set, e.g., H =
{0,1,2,…,N-1}. Then the impulse response h(n) has a finite duration, and this type of filter is
called a finite impulse response (FIR) filter.
∑ h( j ) R ( j − i ) = R
j =0
z sz (i ), 0 ≤ i ≤ N − 1 . (3.2-1)
4
This may be written in matrix form
where we have used the fact that Rz (k ) = Rz (−k ) . Denote Eq. (3.2-2) in a convenient form
Rz h = r sz (3.2-3)
h = Rz −1 r sz . (3.2-4)
This is the FIR Wiener filter of order N − 1. Note that Rz is a Toeplitz matrix and Eq. (3.2-4)
may be solved using a computationally-efficient method such as the Levinson-Durbin algorithm.
5
Mean-Square Error (MSE)
⎧ ⎡ N −1
⎤⎫
n) ⎢s(n) − ∑ h(i )z(n − i ) ⎥ ⎬
MSE = E ⎨s(
⎩ ⎣ i =0 ⎦⎭
(3.2-5)
N −1
n)s(n)} − ∑ h(i ) E {s(
= E {s( n)z(n − i )}.
i =0
⎧⎪ ⎡ ⎤ ⎫⎪
E ⎨z(n − i ) ⎢s(n) − ∑ h( j )z(n − j ) ⎥ ⎬ = 0
⎪⎩ ⎣ j∈H ⎦ ⎪⎭
E {z(n − i ) [s(n) − s(
ˆ n) ]} = 0 (3.2-6)
E {z(n − i )s(
n)} = 0, i ∈ H.
This relation is known as the orthogonality principle for LTI estimators. Applying this relation to
Eq. (3.2-5) gives
6
MSE = E {s(
n)s(n)}
N −1
= E {s (n)} − ∑ h(i ) E {s(n)z(n − i )}
2
i =0
N −1
= Rs (0) − ∑ h(i ) Rsz (i )
i =0 (3.2-7)
= Rs (0) − hT r sz .
Observe that if no filtering is performed and we simply use ŝ(n) = z (n) , the MSE is
{
MSE no filter = E [s(n) − z (n) ]
2
} = R (0) − 2R (0) + R (0).
s sz z (3.2-8)
We can use Eq. (3.2-8) to measure the effectiveness of the Wiener filter. The reduction in MSE
due to Wiener filtering is given by
⎛ MSE no filter ⎞
reduction in MSE = 10log10 ⎜ ⎟ dB . (3.2-9)
⎝ MSE filter ⎠
We can inspect the MSE to determine whether the filter performance is acceptable. If the MSE is
too high, we may wish to use a larger value of N.
7
Example 3.2-1
Suppose we have a signal s(n) with autocorrelation function Rs (k ) = 0.95 k . The signal is
observed in the presence of additive white noise with variance σ v 2 = 2 . Hence Rv (k ) = 2δ (k ) .
s(n) and v(n) are uncorrelated, zero-mean, JWSS random processes. We would like to design
the optimum LTI second-order filter h.
Since the filter is the second-order, we set N = 3 . The matrix equation Eq. (3.2-2) to be solved is
8
Rz (k ) = E {z(n)z(n − k )}
= E {[s(n) + v(n) ][s(n − k ) + v(n − k ) ]}
(3.2-11)
= E {[s(n)s(n − k ) ]} + E {[ v(n)v( n − k ) ]}
= Rs (k ) + Rv (k ),
and
Rz (k ) = 0.95 + 2δ (k )
k
Rsz (k ) = 0.95 .
k
9
⎡3 0.95 0.9025⎤ ⎡ h(0) ⎤ ⎡1 ⎤
⎢0.95 3 0.95 ⎥ ⎢ ⎥ ⎢
h(1) = 0.95 ⎥ . (3.2-13)
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣⎢0.9025 0.95 3 ⎦⎥ ⎣⎢ h(2) ⎦⎥ ⎣⎢0.9025⎦⎥
⎛ MSE no filter ⎞
reduction in MSE = 10log10 ⎜ ⎟ dB
⎝ MSE filter ⎠
⎛ 2 ⎞
= 10log10 ⎜ ⎟ (3.2-15)
⎝ 0.4405 ⎠
≈ 6.5dB.
10
3.3 The Noncausal Wiener Filter
∑ h( j ) R ( j − i ) = R
j∈H
z sz (i ), i ∈ H
and let H = Z , where Z is the set of integers: Z = {", −2, −1,0,1, 2,"}. Then, the Wiener-Hopf
equation becomes
∞
∑ h( j ) R ( j − i ) = R
j =−∞
z sz (i ), for all i ∈ Z . (3.3-1)
H ( z ) Sz ( z ) = Ssz ( z )
so that
Ssz ( z )
H ( z) = . (3.3-2)
Sz ( z )
Eq. (3.3-2) is the solution for the noncausal Wiener filter.
11
Mean-Square Error (MSE)
As with the FIR Wiener filter of Section 3.2, we can derive an expression for the MSE. The MSE
similar to Eq. (3.2-7) can be written as
∞
MSE = Rs (0) − ∑ h(i ) Rsz (i ) . (3.3-3)
i =−∞
Example 3.3-1
Let us apply a noncausal Wiener filter to the filtering problem of Example 3.2-1. In that example
we found that
1 − (0.95) 2
Sz ( z ) = +2
(1 − 0.95 z −1 )(1 − 0.95 z )
2.3955(1 − 0.7931z −1 )(1 − 0.7931z ) 1
= , 0.7931 < z <
(1 − 0.95 z −1 )(1 − 0.95 z ) 0.7931
12
and
1 − (0.95) 2 1
Ssz ( z ) = , 0.95 < z < .
(1 − 0.95 z −1 )(1 − 0.95 z ) 0.95
Ssz ( z ) 0.0975
H ( z) = =
Sz ( z ) 2.3955(1 − 0.7931z −1 )(1 − 0.7931z )
0.1097(1 − (0.7931) 2 ) 1
= , 0.95 < z < .
(1 − 0.7931z −1 )(1 − 0.7931z ) 0.95
h(n) = 0.1097(0.7931) .
n
The mean-square errors associated with the noncausal Wiener filter is according to Eq. (3.3-3),
13
∞
MSE = Rs (0) − ∑ h( n) R
n =−∞
sz ( n)
∞
= (0.95) −
0
∑ 0.1097(0.7931)
n =−∞
n
(0.95)
n
⎡∞ −∞
⎤
= 1 − 0.1097 ⎢ ∑ (0.7931) (0.95) + ∑ (0.7931) − n (0.95) − n ⎥
n n
⎣ n =0 n =−1 ⎦
⎡ ∞ ⎤
= 1 − 0.1097 ⎢ 2∑ (0.7931) n (0.95) n − 1⎥
⎣ n =0 ⎦
⎡ 2 ⎤
= 1 − 0.1097 ⎢ − 1⎥ = 0.2195.
⎣1 − (0.7931)(0.95) ⎦
Since the MSE with no filter is 2 according to Eq. (3.2-14), the improvement of the noncausal
filter over no filtering is
⎛ MSE no filter ⎞
reduction in MSE = 10log10 ⎜ ⎟ dB
⎝ MSE filter ⎠
⎛ 2 ⎞
= 10log10 ⎜ ⎟
⎝ 0.2195 ⎠
≈ 9.6dB.
14
Compared to the second-order filter in Eq. (3.2-15) in Ex. (3.2-1), the noncausal Wiener filter
reduces the estimation error by an additional 3 dB .
∑ h( j ) R ( j − i ) = R
j∈H
z sz (i ), i ∈ H
and let H = N , where N is the set of integers: N = {0,1,2,"}. Then, the Wiener-Hopf equation
becomes
∞
∑ h( j ) R ( j − i ) = R
j =0
z sz (i ), for all i ∈ N . (3.4-1)
Eq. (3.4-1) can not be simplified by taking the z-transform because of its causality restriction. To
solve the Wiener-Hopf equation, the spectral factorization and the causal-part extraction are
necessary.
Let x(n) be a real-valued, zero-mean, WSS random process with power density spectrum S x ( z ) ,
15
where S x ( z ) is rational in z and has no poles on the unit circle. Then S x ( z ) can be factored into
the product
Sx ( z) = Sx + ( z)Sx − ( z ) , (3.4-2)
where
Example 3.4-1
where w(n) is a zero-mean, WSS random process with autocorrelation function Rw (n) = 5(0.6) n .
16
We want to find the spectral factorization of Ss ( z ).
S ( z) 2 + 3z −1 2(1 + 1.5 z −1 )
H ( z) = = = .
W ( z ) 1 − 1.1z −1 + 0.24 z −2 (1 − 0.3z −1 )(1 − 0.8 z −1 )
Using the input-output power spectrum relation, the power spectrum of s(n) is given by
Ss ( z ) = H ( z ) H ( z −1 ) S w ( z )
2(1 + 1.5 z −1 ) 2(1 + 1.5 z ) 3.2 (3.4-4)
= −1 −1 −1
,
(1 − 0.3 z )(1 − 0.8 z ) (1 − 0.3z )(1 − 0.8 z ) (1 − 0.6 z )(1 − 0.6 z )
1 − (0.6) 2
Sw ( z) = 5 .
(1 − 0.6 z −1 )(1 − 0.6 z )
17
Collect the terms in Eq. (3.4-4) that have poles or zeros inside the unit circle to form Ss + ( z ) ,
2 3.2(1 + 1.5 z )
Ss + ( z ) = . (3.4-5)
(1 − 0.3z −1 )(1 − 0.8 z −1 )(1 − 0.6 z −1 )
Collect the terms in Eq. (3.4-4) that have poles or zeros outside the unit circle to form Ss − ( z ) ,
− 2 3.2(1 + 1.5 z −1 )
Ss ( z ) = . (3.4-6)
(1 − 0.3z )(1 − 0.8 z )(1 − 0.6 z )
Ss − ( z ) = Ss + ( z −1 ) .
Consider the impulse response h(n) of a real-valued LTI system with rational z-transform H ( z ) .
In the time domain we can split h(n) into its causal and anticausal parts such that
where
18
causal part of h(n) ≡ [ h(n) ]+ = h(n)1(n),
anticausal part of h(n) ≡ [ h(n) ]− = h(n)1(−n − 1).
H ( z ) = [ H ( z ) ]+ + [ H ( z ) ]− , (3.4-8)
where
[ H ( z )]+ = Z {h(n)1(n)} ,
(3.4-9)
[ H ( z )]− = Z {h(n)1(−n − 1)}.
We now develop a method for determining [ H ( z ) ]+ and [ H ( z ) ]− . By long division in z, z −1 , or
both for a rational H ( z ) , we can always convert it into the form
⎛ N −1 − n ⎞
L M −N ⎜ ∑ bn z ⎟
H ( z ) = ∑ c− n z n + ∑ cn z − n + ⎜ nN=0 ⎟ (3.4-10)
n =1
n =0
⎜⎜ ∑ an z − n ⎟⎟
PA ( z ) PC ( z )
⎝
n =0
⎠
Q( z)
19
H ( z ) consists of two polynomial terms, PA ( z ) and PC ( z ) , and a proper fraction in z −1 , Q ( z ) .
PA ( z ) corresponds to a purely anticausal sequence, and PC ( z ) to a purely causal sequence.
Now we examine Q( z ) . Let K ≤ N be the number of distinct poles of Q( z ) , and denote these
poles and their associated degrees by pk and mk , respectively, k = 1,", K . Then Q( z ) may be
written as
N −1
⎛1⎞ ∑b z n
−n
Q( z ) = ⎜ ⎟ K
n =0
.
⎝ a0 ⎠
∏ (1 − p z
k =1
k
−1 mk
)
Suppose that H ( z ) be stable. Then its ROC includes the unit circle and the ROC of each Qk ( z )
must include the unit circle. Let N A be the total order of the poles of H ( z ) outside the unit
circle, and NC be the total order of the poles of H ( z ) inside the unit circle with the exception
20
of the origin. Then N A + N C = N and we may write
N A −1
K ∑λ z n
−n
QA ( z ) = ∑ Q ( z) = k K
n =0
,
k =1
pk >1 ∑ (1 − p z
k =1
k
−1 mk
)
pk >1
and
NC −1
K ∑μ z n
−n
QC ( z ) = ∑ Qk ( z ) = K
n =0
.
k =1
0< pk <1 ∑
k =1
(1 − pk z −1 ) mk
0< pk <1
[ H ( z )]+ = PC ( z ) + QC ( z ) ,
which is the sum of a purely causal sequence and all terms in the partial fraction expansion of
21
H ( z ) that have poles inside the unit circle.
Example 3.4-2
⎧ 1 ⎫
with ROC = ⎨ z : < z < 4 ⎬ .
⎩ 2 ⎭
22
3z 2 +4
2 − 17 z −1 + 40 z −2 − 16 z −3 6 z 2 − 51z + 128 − 109 z −1 + 197 z −2 − 232 z −3 + 80 z −4
6 z 2 − 51z + 120 − 48z −1
________________________________________
8 − 61z −1 + 197 z −2 − 232 z −3 + 80 z −4
8 − 68z −1 + 160 z −2 − 64 z −3
_____________________________
7 z −1 + 37 z −2 − 168 z −3 + 80 z −4
7 z −1 + 37 z −2 − 168 z −3 + 80 z −4
H ( z ) = 3z + 4 +
2
.
2 − 17 z −1 + 40 z −2 − 16 z −3
23
−5 z −1 − 2
−16 z −3 + 40 z −2 − 17 z −1 + 2 80 z −4 − 168 z −3 + 37 z −2 + 7 z −1
80 z −4 − 200 z −3 + 85 z −2 − 10 z −1
________________________
32z −3 − 48 z −2 + 17 z −1
32z −3 − 80 z −2 + 34 z −1 − 4
____________________
32z −2 − 17 z −1 + 4
−1 4 − 17 z −1 + 32 z −2
H ( z ) = 3z + 4 + (−2) − 5 z +
2
.
2 − 17 z −1 + 40 z −2 − 16 z −3
17 −1
−1 −2 2− z + 16 z −2
4 − 17 z + 32 z 2
Q( z ) = −1 −2 −3
=
2 − 17 z + 40 z − 16 z 1
(1 − z −1 )(1 − 4 z −1 ) 2
2
1 1
= + .
(1 − 4 z −1 ) 2 1 − 1 z −1
2
24
Hence,
1 1
H ( z ) = 3z 2 + 2 − 5 z −1 + + .
(1 − 4 z −1 ) 2 1 − 1 z −1
2
1
[ H ( z )]+ = 2 − 5 z −1 + 1 −1
,
1− z
2
n
⎛1⎞
which corresponds to the sequence 2δ (n) − 5δ (n − 1) + ⎜ ⎟ 1(n) . Also the anticausal part of
⎝2⎠
H ( z ) is
1
[ H ( z )]− = 3z 2 + ,
(1 − 4 z −1 ) 2
25
The Causal Wiener Filter
Assume that S z ( z ) is rational in z and has no poles or zeros on the unit circle. The equation that
we must solve to obtain the causal Wiener filter is as is given in Eq. (3.4-1)
∞
∑ h( j ) R ( j − i ) = R
j =0
z sz (i ), i ≥ 0 .
⎧0, i≥0
⎪
h′(i ) = ⎨ ∞
(3.4-13)
⎪ Rsz (i ) − ∑ h( j ) Rz ( j − i), i < 0.
⎩ j =−∞
Now we can take the z-transform of both sides of Eq. (3.4-12) and obtain
26
H ′( z ) = Ssz ( z ) − H ( z ) S z ( z )
= Ssz ( z ) − H ( z ) S z + ( z ) S z − ( z ).
Dividing by Sz − ( z ) , we have
H ′( z ) Ssz ( z )
−
= − − H ( z ) Sz + ( z ) .
Sz ( z ) Sz ( z )
Extract the causal part of both sides of this equation and find
⎡ H ′( z ) ⎤ ⎡ Ssz ( z ) ⎤
⎢ − ⎥ − ⎡⎣ H ( z ) Sz ( z ) ⎤⎦ + .
+
⎢ − ⎥ = (3.4-14)
⎣ Sz ( z ) ⎦ + ⎣ Sz ( z ) ⎦ +
From Eq. (3.4-13), h′(i ) is purely anticausal. Therefore H ′( z ) contains only poles outside the
1
unit circle. Since the zeros of Sz − ( z ) lie outside the unit circle, the poles of also lie
Sz − ( z )
H ′( z )
outside the unit circle. We conclude that all the poles of are outside the unit circle, and
Sz − ( z )
thus
27
⎡ H ′( z ) ⎤
⎢ − ⎥ = 0.
⎣ Sz ( z ) ⎦ +
H ( z ) is causal by definition, so its poles lie within the unit circle. The poles of Sz + ( z ) are also
within the unit circle. Hence all the poles of H ( z ) Sz + ( z ) are inside the unit circle, and
⎡⎣ H ( z ) S z + ( z ) ⎤⎦ = H ( z ) Sz + ( z ) .
+
Ssz ( z )
We can not conclude anything about . Therefore, Eq. (3.4-14) becomes
Sz − ( z )
⎡ S ( z) ⎤
0 = ⎢ sz− ⎥ − H ( z ) Sz + ( z ) ,
⎣ Sz ( z ) ⎦ +
which produces the causal Wiener filter, described by its system function H ( z )
1 ⎡ Ssz ( z ) ⎤
H ( z) = ⎢ ⎥ . (3.4-15)
Sz + ( z ) ⎣ Sz − ( z ) ⎦ +
Example 3.4-3
28
Consider the same situation of Examples 3.2-1 and 3.3-1 in which s(n) and v(n) are
uncorrelated, zero-mean, JWSS random processes with
1 − (0.95) 2
Ss ( z ) = ,
(1 − 0.95 z −1 )(1 − 0.95 z )
and
Sv ( z) = 2 .
+ (1 − 0.7931z −1 )
S z ( z ) = 1.5477 ,
(1 − 0.95 z −1 )
29
and
(1 − 0.7931z )
Sz − ( z ) = 1.5477 .
(1 − 0.95 z )
Ssz ( z ) 0.0630
=
Sz − ( z ) (1 − 0.95 z −1 )(1 − 0.7931z )
−0.0794 z −1
=
(1 − 0.95 z −1 )(1 − 1.2608 z −1 )
0.2555 0.2555
= −
1 − 0.95 z −1 1 − 1.2608 z −1
and
⎡ Ssz ( z ) ⎤ 0.2555
⎢ − ⎥ = −1
.
⎣ z
S ( z ) ⎦+ 1 − 0.95 z
30
1 ⎡ Ssz ( z ) ⎤
H ( z) = ⎢ ⎥
Sz + ( z ) ⎣ Sz − ( z ) ⎦ +
1 − 0.95 z −1 0.2555 0.1651
= −1
• −1
= ,
1.5477(1 − 0.7931z ) 1 − 0.95 z 1 − 0.7931z −1
and
The MSE associated with this filter is computed using Eq. (3.3-3)
∞
MSE C = 1 − 0.1651∑ (0.7931)i (0.95)i = 0.3302 .
i =0
Compared to no filtering (MSE = 2), the causal Wiener filter reduces the MSE by 7.8 dB . (6.5
dB in FIR filter and 9.6 dB in noncausal filter.)
31