Incorrect Frequency results when using Multiple Signal Classification (MUSIC)

Question

I am using MUSIC (Multiple Signal Classification) to determine Direction Of Arrival (DOA) and frequency of signals impinging on an Antenna Array. I am using a function MUSIC(U,V) which accepts left and right singular vectors of a matrix X. U is used to determine frequencies whereas V is used to determine the DOA (Method I)

However, I would like to reduce SVD computation by making X sparse. Assuming, X has no noise, taking DFT of X results in sparse matrix with peaks only at indexes relating to frequency of a signal and sine of the direction of arrival. The rest are zero. SVD() function in MATLAB is used to calculate reduced singular vectors which on which later inverse DFT is performed. Feeding them to MUSIC(S,V) gives same results as without taking the DFT.

However, Matrix obtained after DFT is not sparse as it contains some value instead of zero. It is because of precision error. If filtering is performed, although there is no noise present, these precision errors are treated as noise and are set to zero. (Method II) This results in inaccurate peaks for frequency. DOA remains unaffected. The original matrix can still be reconstructed without any noticeable difference. However, if SVD is performed after reconstructing the original matrix from svd of filtered matrix, the inaccurate peaks in frequency vanish. (Method V)

Please follow the post below for further clarification. Methods in parentheses are further explained below.

Let: $$X = \begin{pmatrix} x_{11} & x_{12} & \ldots & x_{1L}\\ x_{21} & x_{22} & \ldots & x_{2L}\\ \vdots & \vdots & \ddots & \vdots\\ x_{N1} & x_{N1} &\ldots & x_{NL} \end{pmatrix} $$

$L$ - no of antennas

$N$ - no of samples

So, columns specify antenna number, whereas row specifies sample. For example:

• $x_{11}$ - means first sample of first antenna
• $x_{1L}$ - means first sample Lth antenna
• $x_{51}$ - means fifth sample of first antenna
• $x_{NL}$ - means Nth sample of Lth antenna

Calculating SVD, [U,S,V] = svd(X), is computationally intensive so I would like to reduce X to sparse matrix. One method for doing that is to take double sided DFT and filter noise. In MATLAB:

$$ F_{f} = \frac{1}{\sqrt{N}} * \mathrm{FFT}(I_N) \\ F_d = \frac{1}{\sqrt{L}} * \mathrm{FFT}(I_L) $$

where $I_N$ and $I_L$ are square identity matrix of size $N$ and $L$ respectively.

$$ XDF = F_f \times X \times F_d $$

Lets suppose $X$ had no noise. Also, lets suppose that taking DFT doesn't result in signal leakage. Then if no of signals impinging on antenna array = $D$, with every signal having an orthogonal frequency, XDF will have peaks at $D$ locations, the rest being zero. Since a double sided DFT was performed, row index corresponds to normalised frequencies and column index corresponds to

$$2\pi \frac{\lambda}{c}\sin(\theta)$$

where $\theta$ is a direction of arrival

---

Method I

IF $X$ had originally no noise, and $D$ signals with orthogonal frequencies impinged on antenna array, XDF (matrix of $M$ x $L$) would contain $D$ non-zero values. Then [U,S,V] = svd(X,0) is calculated. Inverse Fourier Transform performed by:

U = Ff' * U
V = Fd * V

Resulting U,V are passed to MUSIC(U,V) and frequencies are determined without error.

Follwoing are graphs of frequency and sparsity of XDF when $D=3$, $N = 30$ and $L = 20$.

---

Method II

We see that despite what we expected, XDF is not sparse. The reason lies in the fact that instead of zero, there are precision errors to the order of $10^{-16}$ which nevertheless appear as a blue dot. Lets suppose a noise filter was applied that was supposed to remove all values smaller than $10^{-6}$, we get the following sparse matrix just as desired, with the precision errors exclusively forced to $0$.

However, it creates the aliasing in frequency domain:

---

Method III

However, if noise of minute order ($10^{-4}$) is re-added to XDF, aliases disappear. Adding noise is improving the MUSIC!

density_ = 0.1;
coeff = 0.0001;
randomr = sprand(N,L,density_);
randomi = sprand(N,L,density_);}
XDF = XDF + 0.0001*(full(randomr) + sqrt(-1)*full(randomi));
[U,S,V] = svd(XDF,0);

---

Method IV

However, if MATLAB function svds (svd for sparse) is used, which calculates using evds (eigen values for sparse) which in turn uses Lanczos Algorithm, we get:

---

Method V

If $$ [U,S] = XDF\times XDF^H \\ [V,S] = XDF^H \times XDF \\ U = F_{f}^{H} \times U \\ V = F_{d} \times V $$

MUSIC(U,V) also gives the same result as above, i.e. method IV. It can be argued that a possible reason for inaccurate results in II are errors induced by setting precision errors to zero. However, X calculated by multiplying U*S*V' is only slightly different from the original X. On the order of $10^{-16}$. Also, if X is again decomposed into its singular vectors, we again get the accurate results as in III.

---

My question is:

1. Problem in II gets solved if we make $L>M$. Make number of columns more than number of rows.
2. In II, if instead of economical SVD, full SVD is used, problems still disappear
3. In II, if [U,S,V] = svd(XDF,0) and next step is U = Ff*U', V = Fd*V;, abs(U) shows that magnitude of every entry in U is equal to $1/\sqrt{N}$. This is not the case for I, III and IV.
4. I tried to try other methods of calculating svd such as Jacobi SVD. However, since XDF is sparse, Jacobi Rotation matrices are zero and division by zero results in undefined U.
5. I have considered using pre-developed C++ libraries such as Eigen, Alglib, but creating their MEX file is really cumbersome. Some libraries such as ARmadillo and redsvd are already based on LAPACK - same as MATLAB so it is useless.

I would really appreciate if anyone could help me regarding in these issues. I would love to narrow down my problem by trying other SVD algorithms but like i said, it is very time consuming to make MEX file for all those C++ libraries. Then, I am not even sure if the problem lies in different methods of SVD calculation. The singular matrices themselves are not incorrect. Since, they successfully reconstruct the initial matrix X.

However, somehow, in II, for sparse matrix it requires full left singular matrix to identify frequencies. Furthermore, this problem is faced only when finding frequencies and not for DOA.

Below is a description of how MUSIC(U,V) works:

I am using a function MUSIC(U,V) to determine frequencies and DOA. $$[U_{s} U_{n}] = U \\ [V_{s} V_{n}] = V$$
if number of incoming signals is $D$, size of $$ U = M \text{ x } L\\ V = L \text{ x } L \\ U_{s} = M \text{ x } D \\ U_{n} = M \text{ x } (L-D) \\ V_{s} = L \text{ x } D \\ V_{n} = L \text{ x } (L-D) $$

Since, signal and noise singular vectors are supposed to be orthogonal to each other, Vector $a$, defined as $$ a = \begin{bmatrix} e^{-j2\pi f0} \\ e^{-j2\pi f1} \\ \vdots \\ e^{-j2\pi f(M-1)} \end{bmatrix} $$ If $f$ is a frequency of one of the arriving signals: $$ a\cdot U_n = 0, \text{ and } \frac{1}{a\cdot U_n} = \infty $$

Hence, what MUSIC(S,V) does is search for peaks over frequency range from $-0.5$ to $0.5$. $f$ is a normalised frequency.