Problem Set I: Linear Algebra Solution

Note: Unless stated otherwise, consider only …nite dimensional linear spaces. Also assume complex values.

Aside from basic MATLAB functions, in this homework you are allowed to use: rank, orth, svd, svds, eig, eigs, pinv, chol. Use doc to familiarize yourself with these functions. Given that these functions are available to you, I want you to use e¢ cient code. Do NOT …nd more advanced or more sophisticated functions.

Notation: For a matrix A, RA represents the range space and NA the null-space. If K is a linear space, PK is the orthogonal projection onto that space. Matrix inversion lemma:
B 1+CD 1CH  1 =B  BC D+CHBC  1CHB

1. Let A be a square invertible matrix, a scalar and u a vector. Use the matrix inversion lemma to simplify:

A +    u uH   1

For certain values of , this matrix may not be invertible; assume this is not the case. Now take the special case where A = I. The answer reduces to something of the form c1I + c2u uH . Find this expression. Also, for this case, determine the condition on for this to be invertible.

    2. Let fuigri=1 be r linearly independent M 1 vectors, and fvigri=1 be r linearly inde-pendent N 1 vectors. Show that the M N matrix A given by:
Xr
A =    uiviH

i=1

has rank exactly r. Hint: Think about RA and NA?.

3. Let ~x; ~u be M    1 column vectors. We want to …nd optimal    such that ~x    ~u.

(a) Find u#. Use your formula to show that:
= hx; ui

juj2

    (b) Let’s generalize this. First: if D = diag fdig, does DA or AD scale the columns of A by di’s? So then what scales the rows?

    (c) Let U be an M N matrix with orthogonal (not necessarily orthonormal) columns fuig. Find the SVD for U, from that …nd U#, and use this to solve the following problem: given x, …nd coe¢ cients f ig such that:
X
x    iui

Note: You should know the answer. Here we are using the pseudo-inverse to con…rm it.

1
    4. Let U be a linear subspace of CM , and P be the M M projection matrix onto U. It can be shown that P = P H and P = P 2 (I’m not asking you to prove this).

        (a) Prove the converse: if P is an M M matrix that is Hermitian with P = P 2, then it is the projection matrix onto a subspace; this space in fact is the _____

space of P . Hint: First prove that, for any x; y, P x ? (I   P ) y.

        (b) Let P be the projection onto U. Consider H = I 2P . The expression Hx is called the Householder transformation of x. Consider the expansion x = x1 + x2, where x1 2 U, x2 2 U?. Express Hx in terms of x1 and x2. Note: In the special case where U is 1-dimensional, say of the form f v : 2 Cg for a unit vector v, then the Householder transformation is the re‡ection across the hyperplane that is ? v.


5.

Sensor Array Signal Model


You will be implementing a number of algorithms for array processing in this course. To prepare for that, this problem has you write the simulation code to generate the signals received at the output of the sensor array, and then apply basic SVD analysis. Note we will assume ALL random signals are 0-mean.

First some preliminary comments regarding Gaussian noise. If n = nI + jnQ is scalar
Gaussian, its variance E
jnj2

= E
jnI j2  + E

jnQj2
(even if nI ; nQ are correlated!).
Thus, for example, the
following will generate an M

N matrix, where every entry is










UNIT variance complex Gaussian, and all entries are independent:

    • = 1=sqrt(2)  (rand n (M; N) + j  rand n (M; N))

Suppose we want to create an M    1 random Gaussian vector with covariance matrix

C, and mean    . If n is a unit variance white (complex) Gaussian vector, then:

    • = C1=2n +

will yield x N ( ; C) where C1=2 could either be the Cholesky factor or Hermitian square root of C.

We consider multiple plane wave sources incident on a sensor array with M isotropic elements. Let f~rigMi=1 be the coordinates of the array. The parameter vector = ( ; ) is called the angle of arrival (AOA) or direction of arrival (DOA), where we use spherical coordinates; is called the polar angle (0 =2) and the azimuthal angle (0 2 ), and a unit vector in direction is given by:

a^ (  ) = sin  cos  a^x + sin  sin  a^y + cos a^z

The wavenumber vector then has the form:

~
2
k (  ) =

a^ (  )



where here we consider the narrowband case only, i.e.    is a constant.  Then the

steering vector is:
2



.



3










s (  ) =

exp

jk~.
(  ) ~r1




6


~.




7

6







7

4






5

6
exp


jk (  )
~rM

7
2
Note this is not a unit vector, but p1M s ( ) is (i.e., the necessary scaling does not depend on ). Assume all coordinates ~ri are expressed as multiples of an underlying spacing d; in this case, s ( ) will depend on the factor d= ; assume d= is given as an external parameter. Each snapshot (instant in time n) is an M 1 vector, comprised of L sources, plus noise:






Xi




1
L











u [n] =
p


i [n] s (  i) + v [n]





M
=1



The Gaussian noise v [n] is white across both space and time, i.e. E
v [n] vH [n]
=
v2IM  M , and for n 6= n0, E
v [n] vH [n0]
= 0M  M . At each time,  i [n]
for i   n
L,
are scalars, representing
the signals associated with each source.  The relationships










