Description
Note that I am looking for good MATLAB programming style; for example avoid for loops
when possible, using vectorization methods. [This doesnít mean for loops are prohibited,
but avoid silly ones].
This assignment all relates to a real WSS signal x (n) that is modeled as:
x (n) = v (n) + v (n 1) + 0:36v (n 2) + 0:6x (n 1) 0:9x (n 2)
where v (n) is 0-mean white noise with
2
v = 4.
1. Before you start coding in MATLAB:
(a) Is x AR, MA or ARMA?
(b) Is this Ölter (with v as input and x as output) the whitening or the innovations
Ölter?
(c) Write a di§erence equation for the inverse Ölter.
(d) Write an explicit formula by hand (not simpliÖed) for the PSD Sx (!) of x.
2. In MATLAB, draw the pole-zero plot for the Ölter (input v, output x).
3. In MATLAB, generate N = 10000 samples of v and then apply the Ölter to generate
N samples of x.
4. Let m0 = 5. In MATLAB, use time-averaging to estimate rx (m) for 0 m m0.
Do this in a for loop indexed by m. Rather than presenting the detailed process for
actually computing rx (m), we will use these estimated values for the remainder of the
problem as the ìtrueîvalues.
5. Consider a vector comprised of M consecutive samples of x, deÖned as:
xM (n) =
2
6
4
x (n)
.
.
.
x (n M + 1)
3
7
5
For which M will the correlation matrix for xM (n) involve rx (m) for jmj m0?
Generate this matrix in MATLAB using the values you found above (remember that
r (m) = r (m) for real random signals!) Hint: Examine the function toeplitz.
6. Check that your correlation matrix is positive deÖnite by computing its eigenval
7. The MATLAB function pwelch can be used to estimate the power spectral density; the
name of the method this function employs is called the modiÖed Welch periodogram.
Invoke pwelch as follows:
[s_est; w] = pwelch(x; hamming(512); 256; 512);
where s_est is the estimated PSD and w is the frequency vector (in normalized digital
radian units).
8. We want to compare the estimated spectrum, s_est, and the actual PSD Sx (!) computed at the points in w, on a linear (not decibel) scale. There is an issue of normalization. A simpliÖed explanation of what pwelch does is that it computes the
magnitude-squared Fourier transform (DFT actually) of blocks of x and averages the
results. But then factors such as the size of the blocks can e§ect the scaling. So for
purposes of comparison here, normalize s_est and the actual Sx (!) you compute so
they each have AVERAGE value 1 (that is, if s is the vector of PSD values, apply the
normalization s = s=mean(s); this approximates the condition 1
2
R
S (!) d! = 1).
Then superimpose graphs of them. The horizontal axis should extend from 0 to (if
you simply call plot MATLAB wonít produce thisñyou will have to adjust the graph
so it goes from 0 to only).
9. If you look at the pole-zero plot of H, you will note there is a pole pair fairly close to
the unit circle, that is in fact close to the frequency where Sx (!) has a peak. Compare
the angle of the poles of the Ölter H with the value ! where Sx (!) has a peak. They
wonít exactly match, but they should be close.
10. Generate a data matrix A whose columns are comprised of sliding blocks of x, that is:
A =
2
6
4
x (m0 + 1) x (m0 + 2) x (N)
.
.
.
.
.
.
.
.
.
x (1) x (2) x (N m0)
3
7
5
[Hint: Use toeplitz again!] Note that A is (m0 + 1) (N m0). By stationarity, each
column of A has the same statistical properties, and we can think of A as representing
a collection of (noisy) measurements with underlying statistical properties we want to
extract. Observe that each column, in fact, viewed as a stand-alone random vector,
has the same correlation matrix R found above. Now, use MATLAB svd to compute
the singular values of A, say (A) (stored in a vector). Compare the vector of values
( (A))2
= (N m0) to the eigenvalues of R. [Remark: When comparing the vectors, the
entries may occur in di§erent orders; use the sort function to overcome this problem]
[Remark: As R is a matrix of second order moments, its eigenvalues represent power;
the square of the singular values of the data matrix A similarly represent power, in a
sense; the normalization can be thought of giving us an average power related to each
column of A]. Also check that the left singular vectors of A are the eigenvectors of
R as follows: take ~u to be one of the singular vectors; do elementwise division of R~u
over ~u and all the entries should be the same (ideally) equal to an eigenvalue of R. In
fact, there may a couple cases with a mismatch but if you look closely the elements of
~u are close to 0 at those points. The moral of the story: the SVD of the data matrix
encapsulates the statistical properties of the data. :-
