MUSIC algorithm for range-azimuth FMCW data processing

I'm using the following 1D music algorithm to estimate targets position (only theta angle).

It works pretty good with my 1D antenna array. I need help to trasform it in the 2D case (azimuth and range) in near field scenario.

Thanks a lot in advance,

Ciao

Luca

% ======= (1) TRANSMITTED SIGNALS ======= %

% Signal source directions
az = [35;39;127]; % Azimuths
el = zeros(size(az)); % Simple example: assume elevations zero
M = length(az); % Number of sources

% Transmitted signals
L = 200; % Number of data snapshots recorded by receiver
m = randn(M,L); % Example: normally distributed random signals

% ========= (2) RECEIVED SIGNAL ========= %

% Wavenumber vectors (in units of wavelength/2)
k = pi*[cosd(az).*cosd(el), sind(az).*cosd(el), sind(el)].';

% Array geometry [rx,ry,rz]
N = 10; % Number of antennas
r = [(-(N-1)/2:(N-1)/2).',zeros(N,2)]; % Assume uniform linear array

% Matrix of array response vectors
A = exp(-1j*r*k);

% Additive noise
sigma2 = 0.01; % Noise variance
n = sqrt(sigma2)*(randn(N,L) + 1j*randn(N,L))/sqrt(2);

% Received signal
x = A*m + n;

% ========= (3) MUSIC ALGORITHM ========= %

% Sample covariance matrix
Rxx = x*x'/L;
d = 0.5;% Distance between elements in Wavelenght
[Q ,D]=eig(Rxx); %Compute eigendecomposition of covariance matrix
[D,I]=sort(diag(D),1,'descend'); %Find r largest eigenvalues
Q=Q (:,I); %Sort the eigenvectors to put signal eigenvectors first
Qs=Q (:,1:M); %Get the signal eigenvectors
Qn=Q(:,M+1:N); %Get the noise eigenvectors
% MUSIC algorithm
% Define angles at which MUSIC “spectrum” will be computed
angles=(-90:1:90);

for k=1:length(angles)
a1(:,k)=exp(-1i*2*pi*(d*(0:N-1)'*sin([angles(k).']*pi/180)));
end
for k=1:length(angles)
%Compute MUSIC “spectrum”
music_spectrum(k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k));

end

plot(angles,abs(music_spectrum))
grid on
title('MUSIC Spectrum')
xlabel('Angle in degrees')

[Shirleyuue Jiang]

Hi, did you succeed in transform 1D to 2D? I have some problems with 2D,can you show me your whole code?My email: [e-mail].

Answer 1

[ruoyan pan]

Hi, I think I have the same problem as yours, did you succeed to solve it？If you did, could you please help me? My email：[e-mail]

Answer 2

[Muhammad]

clc; clear; close all;

% ======= (1) TARGET DEFINITION ======= %

az = [35, 39, 127]; % Azimuth angles (degrees)

range = [3, 3.5, 4]; % Range in wavelengths

el = zeros(size(az)); % Elevation fixed (0 degrees)

M = length(az); % Number of sources

% Cartesian coordinates of targets (in wavelength units)

target_pos = range(:) .* [cosd(el(:)) .* cosd(az(:)), ...

cosd(el(:)) .* sind(az(:)), ...

sind(el(:))]; % [M x 3]

% ======= (2) TRANSMITTED SIGNALS ======= %

L = 200; % Number of snapshots

m = randn(M, L); % Random signals

% ======= (3) ARRAY GEOMETRY (2D URA) ======= %

d = 0.5; % Inter-element spacing in wavelengths

Nx = 5; Ny = 5; % URA dimensions

N = Nx * Ny; % Total number of antennas

[x_grid, y_grid] = meshgrid(0:Nx-1, 0:Ny-1);

r = [x_grid(:), y_grid(:), zeros(N,1)] * d; % [N x 3] antenna positions

% ======= (4) RECEIVED SIGNAL (NEAR-FIELD) ======= %

A = zeros(N, M);

for mIdx = 1:M

dist = sqrt(sum((r - target_pos(mIdx,:)).^2, 2)); % distance from each antenna to target

A(:, mIdx) = exp(-1j*2*pi*dist); % spherical wavefront

end

% Additive complex Gaussian noise

sigma2 = 0.01;

n = sqrt(sigma2) * (randn(N, L) + 1j * randn(N, L)) / sqrt(2);

% Received signal

x = A * m + n; % [N x L]

% ======= (5) MUSIC ALGORITHM ======= %

Rxx = x * x' / L; % Sample covariance matrix

[Q, D] = eig(Rxx);

[D_sorted, idx] = sort(diag(D), 'descend');

Q = Q(:, idx);

Qs = Q(:, 1:M); % Signal subspace

Qn = Q(:, M+1:end); % Noise subspace

% ======= (6) MUSIC SPECTRUM (AZIMUTH & RANGE) ======= %

az_grid = -90:1:90; % Azimuth angle search grid

range_grid = 2:0.05:5; % Range search grid (in wavelengths)

music_spectrum = zeros(length(range_grid), length(az_grid));

for i = 1:length(range_grid)

for j = 1:length(az_grid)

test_pos = range_grid(i) * [cosd(az_grid(j)), sind(az_grid(j)), 0];

dist = sqrt(sum((r - test_pos).^2, 2));

a = exp(-1j * 2 * pi * dist);

music_spectrum(i,j) = 1 / (a' * (Qn * Qn') * a);

end

end

% ======= (7) PLOT MUSIC SPECTRUM ======= %

figure;

imagesc(az_grid, range_grid, 10*log10(abs(music_spectrum)));

xlabel('Azimuth (°)');

ylabel('Range (wavelengths)');

title('2D MUSIC Spectrum (Azimuth vs. Range)');

colorbar;

axis xy;

hope this will help.