across space and time for i [n] can be more complicated. Clearly it is not unreasonable to assume they are NOT white across time; more generally, each could be a WSS signal with some presumed correlation function. Also, the sources are not necessarily uncorrelated with each other; for example, in a multipath environment, re‡ections may cause the same underlying signal to arrive from multiple directions upon the array. However, at least for now, we are going to keep things simple. We will assume the source signals i [n] are uncorrelated with each other, and are white Gaussian across time, with respective variances 2i, 1 i L. The normalization p1M is there so 2i and

2v can be properly viewed as signal and noise powers respectively. If S is the M L matrix with normalized (unit length) steering vectors as columns, b [n] the L 1 vector with entries fbi [n]g, then each snapshot can be written as:

u [n] = Sb [n] + v [n]

To avoid confusion between the signal variances 2i and the singular values of A you will be computing later (which are NOT the same), let’s use S1; S2; to denote the singular values. (If we use si for singular values that looks too much like the steering vectors! Ugh)

    (a) Write code that takes the given parameters as input, and outputs the steering matrix S, and the data matrix A, whose columns are u [n] for 1 n N. I expect EFFICIENCY and CLARITY in the code. The ONLY for loop you are allowed is one loop over the source index 1 i L to generate the S matrix. Once you have S, you must compute A with no for loops.

    (b) Determine the formula for correlation matrix R of u [n], and write code to compute
it. Also, as speci…ed, either K1 AH A or K1 AAH is an unbiased estimate of R, for

^
some choice of K. Via ergodic arguments, this can be used as an estimate R of
^
R. Determine the right answer, and write code to compute R as suggested.

(c) Take d= = 0:5, and consider an array in the form of a cross of two linear arrays aligned along the x- and y-axes, respectively. That is, take locations at (md; 0; 0) and (0; md; 0) where 10 m 10 (Don’t count the origin twice!) Assume three sources, with normalization 21 = 1; the other two sources are 5dB and


3
10dB below the primary source, and the noise is 20dB below the primary source.

Assume AOAs as follows:

1
=
10 ;
= 20


1

2
= 20; 2= 20
3
=
30 ;
= 150


3

^
Take N = 100 snapshots and form the A matrix. Also compute R and R for this case.

1. Eigenanalysis of R can be used to determine "theoretical" values for the sin-gular values and left singular vectors of A (obviously, since A is random, these are not what you would actually get from SVD of A; but we can interpret these as "ideal" representations for the signal and noise subspaces). Do the eigenanalysis to …nd the complete set of "theoretical" singular values and left singular vectors. Graph these computed singular values (as a stem plot), and con…rm that there are three dominant singular values. Compute the projec-tion matrices PS; PN onto the signal space and noise subspaces, respectively, from the singular vectors, and check that the singular vectors lie in the signal space by computing jPN s ( i)j, which should be 0 but for numerical precision errors.

2. Consider the case where the white noise v is not present, i.e., v [n] = 0. This yields noise-free data matrix A0 and correlation matrix R0. Compare R0 and R and comment what impact the noise has on the computed singular values and singular vectors.

3. Now perform SVD on A. Graph the singular values. Check that S4 appears much lower that S3, which would be expected because we have a good SNR.
4. Let’s compute SNR from the singular values. Based on your analysis above,

you should expect 2S4 2v. What do you expect P3i=1 2Si to be? Use this to compute a total SNR and compare to what it should be (i.e., add the signal
powers- since they are uncorrelated this is reasonable- and divide by the noise power).
5. You have a theoretical value for the R matrix of u. On the other hand, you can

 
^    ^
estimate R from the data matrix A. Compute  R    R . In order to interpret

how large this is, it would be reasonable to compare it to the noise power that is present. This is also equivalent to comparing it to a certain eigenvalue of R.
^
Explain, and do this (i.e., compare). Note: Because computing R involves an averaging process, hopefully this error will be SMALLER than the intrinsic noise ‡oor.

    6. Using the ESTIMATED correlation matrix, compute the MVDR spectrum, and use the SVD of A to …nd the MUSIC spectrum. When computing the spectrum, use unit length steering vectors. Evaluate each on a grid ( ; ) of 2 , and obtain both contour and surface plots of each.

    7. Since the MUSIC spectrum and MVDR spectrum are computed for steer-ing vectors, not "general" vectors, …nding their theoretical minimum can be complicated. But we can at least get a lower bound: assuming we evaluate

4
these spectra for UNIT LENGTH steering vectors, obtain a lower bound for each. That is, if you are allowed to plug in ANY unit vector (not necessarily a steering vector), what would be the MINIMUM for each?

    8. Plug in the three source steering vectors into the MVDR and MUSIC spectra, and compare to the theoretical lower bounds found above. You expect a peak at the source steering vectors: I want you see how tall the peak is relative to the theoretical minimum values.

    9. Actually, since you are computing the MVDR and MUSIC spectra on a grid of ( ; ) points, …nd the overal minimum values you computed, and compare to the theoretical lower bounds you found. Are they reasonably close?

    10. You have M = 21 sensors in your array, and used N = 100 snapshots. Repeat the above with N = 25, and comment on what impact, if any, this reduction has on your results; in practice, people usually prefer to use N > 2M; here we have N M (though N M still holds). Repeat, again, if N = 10. It may be that something "fails" completely if N < M: if yes, report it. One could argue that, since N > L, we may still be able to identify the sources.

    11. Reset back to N = 100. This time, however, make 2 = 10 , 2 = 10 . This brings the …rst two sources close together. Repeat your experiment and comment on any e¤ect this may have had.








































5