2
$\begingroup$

I'm working on simulating the MUSIC algorithm in MATLAB to determine the azimuth angle of arrival of a simulated signal. My virtual antenna array consists of 4 antennas arranged in a plane.

While I have the algorithm "working" in the sense that MUSIC for the most part picks a good estimate of the angle of arrival, I noticed that my correlation matrix is of a higher rank than what the algorithm claims. This leads to one fewer non-zero eigenvalue, thus one fewer eigenvector corresponding to the noise subspace. Moreover, this additional non-zero eigenvalue is not close to zero, and taking the corresponding eigenvector as part of the noise subspace greatly degrades MUSIC's performance.

I've noticed the same behavior when I look at the 2D case; i.e., when searching for azimuth and elevation angles. I've copied my code below. Can anyone shed some light on what is going on here?

% Assume 0 degrees due north, with positive clockwise orientation
% Positive phase difference indicates that antennna 1 recieves the waveform
% before the respective antenna (2,3, or 4)
% 
% Fix antenna 1 at the origin. From the x-axis, rotate antenna 2 60 degrees,
% rotate antenna 3 90 degrees, and rotate antenna 4 120 degrees. The
% optimal antenna spacing is still to be determined.
% 
%                         | 0 deg         --> +
%                         |
%                         |
%                         |
%             ---------- A1 ----------
%                   -60 ' | `  60
%                      '  |  `
%                     '   A3  `
%                  A4'         `A2


num_antennas = 4;
num_sig_sources = 1;
RF_freq = 600e6;        % freq of incoming signal
lambda = 3e8/RF_freq;   % wavelength of incoming signal
fs = 100e6;             % sampling frequency
samples = 200;          % number of snapshots
cos_freq = 1e6;         % freq to downconvert to

ant_spacings = [.25, .98, .91]; % antenna spacing

% generate_lookup_table just calculates a set of theoretical phase
% differences corresponding to some random AoA for the purpose of signal
% generation.
n = randi([0,360],1);
[AoAs, PDs1to2, PDs1to3, PDs1to4] = generate_lookup_table(-90, 90, .5, lambda, ant_spacings(1), ant_spacings(2), ant_spacings(3));
AoA = AoAs(n);

% generate signal
t = linspace(0,samples/fs,samples);
recived_signal_ant_1 = cos(2*pi.*t.*cos_freq);
recived_signal_ant_2 = cos(2*pi.*t.*cos_freq + PDs1to2(n));
recived_signal_ant_3 = cos(2*pi.*t.*cos_freq + PDs1to3(n));
recived_signal_ant_4 = cos(2*pi.*t.*cos_freq + PDs1to4(n));


% build sample correlation matrix
cov_mat = zeros(4,4);
for i = 1:samples
    X = [recived_signal_ant_1(i);
         recived_signal_ant_2(i); 
         recived_signal_ant_3(i);
         recived_signal_ant_4(i)];
    cov_mat = cov_mat + X*X';
end
cov_mat = cov_mat*(1/samples);

% grab eigenvalues and eigenvectors
[eig_vecs, eig_vals] = eig(cov_mat);

% generate matrix w/ noise eigenvector columns. An eigenvector corresponds
% to noise if the associated eigenvalue is close to 0. 
col_num = 1;
for i = 1:4
    if eig_vals(i,i) < 1e-4
        noise_vecs(:,col_num) = eig_vecs(:,i);
        col_num = col_num + 1;
    end
end

