doa估计算法

clear 
i=sqrt（-1）; 
j=i;
m=8;%均匀线阵列数目
M=m-1;
p=1;%信号源数目 
lma=150;%入射信号波长
d=lma/2;%阵元间距
N=200;%快拍数
angle1=60.54;
th=[angle1];
degrad=pi/180; %求弧度所要公式
load sig3
%%.........................ESPRIT.............%%
tt=1:N; %信号的采样点数
S=[s1（tt）.‘]; 
bbb=zeros（110）;
for kk=1:10
snr=[-20 -15 -10 -5 0 5 10 15 20 25].‘;%信噪比（dB）
aaa=zeros（1100）;
for k=1:100
tmp=[0:m-1]‘;
A=exp（tmp*-i*pi*2*d*sin（th‘*degrad）/lma）;%方向矩阵中自然底数的方值
X=A*S;
X=X+awgn（Xsnr（kk））;
Rxx=X*X‘;%相关矩阵
x1=X（1:M:）;%两个子阵接收的数据
x2=X（2:M+1:）;
X1=[x1;x2];
Rxx1=X1*X1‘;%相关矩阵
[QDU]=svd（X1）;%求矩阵奇异值分解，得到Q为特征值对应的向量矩阵
V2=Q（:1:p）;%2M行，p列是取大特征值对应的特征向量
Vs1=V2（1:M:）;% 信号子空间中对应子阵列1的信号子空间     
Vs2=V2（M+1:2*M:）;%对应阵列2的信号子空间
W=Vs1\Vs2;%A\B = A1*B。矩阵乘法不满足交换律，因此需要有左右除法之分。矩阵求逆的命令是 inv ，即 A1 = inv（A）。
[IJ]=eig（W）;%对w矩阵你进行特征值分解，J是包含特征值的对角阵。信号的波达方向就在这些对角线元素上
J=（diag（J））;%取出矩阵J中对角线上元素
ESPRIT_doa=-asin（angle（J）/pi）*180/pi;%将特征值转换成角度  
aaa（:k）=ESPRIT_doa;
end
RMSE_doa=sqrt（sum（（aaa（1:）-th）.^2）/100）;%做100次试验的均方根误差
bbb（:kk）=RMSE_doa;
end
disp（‘ESPRIT_doa‘）;
disp（ESPRIT_doa）;  
plot（snrbbb（1:）‘k*-‘）;
save ESPRIT_snr_rmse.mat;
%%.................MUSIC...............%%
bbb=zeros（110）;
for kk=1:10
snr=[-20 -15 -10 -5 0 5 10 15 20 25].‘;%信噪比（dB）
aaa=zeros（1100）;
for k=1:100 
tmp=[0:m-1]‘;
A=exp（tmp*-i*pi*2*d*sin（th‘*degrad）/lma）;%方向矩阵中自然底数的方值
X=A*S;
X=X+awgn（Xsnr（kk））;
Rxx=X*X‘/N;%相关矩阵 m*m矩阵  
[V1D1]=eig（Rxx）; %求协方差矩阵的特征矢量及特征值，D是包含特征值的对角矩阵。V是其列是对应的右特征向量的矩阵。
Vn=V1（:1:m-p）; %取矩阵V中的噪声特征值对应的特征向量，后m-p个特征向量便是噪声特征值对应的，利用它构造噪声子空间。   
t