Personal tools
  •  
You are here: Home Electrical and Computer Engineering Stochastic Processes Programming Assignments

Programming Assignments

Document Actions
  • Content View
  • Bookmarks
  • CourseFeed

Introduction   ::   System Model   ::   Expectation   ::   Part 1   ::   Filters   ::   Part 2   ::   Submission

System Model

For brevity, we will employ an operator notation frequently used in the literature. If

$\displaystyle G(z) = \sum_{i=0}^{L_g} g_i z^{-i}.
$

we will use the notation

$\displaystyle s(t) = G(z) u(t)
$

as a shorthand for the filtering operation

$\displaystyle s(t) = \sum_{i=0}^{L_g} g_i u(t-i).
$

A known discrete-time input signal $ u(t)$ is applied to a system which is assumed to be FIR, having transfer function

$\displaystyle G(z) = \sum_{i=0}^{L_g} g_i z^{-i},
$

where we will also assume that the order of the filter $ L_g$ is known. The filtered signal $ s(t) = G(z) u(t)$ is corrupted by an additive noise signal which is assumed to be autoregressive (AR):

$\displaystyle \nu(t) = H(z) e(t)
$

where $ e(t)$ is a zero-mean, stationary, ergodic, white-noise signal with variance $ \sigma_e^2$ , and $ H(z)$ has the all-pole form

$\displaystyle H(z) = \frac{1}{1-\sum_{i=1}^{L_H} a_i z^{-i}}.$ (1)

The measured output signal is thus

$\displaystyle y(t) = s(t) + \nu(t) = G(z) u(t) + \nu(t) = G(z) u(t) + H(z) e(t).
$

\begin{center}\vbox{\input{sysidfig.pstex_t}
}\end{center}

The system identification problem to be addressed here is this: given the input/output measurements $ \{u(t),y(t)\}$ , estimate $ G(z)$ and $ H(z)$ .

The AR Model and Its Prediction

The model for the noise $ \nu(t)$ can be written another way. The IIR filter $ H(z)$ could be written (say, using long division) as

$\displaystyle H(z) =\sum_{i=0}^\infty h_i z^{-i}.
$

with $ h_0 = 1$ . (That is why we call it IIR -- it has an infinite number of coefficients in its impulse response.) So

$\displaystyle \nu(t) = \left(\sum_{i=0}^{\infty} h_i z^{-i}\right) e(t) =
e(t) + \sum_{i=1}^\infty h_i e(t-i).
$

For the development below, we will find it convenient to develop a predictor for $ \nu(t)$ . Given the sequence of measurements $ \nu(s)$ for $ s \leq t-1$ , what is the best estimate of $ \nu(t)$ ? We will denote this estimate by $ \nuhat(t\vert t-1)$ . If by ``best'' we mean best in the minimum mean-squared error sense, we want to find an estimator $ \nuhat(t\vert t-1)$ which minimizes

$\displaystyle E[(\nu(t) - \nuhat(t\vert t-1))^2].
$

We know that this will be the conditional mean:

$\displaystyle \nuhat(t\vert t-1) = E[\nu(t)\vert \nu(s), s=t-1,t-2,\ldots].
$

Note that if we know all the data $ \{\nu(s), s\leq t-1\}$ , we equivalently know all the data $ \{e(s), s \leq t-1\}$ . Thus $ \nuhat(t\vert t-1)$ can equivalently be written

$\displaystyle \nuhat(t\vert t-1) = E[\nu(t) \vert e(s), s=t-1,t-2, \ldots].
$

From the form

$\displaystyle \nu(t) = e(t) + \sum_{i=1}^\infty h_i e(t-i)
$

we see immediately that this conditional expectation is

\begin{equation*}\begin{aligned}
\nuhat(t\vert t-1) &= E[e(t) + \sum_{i=1}^\inft...
...=t-1,t-2,\ldots] \\
&=\sum_{i=1}^\infty h_i e(t-i)
\end{aligned}\end{equation*}

since $ e(t)$ has zero mean. This can be written in our operator notation as

$\displaystyle \nuhat(t\vert t-1) = [H(z)-1] e(t).
$

Since $ e(t) = H^{-1}(z) \nu(t)$ , we can write

$\displaystyle \nuhat(t\vert t-1) = \frac{[H(z)-1]}{H(z)} \nu(t) = (1-H^{-1}(z)) \nu(t).$ (2)

Let the inverse transfer function $ H^{-1}(z)$ be written in the form

$\displaystyle H^{-1}(z) = \sum_{i=0}^\infty \htilde_i z^{-1},$ (3)

