An Introduction to Statistical Signal Processing

Comments. Two random variables U and V are said to be orthogonal if E(UV ) = 0. Therefore equation (5.53) is known as the orthogonality principle because it states that the optimal filter causes the estimation error to be orthogonal to the observations. Note that (5.53) implies not only that the estimation error is orthogonal to the observations, but that it is also orthogonal to all linear combinations of the observations. Relation (5.56) with K = (−∞, n) is known as the Wiener-Hopf equation. To be useful in practice, we must be able to find a pulse response that solves one of these equations. We shall later find solutions for some simple cases. A more general treatment is beyond the intended scope of this book. Our emphasis here is to demonstrate an example in which determination of an optimal filter for a reasonably general problem requires the solution of an equation given in terms of second-order moments.

Proof. Suppose that we have a filter h that satisfies the given conditions. Let g be any other linear filter with the same input observations and let

˜n be the resulting estimate. We will show that the given conditions imply

X

that g can yield an expected squared error no better than that of h. Let

− ˜

*˜n = Xn Xn be the estimation error using g so that

Add and subtract the estimate using h satisfying the conditions of the theorem and expand the square to obtain

i: n−i K

The first term on the right is the expected squared error using the filter h, say E(*2n). The last term on the right is the expectation of something squared and is hence nonnegative. Thus we have the lower bound

where we have brought one of the sums out, used di erent dummy variables for the two sums, and interchanged some expectations and sums. From

(5.54), however, the bracketed term is zero for each i in the index set being summed over, and hence the entire sum is zero, proving that

E(˜*2n) ≥ E(*2n) ,

which completes the proof of the theorem.

Note that from (5.52) through (5.56) we can write the mean square error

k: n−k K

for weakly stationary processes.

Older proofs of the result just given use the calculus of variations, that is, calculus minimization techniques. The method we have used, however, is simple and intuitive and shows that a filter satisfying the given equations actually yields a global minimum to the mean squared error and not only a local minimum as usually obtained by obtained by calculus methods. A popular proof of the basic orthogonality principle is based on Hilbert space methods and the projection theorem, the generalization of the standard geometric result that the shortest line from a point to a plane is the projection of the point on the plane — the line passing through the point which meets the plane at a right angle (is orthogonal to the plane). The projection method also proves that the filter of (5.53) yields a global minimum.

We consider four examples in which the theorem can be applied to construct an estimate. The first two are fairly simple and su ce for a brief reading.

[5.4] Suppose that the processes are jointly weakly stationary, that we are given the entire two-sided realization of the random process {Yn}, and that there are no restrictions on the linear filter h. Equation (5.56) then becomes

RX,Y (k) = hiRY (k − i) ; all k Z .

i Z

This equation is a simple convolution and can be solved by standard Fourier techniques. Take the Fourier transform of both sides and define the transform of the cross-correlation function RX,Y to be the cross-spectral density SX,Y (f). We obtain SX,Y (f) = H(f)SY (f) or

H(f) = SX,Y (f) ,

SY (f)

which can be inverted to find the optimal pulse response h:

h(k) = 1/2 SX,Y (f) ei2πkf df ,

−1/2 SY (f)

which yields an optimum estimate

Thus we have an explicit solution for the optimal linear estimator for this case in terms of the second-order properties of the given processes. Note, however, that the resulting filter is not causal in general. Another important observation is that the filter itself does not depend on the sample time n at which we wish to estimate the X process; e.g., if we want to estimate Xn+1, we apply the same filter to the shifted observations; that is,

Thus in this example not only have we found a means of estimating Xn for a fixed n, but the same filter also works for any n. When one filter works for all estimate sample times by simply shifting the observations, we say that it is a time-invariant or stationary estimator. As one might guess, such time invariance is a consequence of the weak stationarity of the processes.

The most important application of example [5.4] is to “infinite smoothing,” where Yn = Xn + Vn. {Vn} is a noise process that is uncorrelated with the signal process {Xn}, i.e., RX,V (k) = 0 for all k. Then RX,Y = RX and hence RY = RX + RV , so that

[5.5] Again assume that the processes are jointly weakly stationary. Assume that we require a causal linear filter h but that we observe

the infinite past of the observation process. Assume further that the observation process is white noise; that is, RY (k) = N20 δk. Then

K = {n, n − 1, n − 2, . . . }, and equation (5.56) becomes the WienerHopf equation

This equation easily reduces because of the delta function to

2

hk = N0 RX,Y (k) , k Z+ .

Thus we have for this example the optimal estimator

As with the previous example, the filter does not depend on n, and hence the estimator is time-invariant.

