%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*