Download LMS Adaptive Filter for System Identification in Matlab and more Assignments Stochastic Processes in PDF only on Docsity!
ECE 6010
Programming Assignment
System and Autoregressive Identification
1 Introduction
System identification is the means by which systems are modeled mathematically based on measured data. It is often a precursor to other engineering tasks, such as control or signal processing. The system models are frequently described parametrically. System identification is sometimes done using only measurements of the output of the system. In this exercise, however, it is assumed that both input and output data are available. (Usually having both input and output data significantly simplifies the system identification, often resulting in straightforward linear equations to solve.) It is usually assumed that the system measurements are made in the presence of noise. A common assumption is that the noise signals are white. However, to make things more interesting, in this assignment the noise is assumed to be from an AR process, whose parameters are also to be identified. Your task is to take measured data and from it, determine the system and noise parameters.
2 System Model
For brevity, we will employ an operator notation frequently used in the literature. If
G(z) =
∑^ Lg
i=
giz−i.
we will use the notation s(t) = G(z)u(t)
as a shorthand for the filtering operation
s(t) =
∑^ Lg
i=
giu(t − i).
A known discrete-time input signal u(t) is applied to a system which is assumed to be FIR, having transfer function
G(z) =
∑^ Lg
i=
giz−i,
where we will also assume that the order of the filter Lg 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):
ν(t) = H(z)e(t)
where e(t) is a zero-mean, stationary, ergodic, white-noise signal with variance σ e^2 , and H(z) has the all-pole form
H(z) =
∑LH
i=1 aiz −i
The measured output signal is thus
y(t) = s(t) + ν(t) = G(z)u(t) + ν(t) = G(z)u(t) + H(z)e(t).
G(z)
ν(t)
y(t)
u(t) s(t)
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).
2.1 The AR Model and Its Prediction
The model for the noise ν(t) can be written another way. The IIR filter H(z) could be written (say, using long division) as
H(z) =
∑^ ∞
i=
hiz−i.
with h 0 = 1. (That is why we call it IIR — it has an infinite number of coefficients in its impulse response.) So
ν(t) =
i=
hiz−i
e(t) = e(t) +
∑^ ∞
i=
hie(t − i).
For the development below, we will find it convenient to develop a predictor for ν(t). Given the sequence of measurements ν(s) for s ≤ t−1, what is the best estimate of ν(t)? We will denote this estimate by ˆν(t|t−1). If by “best” we mean best in the minimum mean-squared error sense, we want to find an estimator ˆν(t|t − 1) which minimizes E[(ν(t) − νˆ(t|t − 1))^2 ].
We know that this will be the conditional mean:
ˆν(t|t − 1) = E[ν(t)|ν(s), s = t − 1 , t − 2 ,.. .].
Note that if we know all the data {ν(s), s ≤ t − 1 }, we equivalently know all the data {e(s), s ≤ t − 1 }. Thus ˆν(t|t − 1) can equivalently be written
νˆ(t|t − 1) = E[ν(t)|e(s), s = t − 1 , t − 2 ,.. .].
From the form
ν(t) = e(t) +
∑^ ∞
i=
hie(t − i)
we see immediately that this conditional expectation is
νˆ(t|t − 1) = E[e(t) +
∑^ ∞
i=
hie(t − i)|e(s), s = t − 1 , t − 2 ,.. .]
∑^ ∞
i=
hie(t − i)
since e(t) has zero mean. This can be written in our operator notation as
νˆ(t|t − 1) = [H(z) − 1]e(t).
Since e(t) = H−^1 (z)ν(t), we can write
ˆν(t|t − 1) =
[H(z) − 1] H(z)
ν(t) = (1 − H−^1 (z))ν(t). (2)
Let the inverse transfer function H−^1 (z) be written in the form
H−^1 (z) =
∑^ ∞
i=
˜hiz−^1 , (3)
2.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 ≤ 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 ˆy(t|t − 1). Since the input u(t) is known, this prediction must be made on the best prediction of ν(t) given the information up to time t − 1, which we have denoted as ˆν(t|t − 1): ˆy(t|t − 1) = G(z)u(t) + ˆν(t|t − 1).
We have seen above that ˆν(t|t − 1) can be written as (1 − H −^1 (z))ν(t). We thus have
ˆy(t|t − 1) = G(z)u(t) + (1 − H−^1 (z))ν(t) = G(z)u(t) + (1 − H−^1 (z))(y(t) − G(z)u(t))
or yˆ(t|t − 1) = H−^1 (z)G(z)u(t) + (1 − H−^1 (z))y(t).
The prediction error is y(t) − ˆy(t|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
y(t) − ˆy(t|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
e(t) = H−^1 (z)(y(t) − G(z)u(t)) (6)
We can also write 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.
2.2.1 Representation 1
Suppose that G(z) were known. Then in (6) we could form the signal
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
e(t) = H−^1 (z)q(t) = q(t) − a 1 q(t − 1) − a 2 q(t − 2) −... − aLH q(t − LH ) (8)
The identification problem at this point is then: Determine the coefficients of the linear predictor with coefficients {ai} that minimizes E[e(t)^2 ]. The idea is suggested by the following figure.
G(z)
y(t)
u(t)
q(t) z−^1 Predictor
q(t − 1)^ e(t)
g
Chooise H−^1 (z) to minimize E[e^2 (t)]
ˆq(t|t − 1)
Exercise 3 Let
a =
a 1 a 2 .. . aLH
Show that the predictor which minimizes E[e(t)^2 ] in (8) has coefficients determined by
Rq a = rq (9)
where Rq = E[q(t)q(t)T^ ] is the autocorrelation matrix of q(t), and where
rq = E[q(t)q(t − 1)
is the cross-correlation vector between q(t) and q(t − 1), where
q(t) =
q(t) q(t − 1) .. . q(t − LH + 1)
2.2.2 Representation 2
Now let us look at the representation in (7). Suppose that H −^1 (z) were known, and let
y˜(t) = H−^1 (z)y(t)
u˜(t) = H−^1 (z)u(t).
Then (7) can be written as e(t) = ˜y(t) − G(z)˜u(t). (10)
A system identification problem can be stated: Find G(z) to minimize
E[e^2 (t)].
The idea is suggested by the following figure.
H−^1 (z)
G(z)
H−^1 (z) u(t)
y(t)
˜u(t)
˜y(t)
e(t)
Choose G(z) to minimize E[e(t)^2 ]
Exercise 4 Let
g =
g 0 g 1 .. . Lg
Let y(t) be a random process, with measured sample values y(1), y(2),... , y(n). The cross correlation vector between a random process y(t) and x(t), r = E[y(t)x(t)] can be estimated as
ˆr =
N
∑^ N
i=
y(t)x(t).
Exercise 5 Let X =
[
x(1) x(2) x(3)... x(N )
]
show that the estimated autocorrelation matrix Rx in (12) can be written as
Rˆx = 1 N − 1
XXT^.
4 Programming Assignment Part 1
The file sysid1 on the class website contains input/output data for an unknown system. The first column contains input u(t); the second column contains output y(t). It is known that the filter has five coefficients,
G(z) = g 0 + g 1 z−^1 + g 2 z−^2 + g 3 z−^3 + g 4 z−^4
(that is, Lg = 4) and the AR process has
H−^1 (z) = 1 + a 1 z−^1 + a 2 z−^2 + a 3 z−^3
(that is, Lh = 3). Write a Matlab program that reads in the data and, using the iterative procedure outlined in section 2.2.3, estimates G(z) and H−^1 (z). Also, determine an estimate for σ^2 e. Hints:
- You are strongly encouraged to choose your own G(z) and H −^1 (z) filters and create your own in- put/output data for debugging purposes. You can verify that your program correctly identifies the parameters in your simulated system, then use your program on the unknown data.
- Helpful Matlab commands: filter, rand, randn.
5 Adaptive filters
Solving for the unknown filter coefficients requires setting up the estimated autocorrelation and cross- correlation information, then solving a set of linear equations. These will work provided that the system is stationary, that sufficient data are available to make accurate estimates, and that the process does not take too long. (In real-time circumstances, the computation time might exceed the available computational resources.) Overcoming these demands can be accomplished in part by the use of adaptive filters. These are filters that adapt themselves to minimize some mean-squared error criterion. Here we provide a very brief introduction to the topic, by introducing what are known as LMS adaptive filters. (LMS stands for least mean squares.) Let u(t) be the input to an FIR filter, with output
x(t) =
∑^ L
i=
wiu(t − i).
The coefficients wi are the FIR filter coefficients, also referred to as filter weights. This filter operation can be represented as x(t) = wT^ u(t)
where
w =
w 0 w 1 .. . wL
is the vector of filter coefficients and
u(t) =
u(t) u(t − 1) .. . u(t − L)
is the vector of filter inputs/memory. Suppose that some desired signal d(t) is available, and we with to filter the input signal u(t) so that the filter output x(t) matches the desired signal as closely as possible. That is, we form
e(t) = x(t) − d(t),
and choose w to minimize e(t). More precisely, we desire to minimize the mean-squared error in e(t):
min w E[e^2 (t)].
Exercise 5 Show that E[e^2 (t)] is minimized when
Rw = p
where R = E[u(t)u(t)T^ ] and p = E[d(t)u(t)].
By this point, this form of the optimal solution should be familiar. So far, the filter is not adaptive. We make an adaptive filter as possible. We form an error measure as a function of the filter coefficients: J(w) = E[(d(t) − wu(t))^2 ]
Rather than minimizing all in one step, as before, we compute the gradient of J(w) with respect to w, and update the current w by moving in the direction of the negative gradient. That is, we “slide downhill” on the surface of J(w). The idea of sliding downhill is conveyed in the following figure:
Initial point
This figure shows the contours of a function (of two variables). At each point in the plane, the contours are orthogonal to the direction of the gradient, with the gradient pointing in the direction of greatest increase. Thus, starting at an initial point and moving some small distance from the point in the direction of the negative gradient at that point decreases the function value. Then starting at the new point and moving in the direction of the negative gradient at that point again decreases the function value. A series of such steps will eventually reach a point of local minimum. An update rule such as this is referred to as steepest descent.
The adaptive filter configuration is often portrayed as in the figure here, where W (z) represents the adaptive filter, W (z) = w 0 + w 1 z−^1 + · · · + wLz−L.
x(t)
− (^) e(t)
u(t)
d(t)
W (z)
The dashed line feeding back through the filter box is intended to suggest variability, such as in a variable resistor.
5.1 Adaptive predictors
The adaptive filter can be used as an adaptive predictor, that is, a predictor which learns the taps that it needs to predict its input signal.
z−^1 x(t)
− (^) e(t)
d(t)
u(t) u(t − 1) W (z)
In this configuration, the input to the filter is u(t − 1). The output of the filter is interpreted as
x(t) = ˆu(t|t − 1),
that is, the prediction of u(t) based on previous information. The desired signal is d(t) = u(t): The filter adapts itself until d(t) is as close as possible to u(t), using the information in u(t − 1), u(t − 2),... , u(t − L). Thus, the predictor in the figure in section 2.2.1 can be replaced with an adaptive filter.
5.2 Adaptive identification
The problem of identifying G(z) to minimize the error in the representation
e(t) = ˜y(t) − G(z)˜u(t)
can also be expressed using an adaptive filter, as shown here:
− (^) e(t)
˜u(t)
d(t) = ˜y(t)
G(z)
6 Programming Assignment Part 2
- Implement an LMS adaptive filter function in Matlab. You will want to write it so that the filter memory is passed in (so that you may have multiple matched filters in your code. It might be implement it so that it is declared as
[x,e, uvect,w] = function lms(u,d,uvectin,win,mu)
where:
- u and d represent the filter input and desired signal, respectively,
- uvectin represents u(t − 1) (the memory contents of the filter from the last time,
- win represents w[k], the filter weights from the last time,
- mu represents μ, the adaptive filter step size,
- x and e represent the filter output and filter error, respectively,
- uvect represents u(t) (to be used next time around)
- w represents w[k+1]^ (to be used next time around).
- Use your adaptive filter function to replace the prediction block and filtering block of the system identification. (You will need two adaptive filters.) Try various values of μ to see which give good convergence.
- Make a plot of the filter output error as a function of time for the two filters for various values of μ.
7 What to turn in
- A discussion of what you did, what you learned, and what troubles you had.
- Printouts of your program listings.
- Results of the system identification for the two different programming exercises. (That is, what are the filter taps? What is σ^2 e )
- Plots as appropriate and necessary.
Careful writing is important, as is getting the right answer! (That is, make sure you identify the system correctly!)