The case of a white observation process is indeed special, but it suggests a general approach to solving the Wiener-Hopf equation, which we sketch next.

[5.6] Assume joint weak stationarity and a causal filter on a semi-infinite observation sequence as in example [5.5], but do not assume that the observation process is white. In addition, assume that the observation process is physically realizable so that a spectral factorization of the form of (5.40) exists; that is,

SY (f) = |G(f)|2

for some causal stable filter with transfer function G(f). As previously discussed, for practical purposes, all random processes have spectral densities of this form. We also assume that the inverse filter, the filter with transfer function 1/G(f), is causal and stable. Again, this holds under quite general conditions. Observe in particular that you can’t run into trouble with G(f) being zero on a frequency interval of nonzero length because the condition in the spectral factorization theorem would be violated.

Unlike the earlier examples, this example does not have a trivial solution. We sketch a solution as a modification to the solution of example

[5.5]. The given observation process may not be white, but suppose that we pass it through a linear filter r with transfer function R(f) = 1/G(f) to obtain a new random process, say {Wn}. Since the inverse filter 1/G(f)is assumed stable, then the W process has power spectral density SW (f) = SY (f)|R(f)|2 = 1 for all f; that is, {Wn} is white. One says that the W process is a whitened version of the Y process, sometimes called the innovations process of {Yn}. Intuitively, the W process contains the same information as the Y process from it by passing it through the filter G(f) (at least in principle). Thus we can get an estimate of Xn, from the W process that is just as good as (and no better than) that obtainable from the Y process. Furthermore, if we now filter the W process to estimate Xn, then the overall operation of the whitening filter followed by the estimating linear filter is also a linear filter, producing the estimate from the original observations. Since the inverse filter is causal, a causal estimate based on the W process is also a causal estimate based on the Y process.

Because W is white, the estimate of Xn from {Wn} is given immediately by the solution to example [5.5]; that is, the filter h with the W process as input is given by hk = RX,W (k) for k ≥ 0. The cross-correlation of the X and W processes can be calculated using the standard linear filter I/O techniques. It turns out that the required cross-correlation is the inverse Fourier transform of a cross-spectral density given by

where the asterisk denotes the complex conjugate. (See problem 5.21.) Thus the optimal causal linear estimator given the whitened process is

and the overall optimal linear estimate has the form shown in Figure 5.3. Although more complicated, we again have a filter that does not depend

on n.

This approach to solving the Wiener-Hopf equation is called the “prewhitening” (or “innovations” or “shaping filter”) approach and it can be made rigorous under quite general conditions. That is, for all practical purposes, the optimal filter can be written in this cascade form as a whitening filter followed by a LLSE filter given the whitened observations as long as the processes are jointly weakly stationary and the observation process’s power spectral density.

When the observation interval is finite or when the processes are not jointly weakly stationary, the spectral factorization approach becomes quite

Figure 5.3: Prewhitening method

complicated and cumbersome, and alternative methods, usually in the time domain, are required. The final example considers such an estimator.

[5.7] Suppose that the random process we wish to estimate satisfies a difference equation of the form

where the process {Un} is a zero-mean process that is uncorrelated with a possibly time-varying second moment E(Un2) = Γn and X0 is an initial random variable, independent of the {Un}. {Φn} is a known sequence of constants. In other words, we know that the random process is defined by a time-varying linear system driven by noise that is uncorrelated but not necessarily stationary. Assume that the observation process has the form

a scaled version of the X process plus observation noise, where Hn is a known sequence of constants. We also assume that the observation noise has zero mean and is uncorrelated but not necessarily stationary, say E(Vn2) = Ψn. We further assume that the U and V processes are uncorrelated: E(UnVk) = 0 for all n and k. Intuitively, the random processes are such that new values are obtained by scaling old values and adding some perturbations. Additional noise influences our observations. Suppose that we observe Y0, Y1, . . . , Yn−1, what is the best linear estimate of Xn?

In a sense this problem is more restrictive than the Wiener-Hopf formulation of example [5.6] because we have assumed a particular structure for the process to be estimated and for the observations. On the other hand,

it is not a special case of the previous model because the time-varying parameters make it nonstationary and because we restrict the observations to a finite time window (not including the current observation), often a better approximation to reality. Because of these di erences the spectral techniques of the standard Wiener-Hopf solution of example [5.6] do not apply without significant generalization and modification. Hence we consider another approach, called recursive estimation of Kalman-Bucy filtering, whose history may be traced to Gauss’s formula for plotting the trajectory of heavenly bodies. The basic idea is the following: Instead of considering how to operate on a complete observation record in order to estimate something at