where $ \htilde_0 = 1$ . (Note: It should be clear that $ \htilde_i \neq
h_i$ in general.) In the case that $ H(z)$ is in fact an all-pole filter as in ( 1 ), we have

$\displaystyle H^{-1}(z) = 1 - \sum_{i=1}^{L_H} a_i z^{-1},
$

so that $ \htilde_i = a_i, i=1,2,\ldots, L_H$ . Substituting ( 3 ) into ( 2 ) we obtain in the general case

$\displaystyle \nuhat(t\vert t-1) = -\sum_{i=1}^\infty \htilde_i \nu(t-i).
$

In the case that $ H(z)$ is an all-pole filter, we obtain

$\displaystyle \nuhat(t\vert t-1) = +\sum_{i=1}^{L_H} a_i \nu(t-i).
$

That is: the best predictor of $ \nu(t)$ given previous measurements of $ \nu(t)$ has the form of an FIR filter, as shown in the block diagram here:
\begin{center}\vbox{\input{sysidfig2.pstex_t}
}\end{center}
From a system identification perspective, we would want to choose the coefficients $ \{a_i\}$ in this filter so that the prediction error $ \nu(t) - \nuhat(t\vert t-1)$ is as small as possible. That is, we would choose $ \{a_1, a_2, \ldots, a_{L_H}\}$ so that

$\displaystyle E[(\nu(t) - \sum_{i=1}^{L_H} a_i \nu(t-i))^2]
$

is as small as possible. We can write this in a vector form as follows. Let

$\displaystyle \nubf(t) = \begin{bmatrix}\nu(t) \\ \nu(t-1) \\ \vdots \\ \nu(t-L_H+1)
\end{bmatrix}$     and  $\displaystyle \abf =\begin{bmatrix}a_1 \\ a_2 \\ \vdots \\ a_{L_H}
\end{bmatrix}.
$

Then we want to minimize

$\displaystyle E[(\nu(t) - \abf^T \nubf(t-1))^2]$ (4)


Exercise 1:
Show that minimizing the expression in ( 4 ) leads to the Widrow-Hopf normal equations

$\displaystyle R_h \abf = \pbf_h,$ (5)

where $ R_h = E[\nubf(t-1) \nubf(t-1)^T]$ and $ \pbf_h = E[\nu(t)
\nubf(t-1)]$ .

Also show that the minimum mean squared prediction error is

$\displaystyle E(\nu^2(t)) - \abf^T E[\nu(t) \nubf(t-1)],
$

where $ \abf$ is the solution to ( 5 ). This can be used as an estimate of the variance of the noise $ e(t)$ , $ \sigma_e^2$ .

Prediction of $ y(t)$ and System Identification

An important step in the overall system identification process is a one-step-ahead predictor for $ y(t)$ . Suppose that $ y(s)$ is known for $ s \leq t-1$ . What is the best estimate (prediction) that can be made for $ y(t)$ using all of this information? We will denote this predictor as $ \yhat(t\vert t-1)$ . Since the input $ u(t)$ is known, this prediction must be made on the best prediction of $ \nu(t)$ given the information up to time $ t-1$ , which we have denoted as $ \nuhat(t\vert t-1)$ :

$\displaystyle \yhat(t\vert t-1) = G(z) u(t) + \nuhat(t\vert t-1).
$

We have seen above that $ \nuhat(t\vert t-1)$ can be written as $ (1-H^{-1}(z))
\nu(t)$ . We thus have

$\displaystyle \yhat(t\vert t-1) = G(z) u(t) + (1-H^{-1}(z))\nu(t) =
G(z) u(t) + (1-H^{-1}(z))(y(t) - G(z) u(t))
$

or

$\displaystyle \yhat(t\vert t-1) = H^{-1}(z) G(z) u(t) + (1-H^{-1}(z)) y(t).
$

The prediction error is

$\displaystyle y(t) - \yhat(t\vert t-1) = -H^{-1}(z) G(z)u(t) + H^{-1}(z) y(t)
$


Exercise 2
Show that the prediction error can be written as

$\displaystyle y(t) - \yhat(t\vert t-1) = e(t).
$


The results of this exercise are very important: The prediction error for the optimal predictor form a white noise sequence. The signal $ e(t)$ that cannot be predicted from past observations represents new information. For this reason, $ e(t)$ is called the innovation at time $ t$ .

Let us take the prediction error and write it in two ways. We have

$\displaystyle e(t) = H^{-1}(z)(y(t) - G(z) u(t))$ (6)

We can also write

$\displaystyle e(t) = (H^{-1}(z) y(t)) - G(z)(H^{-1} u(t)).$ (7)

