Understanding the code, 2D-MUSIC Estimation DoA

100 views (last 30 days)
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

Accepted Answer

Tanmay Das
Tanmay Das on 14 Sep 2021
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

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!