Understanding the code, 2D-MUSIC Estimation DoA

43 次查看（过去 30 天）
Giovanni Latella
Giovanni Latella 2021-9-10
评论： 某 刘 2023-10-16
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

请先登录，再回答此问题。

采纳的回答

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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Beamforming and Direction of Arrival Estimation 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by