2

one time, suppose that we already have a good estimate Xn for Xn and that

we make a single new observation Yn. How can we use this new information

2 to update our old estimate in a linear fashion to form a new estimate Xn+1

of Xn+1? For example, can we find sequences of numbers an and bn so that

is a good estimate? One way to view this is that instead of constructing a filter h described by a convolution that operates on the past to produce an estimate for each time, we wish a possibly time-varying filter with feedback that observes its own past outputs or estimates and operates on this and a new observation to produce a new estimate. This is the basic idea of recursive filtering, which is applicable to more general models that that considered here. In particular, the standard developments in the literature consider vector generalizations of the above di erence equations. We sketch a derivation for the simpler scalar case.

We begin by trying to apply directly the orthogonality principle of (5.53) through (5.55). If we fix a time n and try to estimate Xn by a linear filter as

then the LLSE filter is described by the time-dependent pulse response, say h(n−1), which, from (5.55), solves the equations

where the superscript reflects the fact that the estimate is based on observations through time n − 1 and the fact that for this very nonstationary problem, the filter will likely depend very much on n. To demonstrate this,

322 CHAPTER 5. SECOND-ORDER MOMENTS

consider the estimate for Xn+1. In this case we will have a filter of the form

Note that (5.64) is di erent from (5.62), and hence the pulse responses satisfying the respective equations will also di er. In principle these equations can be solved to obtain the desired filters. Since they will in general depend on n, however, we are faced with the alarming possibility of having to apply for each time n a completely di erent filter h(n) to the entire record of observations Y0, . . . , Yn up to the current time, clearly an impractical system design. We shall see, however, that a more e cient means of recursively computing the estimate can be found. It will still be based on linear operations, but now they will be time-varying.

We begin by comparing (5.62) and (5.64) more carefully to find a relation between the two filters h(n) and h(n − 1). If we consider l < n, then (5.59) implies that

RX,Y (n + 1, l) = E(Xn+1Yl) = E((ΦnXn + Un)Yl)

= ΦnE(XnYl) + E(UnYl) = ΦnE(XnYl) = ΦnRX,Y (n, l) .

Re-indexing the sum of (5.64) using this relation and restricting ourselves to l < n then yields

But for l < n we also have that

11

=Hn E(YnYl) = Hn RY (n, l)

which has the desired form. It remains, however, to find a means of computing the numbers h(1n). Since this really depends on only one argument n, we now change notation for brevity and henceforth denote this term by

xn; that is,

xn = h(1n) .

To describe the estimator completely we need to find a means of computing xn and an initial estimate. The initial estimate does not depend on any observations. The LLSE estimate of a random variable without observations is the mean of the random variable (see, e.g., problem 4.23). Since by

2

assumption the processes all have zero mean, X0 = 0 is the initial estimate. Before computing xn, we make several remarks on the estimator and its

properties. First, we can rewrite the estimator as

2

∆

estimate of Xn given Y0, Y1, . . . , Yn−1, then Yn −HnXn = Yn −Yn = νn can be interpreted as the “new” information in the observation Yn in the sense

This can be interpreted as saying that the new estimate is formed from the old estimate by using the same transformation Φn used on the actual samples and then by adding a term depending only on the new information.

It also follows from the orthogonality principle that the sequence νn is

2

uncorrelated: Since HnXn is the LLSE estimate for Yn based on Y0, . . . , Yn−1, the error νn must be orthogonal to past Yl from the orthogonality principle. Hence νn must also be orthogonal to linear combinations of past Yl and hence also to past νl. It is straightforward to show that Eνn = 0 for all n and hence orthogonality of the sequence implies that it is also uncorrelated (problem 5.24). Because of the various properties, the sequence {νn} is called the innovations sequence of the observations process. Note the analog with example [5.6], where the observations were first whitened to form innovations and then the estimate was formed based on the whitened version.

Observe next that the innovations and the estimation error are simply related by the formula

or

νn = Hn*n + Vn ,

a useful formula for deriving some of the properties of the filter. For example, we can use this formula to find a recursion for the estimate error:

*0 = X0

This formula implies that

E(*n) = 0 ; n = 0, 1, . . . ,

and hence the estimate is unbiased (i.e., an estimate having an error which has zero mean is defined to be unbiased). It also provides a recursion for finding the expected squared estimation error:

E(*2n+1) = (Φn − κnHn)2E(*2n) + E(Un2) + κ2nE(Vn2) ,

where we have made use of the assumptions of the problem statement, viz., the uncorrelation of Un and Vn sequences with each other, with Y0, . . . , Yn−1,