Let us consider what we can learn from each of these.


Representation 1

Suppose that $ G(z)$ were known. Then in ( 6 ) we could form the signal

$\displaystyle q(t) = y(t) - G(z) u(t),
$

With the input $ u(t)$ known and the output $ y(t)$ known, $ q(t)$ can be computed. Now rewrite ( 6 ) as

$\displaystyle e(t) = H^{-1}(z) q(t) = q(t) - a_1 q(t-1) - a_2 q(t-2) - \ldots - a_{L_H} q(t-L_H)$ (8)

The identification problem at this point is then: Determine the coefficients of the linear predictor with coefficients $ \{a_i\}$ that minimizes $ E[e(t)^2]$ . The idea is suggested by the following figure.
\begin{center}\vbox{\input{sysidfig3.pstex_t}
}\end{center}

Exercise 3
Let

$\displaystyle \abf =\begin{bmatrix}a_1 \\ a_2 \\ \vdots \\ a_{L_H}
\end{bmatrix}.
$

Show that the predictor which minimizes $ E[e(t)^2]$ in ( 8 ) has coefficients determined by

$\displaystyle R_q \abf = \rbf_q$ (9)

where

$\displaystyle R_q = E[\qbf(t) \qbf(t)^T]
$

is the autocorrelation matrix of $ \qbf(t)$ , and where

$\displaystyle \rbf_q = E[q(t) \qbf(t-1)
$

is the cross-correlation vector between $ q(t)$ and $ \qbf(t-1)$ , where

$\displaystyle \qbf(t) = \begin{bmatrix}q(t) \\ q(t-1) \\ \vdots \\ q(t-L_H+1)
\end{bmatrix}.
$


Representation 2

Now let us look at the representation in ( 7 ). Suppose that $ H^{-1}(z)$ were known, and let

$\displaystyle \ytilde(t) = H^{-1}(z) y(t)
$

$\displaystyle \utilde(t) = H^{-1}(z) u(t).
$

Then ( 7 ) can be written as

$\displaystyle e(t) = \ytilde(t) - G(z) \utilde(t).$ (10)

A system identification problem can be stated: Find $ G(z)$ to minimize

$\displaystyle E[e^2(t)].
$

The idea is suggested by the following figure.
\begin{center}\vbox{\input{sysidfig4.pstex_t}
}\end{center}


Exercise 4
Let

$\displaystyle \gbf = \begin{bmatrix}g_0 \\ g_1 \\ \vdots \\ L_g \end{bmatrix}$

be the vector of coefficients in $ G(z)$ and let

$\displaystyle \ubftilde(t) = \begin{bmatrix}\utilde(t) \\ \utilde(t-1) \\
\vdots \\ \utilde(t-L_g)
\end{bmatrix}.
$

Show that the filter $ G(z)$ which minimizes $ E[e^2(t)]$ , with $ e(t)$ as in ( 10 ), has coefficients determined by

$\displaystyle R_u \gbf = \rbf_g$ (11)

where

$\displaystyle R_u = E[\ubftilde(t) \ubftilde(t)^T]
$

is the autocorrelation matrix of $ \ubftilde(t)$ and where

$\displaystyle \rbf_g = E[\ytilde(t) \ubftilde(t)].
$

is the cross-correlation matrix between $ \ytilde(t)$ and $ \ubftilde(t)$ .


Putting the pieces together

The above sections suggest that:

  • Knowing $ G(z)$ , we can estimate $ H^{-1}(z)$ , and
  • Knowing $ H^{-1}(z)$ , we can estimate $ G(z)$ .
However, to begin with we don't know either $ G(z)$ or $ H^{-1}(z)$ . Nevertheless, all is not lost. We start with a guess of $ G(z)$ , then iterate as follows:

Pick an initial $ G(z)$ .

  1. Solve for an updated $ H^{-1}(z)$ using ( 9 ).
  2. Solve for an updated $ G(z)$ using ( 11 ).
  3. Repeat from step 1 until convergence.

A suggested $ G(z)$ is something like

$\displaystyle G(z) = 0.1 z^{-L_g/2}
$

(put a 0.1 in the middle of the impulse response).
Copyright 2008, by the Contributing Authors. Cite/attribute Resource . admin. (2006, June 13). Programming Assignments. Retrieved January 07, 2011, from Free Online Course Materials — USU OpenCourseWare Web site: http://ocw.usu.edu/Electrical_and_Computer_Engineering/Stochastic_Processes/sysid_2.htm. This work is licensed under a Creative Commons License Creative Commons License