二维DOA估计中的ROOTMUSIC算法

%ROOT_MUSIC ALOGRITHM
%DOA ESTIMATION BY ROOT_MUSIC
clear all;
close all;
clc;
tic
N=100;
source_number=2;%信元数
sensor_number=8;%阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/6].‘;%信号频率
l=（（2*pi*3e8）/w（1）+（2*pi*3e8）/w（2））/2;%信号波长  
d=0.5*l;%阵元间距
snr=0;%信噪比

theta1=50;theta2=60;%两个信号的入射角度
fai1=60;fai2=50;%方位角

Ax=[exp（-j*（0:sensor_number-1）*d*2*pi*cos（fai1*pi/180）*sin（theta1*pi/180）/l）;exp（-j*（0:sensor_number-1）*d*2*pi*cos（fai2*pi/180）*sin（theta2*pi/180）/l）].‘;%阵列流型

s=10.^（（snr/2）/10）*exp（j*w*[0:N_x-1]）;%仿真信号
%x=awgn（ssnr）;
for n=1:N
x=Ax*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪声后的阵列接收信号

R=（x*x‘）/N_x;
%对协方差矩阵进行特征分解
[VD]=eig（R）;
Un=V（:1:sensor_number-source_number）;
Gn=Un*Un‘;

%找出多项式的系数，并按阶数从高至低排列
a = zeros（2*sensor_number-11）‘;%1*15的0矩阵
for i=-（sensor_number-1）:（sensor_number-1）
    a（i+sensor_number） = sum（ diag（Gni） ）;%Gn的第i个对角元的和，放进行向量a中。
end
%使用ROOTS函数求出多项式的根
a1=roots（a）;%求多项式的根
%除掉大于1的增根（方程式在此条件下变形的）
a2=a1（abs（a1）<1）;
%disp（‘a1‘）;
%disp（a1）;
%disp（‘a2‘）;
%disp（a2）;
%挑选出最接近单位圆的N个根
[lamdaI]=sort（abs（abs（a2）-1））;
%disp（‘lamda‘）;
%disp（lamda）;
%disp（‘I‘）;
%disp（I）;
f1=a2（I（1:source_number））;
angle1=angle（f1）/pi;
%计算信号到达方向角
%source_doa=[asin（angle（f（1））/pi）*180/pi asin（angle（f（2））/pi）*180/pi];
%source_doa=sort（source_doa）;

%disp（‘source_doa‘）;
%disp（source_doa）;
Ay=[exp（-j*（0:sensor_number-1）*d*2*pi*sin（fai1*pi/180）*sin（theta1*pi/180）/l）;exp（-j*（0:sensor_number-1）*d*2*pi*sin（fai2*pi/180）*sin（theta2*pi/180）/l）].‘;
s2=10.^（（snr/2）/10）*exp（j*w*