% walk t
theta = linspace(-90,90,200);
spectrum_func = zeros(1,100);
for i = 1:length(theta)

    % steering vector - only taking real part as of now (Re(e^i*phasediff) = cos(phase_diff))
    steering_vec = [1;
        cos((-2*pi*ant_spacings(1)/lambda)*sin(theta(i)*(pi/180) - pi/3));
        cos((-2*pi*ant_spacings(2)/lambda)*sin(theta(i)*(pi/180) - pi/2));
        cos((-2*pi*ant_spacings(3)/lambda)*sin(theta(i)*(pi/180) - 2*pi/3))];
    spectrum_func(i) = 1/(((steering_vec'*noise_vecs)*noise_vecs')*steering_vec);
end

% grab theta corresponding to highest peak, round to nearest half integer
[mx,argmax] = max(spectrum_func);
estimated_AoA = floor(2*theta(argmax)+.5)/2;
fprintf("Estimated AoA: %.1f \t Hard-coded AoA: %.1f.\n", estimated_AoA, AoA)

plot(theta,spectrum_func)
title('MUSIC Psuedospectrum')
xlabel('Angle of Arrival')





EDIT: I've added the source code for generate_lookup_table():


function [AoAs, PDs1to2, PDs1to3, PDs1to4] = generate_lookup_table(start, stop, step, lambda, d1to2, d1to3, d1to4)
% d1to2 = spacing between antennas 1 and 2

num_pts = (stop-start)/step + 1;
AoAs = zeros(1,num_pts);
PDs1to2 = zeros(1,num_pts);
PDs1to3 = zeros(1,num_pts);
PDs1to4 = zeros(1,num_pts);

for i = 1:num_pts 
    AoAs(i) = start + (i-1)*step;

    % phase differences
    value_to_push_to_2 = (-2*pi*d1to2/lambda)*sin(AoAs(i)*(pi/180) - pi/3);
    value_to_push_to_3 = (-2*pi*d1to3/lambda)*sin(AoAs(i)*(pi/180) - pi/2);
    value_to_push_to_4 = (-2*pi*d1to4/lambda)*sin(AoAs(i)*(pi/180) - 2*pi/3);

    % increment/decrement until in -180 to 180
    if value_to_push_to_2 > 2*pi
        n = ceil(value_to_push_to_2/(2*pi) - 1);
        value_to_push_to_2 = value_to_push_to_2 - 2*pi*n;
    end

    if value_to_push_to_2 > pi
        value_to_push_to_2 = value_to_push_to_2 - 2*pi;
    end

    if value_to_push_to_3 > 2*pi
        n = ceil(value_to_push_to_3/(2*pi) - 1);
        value_to_push_to_3 = value_to_push_to_3 - 2*pi*n;
    end

    if value_to_push_to_3 > pi
        value_to_push_to_3 = value_to_push_to_3 - 2*pi;
    end

    if value_to_push_to_4 > 2*pi
        n = ceil(value_to_push_to_4/(2*pi) - 1);
        value_to_push_to_4 = value_to_push_to_4 - 2*pi*n;
    end

    if value_to_push_to_4 > pi
        value_to_push_to_4 = value_to_push_to_4 - 2*pi;
    end

    if value_to_push_to_2 < -2*pi
        n = ceil(-1 - value_to_push_to_2/(2*pi));
        value_to_push_to_2 = value_to_push_to_2 + 2*pi*n;
    end

    if value_to_push_to_2 < -pi
        value_to_push_to_2 = value_to_push_to_2 + 2*pi;
    end

    if value_to_push_to_3 < -2*pi
        n = ceil(-1 - value_to_push_to_3/(2*pi));
        value_to_push_to_3 = value_to_push_to_3 + 2*pi*n;
    end

    if value_to_push_to_3 < -pi
        value_to_push_to_3 = value_to_push_to_3 + 2*pi;
    end
    if value_to_push_to_4 < -2*pi
        n = ceil(-1 - value_to_push_to_4/(2*pi));
        value_to_push_to_4 = value_to_push_to_4 + 2*pi*n;
    end

    if value_to_push_to_4 < -pi
        value_to_push_to_4 = value_to_push_to_4 + 2*pi;
    end

    if value_to_push_to_4 < -2*pi
        n = ceil(-1 - value_to_push_to_4/(2*pi));
        value_to_push_to_4 = value_to_push_to_4 + 2*pi*n;
    end

    if value_to_push_to_4 < -pi
        value_to_push_to_4 = value_to_push_to_4 + 2*pi;
    end

    PDs1to2(i) = value_to_push_to_2;
    PDs1to3(i) = value_to_push_to_3;
    PDs1to4(i) = value_to_push_to_4;
end
end
```
$\endgroup$
  • $\begingroup$ Do you mind providing the source for generate_lookup_table()? $\endgroup$ – Envidia Nov 15 '19 at 17:43
  • $\begingroup$ Are you sure you are simulation the signal response based on an antenna array correctly? I don't see where you use a 4-element antenna response sum to gather the signal returns. What I do see is that you simply assign a sinusoid with the corresponding phase shifts to each element based on the angle of arrival, which does not completely describe the response of an N-element array. $\endgroup$ – Envidia Nov 18 '19 at 21:33
2
$\begingroup$

You are really close! Change your signal and steering vectors to be complex. Specifically for the steering vectors, these coefficients are meant to act as phase shifts. Using a real sinusoid will introduce a phase shift term in the opposite angle direction, which you don't want. Doing this alone you will see an improvement in your pseudospectrum.

In regards to the extra non-zero eigenvalue:

Your choice of eigenvectors that span the noise subspace are based on your definition of their associated eigenvalues that are near zero. This is subjective of course but can be appropriate depending on the circumstances. Be careful when going this route.

The more general solution is to sort the eigenvectors based on their eigenvalues in decreasing order. MATLAB's svd() will do this for you. Then choose the M+1 to N eigenvectors to span the noise subspace where M is the number of expected signals and N is the number of antenna elements. In your case, M = 1 and N = 4 so you will pick the eigenvectors at columns 2, 3, and 4.

Some quick results:

After changing your signals and steering vectors to be complex and using svd() as described above, a run gives the following eigenvalues:

$$ \ \begin{bmatrix} 4.000 & 0 & 0 & 0 \\ 0 & 7.141*10^{-14} & 0 & 0 \\ 0 & 0 & 6.0258*10^{-14} & 0 \\ 0 & 0 & 0 & 2.5527*10^{-14} \\ \end{bmatrix} $$

These will produce a nice pseudospectrum after choosing the noise eigenvectors. Given that you were close, below is an edited version of your code:

% Assume 0 degrees due north, with positive clockwise orientation
% Positive phase difference indicates that antennna 1 recieves the waveform
% before the respective antenna (2,3, or 4)
% 
% Fix antenna 1 at the origin. From the x-axis, rotate antenna 2 60 degrees,
% rotate antenna 3 90 degrees, and rotate antenna 4 120 degrees. The
% optimal antenna spacing is still to be determined.
% 
%                         | 0 deg         --> +
%                         |
%                         |
%                         |
%             ---------- A1 ----------
%                   -60 ' | `  60
%                      '  |  `
%                     '   A3  `
%                  A4'         `A2

clearvars;
close all;

%%

num_antennas = 40;
num_sig_sources = 1;
RF_freq = 600e6;        % freq of incoming signal
lambda = 3e8/RF_freq;   % wavelength of incoming signal
fs = 100e6;             % sampling frequency
samples = 2000;          % number of snapshots
cos_freq = 1e6;         % freq to downconvert to

ant_spacings = [.25, .98, .91]; % antenna spacing

% generate_lookup_table just calculates a set of theoretical phase
% differences corresponding to some random AoA for the purpose of signal
% generation.
n = randi([0,360],1);
[AoAs, PDs1to2, PDs1to3, PDs1to4] = generate_lookup_table(-90, 90, .5, lambda, ant_spacings(1), ant_spacings(2), ant_spacings(3));
AoA = AoAs(n);

% generate signal
t = linspace(0,samples/fs,samples);
% recived_signal_ant_1 = cos(2*pi.*t.*cos_freq);
% recived_signal_ant_2 = cos(2*pi.*t.*cos_freq + PDs1to2(n));
% recived_signal_ant_3 = cos(2*pi.*t.*cos_freq + PDs1to3(n));
% recived_signal_ant_4 = cos(2*pi.*t.*cos_freq + PDs1to4(n));

% Change signals to be complex
recived_signal_ant_1 = exp(1i*2*pi.*t.*cos_freq);
recived_signal_ant_2 = exp(1i*(2*pi.*t.*cos_freq + PDs1to2(n)));
recived_signal_ant_3 = exp(1i*(2*pi.*t.*cos_freq + PDs1to3(n)));
recived_signal_ant_4 = exp(1i*(2*pi.*t.*cos_freq + PDs1to4(n)));


% build sample correlation matrix
cov_mat = zeros(4,4);
for i = 1:samples
    X = [recived_signal_ant_1(i);
         recived_signal_ant_2(i); 
         recived_signal_ant_3(i);
         recived_signal_ant_4(i)];
    cov_mat = cov_mat + X*X';
end

cov_mat = cov_mat*(1/samples);

% grab eigenvalues and eigenvectors
% [eig_vecs, eig_vals] = eig(cov_mat);

% generate matrix w/ noise eigenvector columns. An eigenvector corresponds
% to noise if the associated eigenvalue is close to 0. 
% col_num = 1;
% for i = 1:4
%     if eig_vals(i,i) < 1e-4
%         noise_vecs(:,col_num) = eig_vecs(:,i);
%         col_num = col_num + 1;
%     end
% end

% Use svd() instead
[eig_vecs, eig_vals] = svd(cov_mat);
noise_vecs = eig_vecs(:,num_sig_sources+1:end);

% walk t
theta = linspace(-90,90,200);
spectrum_func = zeros(1,100);
for i = 1:length(theta)

    % steering vector - only taking real part as of now (Re(e^i*phasediff) = cos(phase_diff))
    steering_vec = [1;
%         cos((-2*pi*ant_spacings(1)/lambda)*sin(theta(i)*(pi/180) - pi/3));
%         cos((-2*pi*ant_spacings(2)/lambda)*sin(theta(i)*(pi/180) - pi/2));
%         cos((-2*pi*ant_spacings(3)/lambda)*sin(theta(i)*(pi/180) - 2*pi/3))];

    % Make the steering vectors complex
    exp((1i*(-2*pi*ant_spacings(1)/lambda)*sin(theta(i)*(pi/180) - pi/3)));
        exp((1i*(-2*pi*ant_spacings(2)/lambda)*sin(theta(i)*(pi/180) - pi/2)));
        exp((1i*(-2*pi*ant_spacings(3)/lambda)*sin(theta(i)*(pi/180) - 2*pi/3)))];
    spectrum_func(i) = 1/((steering_vec'*noise_vecs)*noise_vecs'*steering_vec);
end

% grab theta corresponding to highest peak, round to nearest half integer
[mx,argmax] = max(spectrum_func);
estimated_AoA = floor(2*theta(argmax)+.5)/2;
fprintf("Estimated AoA: %.1f \t Hard-coded AoA: %.1f.\n", estimated_AoA, AoA)

plot(theta,20*log10(abs((spectrum_func))))
title('MUSIC Psuedospectrum')
xlabel('Angle of Arrival')

```
| improve this answer | |
$\endgroup$
  • $\begingroup$ Ah, I thought this might have had something to do with it. Is there any way you could elaborate on the 'using a real sinusoid will introduce a phase shift term in the opposite angle direction' aspect? $\endgroup$ – OGBerglemir Nov 19 '19 at 20:00
  • $\begingroup$ @OGBerglemir When you multiply two complex exponentials, you can add the general phase terms together. If you recall, a real sinusoid is composed of two complex exponentials, which rotate in different directions. Multiplying by a real sinusoid uses both components, so you will shift a signal's phase using both, and not one as intended. $\endgroup$ – Envidia Nov 19 '19 at 20:07
  • $\begingroup$ Ah, understood. $\endgroup$ – OGBerglemir Nov 19 '19 at 21:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.