Has anyone done DOA estimation for L-shaped arrays? there was a problem with the Angle matching part

%二维L型DOA估计 2-D root-music算法
clear all
close all
clc;

lamda = 1; %波长
derad = pi/180; %角度——>弧度
radeg = 180/pi; %弧度——>角度
twpi = 2*pi;
kelm = 9; %阵元数
theta = [25 50]; % 仰角 (0-90度)
fe = [10 30]; % 方位角 (0-90度)
iwave = length(theta); %信源数
SNR=-10:5:30;
%% 蒙特卡罗实验
num_MC=1000;%蒙特卡洛次数
for idx_SNR = 1 : length(SNR)
for idx_MC = 1 : num_MC
snr=SNR(idx_SNR);
%snr=10;
n = 400; %采样数
dd = 0.5*lamda; %阵元间距
dx = 0:dd:(kelm-1)*dd; %X轴阵元分布
dy = 0:dd:(kelm-1)*dd; %Y轴阵元分布
Ax = exp(-1j*twpi*dx.'*(sin(theta*derad).*cos(fe*derad))); %X轴上阵元对应的方向矩阵
Ay = exp(-1j*twpi*dy.'*(sin(theta*derad).*sin(fe*derad))); %Y轴上阵元对应的方向矩阵
A = [Ax;Ay];
S = randn(iwave,n); %信源信号
X0 = Ax*S; %X轴上接收信号
X = awgn(X0,snr,'measured'); %添加噪声
Y0 = Ay*S; %Y轴上接收信号
Y = awgn(Y0,snr,'measured'); %添加噪声
Rxx = X*X'; %计算协方差矩阵
[EVx,Dx] = eig(Rxx); %特征值分解
EVAx = diag(Dx)';
[EVAx,Ix] = sort(EVAx); %特征值从小到大排序
EVAx = fliplr(EVAx); %左右翻转，从大到小排序
EVx = fliplr(EVx(:,Ix)); %对应特征向量排序
Ryy = Y*Y'; %计算协方差矩阵
[EVy,Dy] = eig(Ryy); %特征值分解
EVAy = diag(Dy)';
[EVAy,Iy] = sort(EVAy); %特征值从小到大排序
EVAxy = fliplr(EVAy); %左右翻转，从大到小排序
EVy = fliplr(EVy(:,Iy)); %对应特征向量排序
Z = [X;Y]; %添加噪声
Rzz = Z*Z'/n; %自相关函数
[EVz,Dz] = eig(Rzz); %求矩阵的特征向量和特征值
[EVAz,Iz] = sort(diag(Dz).'); %特征值按升序排列
EVz = fliplr(EVz(:,Iz)); %左右翻转，特征值按降序排列
Unx = EVx(:,iwave+1:kelm); %提取噪声子空间
Uny = EVy(:,iwave+1:kelm); %提取噪声子空间
Un = EVz(:,iwave+1:end); %噪声子空间
%利用root-music进行求根
syms z1
pz1_1 = z1.^([0:kelm-1]');
pz1_2 =(z1^(-1)).^([0:kelm-1]);
fz1 = z1.^(kelm-1)*pz1_2*Unx*Unx'*pz1_1; %构造多项式
a1 = sym2poly(fz1); %符号多项式——>数值多项式
zx = roots(a1); %求根
rx = zx.';
%利用root-music进行求根
syms z2
pz2_1 = z2.^([0:kelm-1]');
pz2_2 =(z2^(-1)).^([0:kelm-1]);
fz2 = z2.^(kelm-1)*pz2_2*Uny*Uny'*pz2_1; %构造多项式
a2 = sym2poly(fz2); %符号多项式——>数值多项式
zy = roots(a2); %求根
ry = zy.';
rx1 = rx(abs(rx)<1); % 找出在单位圆里的根
ry1 = ry(abs(ry)<1);

[lamdax,Ix] = sort(abs(abs(rx1)-1)); % 挑选出最接近单位圆根
[lamday,Iy] = sort(abs(abs(ry1)-1));

gamax = -angle(rx1(Ix(1:iwave)))/pi; %把最接近的K个根
gamay = -angle(ry1(Iy(1:iwave)))/pi;

%利用接收信号与噪声子空间的正交性来匹配
for iu = 1:1:iwave
for iv = 1:1:iwave
gx = gamax(1,iu);
gy = gamay(1,iv);

Agx = exp(-1j*twpi*dx.'*( gx )); %X轴上阵元对应的方向矩阵
Agy = exp(-1j*twpi*dy.'*( gy )); %Y轴上阵元对应的方向矩阵
Ag = [Agx;Agy];

SP(iu,iv) = Ag'*Un*Un'*Ag;
end
end
SP1 = 1./SP; %变为求最大值
SP2 = reshape(SP1,1,iwave*iwave);

ddddd = 1;
Index_r = [];
for ix = 1:1:iwave
if ddddd <= iwave
[max_val,max_idx]= max(SP2);% 找到A中的最大值及其索引
col = floor((max_idx-1)/iwave) + 1; %列
row = mod((max_idx-1),iwave) + 1; %行
roco = [row;col];
Index_r = [Index_r,roco]; %正确的配对顺序
SP1(row,col) = 0;
SP2 = reshape(SP1,1,iwave*iwave);
ddddd = ddddd+1;
end
end
%利用配对顺序求出DOA估计的角度
for iu = 1:1:iwave
ax = Index_r(1,iu);
ay = Index_r(2,iu);
gx = gamax(1,ax);
gy = gamay(1,ay);
out1(1,iu) = asin( sqrt(gx.^2 + gy.^2) );
if gy/gx < 0
out2(1,iu) = atan( gy ./ gx ) + pi;
else
out2(1,iu) = atan( gy ./ gx );
end
end
DOA1=out1/derad;
DOA2=out2/derad;
DOA3(idx_MC,:)=DOA1;
DOA4(idx_MC,:)=DOA2;
errors(idx_MC)=sum((sort(theta(:)) - sort(DOA1(:))).^2+(sort(fe(:)) - sort(DOA2(:))).^2);
end
RMSE(idx_SNR) =(sqrt(sum(errors)/iwave/num_MC));
end
figure;
semilogy(SNR, RMSE, '-^', 'Linewidth', 1); % 绘制 CRB
ylim([0.01 10])