Understanding the code, 2D-MUSIC Estimation DoA

Hi all, I have developed the following code for estimating the direction of arrival of a signal on a sensor array (URA) with the MUSIC algorithm. The problem is that I couldn't estimate the two angles, azimuth and zenith, and I copied the max_matrix function found online. The code works perfectly, but I cannot understand the meaning of this function, despite being commented. Please help me. Thank you.

clear , clf
M = 10; N = 5;   % URA
d1 = 0.5; d2 = 0.5;  % Distance between antennas


f1 = 450.0;
phd = 200; 
thd = 60;
fc = 150e4;        
c = physconst('LightSpeed');
lamb = c/fc;
k = (2*pi)/lamb;

x = exp(1i*(k-(2*pi*f1)));    % Wave plane
noise = 0.001/sqrt(2)*(randn(size(x)) + 1i*randn(size(x)));

pi2 = 2*pi; j2pi = 1j*pi2;    
a_URA = @(ph,th,m,n)exp(j2pi*sin(th)*(d1*m*cos(ph)+d2*n*sin(ph)));  % Steering vector
mm = (0:M-1).'; nn = 0:N-1;   % Row and column
ph = deg2rad(phd);  % Azimuth in radians
th = deg2rad(thd);  % % Zenith in radians
yy = x*a_URA(ph,th,mm,nn) + noise;  % MxN matrix of signals received by each antenna
y_ = yy(:);
R = y_*y_';    % Correlation matrix


[phh,thh,P,phph,thth] = DoA2_MUSIC_URA(R,M,N,d1,d2)
PdB = log10(abs(P));    % Spatial spectrum in dB
mesh(rad2deg(phph),rad2deg(thth),PdB); hold on
plot3(rad2deg(phh),rad2deg(thh),max(max(PdB)),'^')   % Pick at DoA
xlabel('Azimuth Angle')
ylabel('Zenith Angle')


function [phh,thh,P,phph,thth] = DoA2_MUSIC_URA (R,M,N,d1,d2,phph,thth)    % MUSIC 2D-DoA in URA MxN 
mm = (0:M-1).'; nn = 0:N-1; j2pi = 1j*2*pi;
a_URA = @(ph,th,m,n)exp(j2pi*sin(th)*(d1*m*cos(ph)+d2*n*sin(ph)));
thth = (0:90)*(pi/180);     
phph = (0:360)*(pi/180);    % thth/phph = range of Azimuth/Zenith (phi/theta) [rad]
[V,Lam] = eig(R);   % Eigendecomposition
[lambdas,idx] = sort(abs(diag(Lam)));
[~,Nn] = max(diff(log10(lambdas+1e-3)));  % Nn
Vn = V(:,idx(1:Nn));     % MNxNn matrix consisting of presumed noise eigenvectors
VVn = Vn*Vn';
for m = 1:length(phph)
    for n = 1:length(thth)
        ph = phph(m); th = thth(n); 
        a = a_URA(ph,th,mm,nn);   
        a = a(:);      % Steering vector
        P(n,m) = 1/(a'*VVn*a);   % Spatial spectrum for MUSIC
    end
end
[nmax,mmax] = max_matrix(P);       
phh = phph(mmax); thh = thth(nmax);   % Estimation in radians of the Azimuth / Zenith angles
end

function [m,n,amax] = max_matrix(A)
% m/n = Row/Column number of maximum entry, amax, of matrix A
[amax,idx] = max(A(:)); M = size(A,1);
m = rem(idx,M); n = (idx-m)/M+1;     
end

[某 刘]

What does f1 mean？

Answer 1

[Tanmay Das]

Hi,

I am assuming that you are facing trouble in understanding the code of max_matrix function. A(:) reshapes all elements of A into a single column vector. This has no effect if A is already a column vector. max() function returns the maximum value(along with its index) of a column/row vector. if the vector has complex values, it return the element that has the largest magnitude. So, amax consists of the element with largest magnitude in column vector A(:) and idx consists the index of that particular element.

Now, rem(a,b) returns the remainder after division of a by b. The operation is done because now it is needed to convert the idx( i.e., with respect to column vector A(:)) into row(m) and column(n) indices of matrix A. Recall that A(:) was formed by stacking each column of matrix A below another. For example:

A = [1 2 3;... % 3 x 3 matrix
    4 5 6;...
    7 8 9];
B = A(:);      % 9 x 1 column vector [1;4;7;2;5;8;3;6;9]
% Suppose we need to find out the row and column indices of 
% 5th element of B i.e. 5
idx = 5;
M = size(A,1); % number of columns in A i.e. 3

m = rem(idx,M); % we get m = rem(5,3) => m = 2
n = (idx-m)/M+1; % we get n = [(5 - 2)/3] + 1 => n = 2