mucis算法实现DOA估计

clear all
close all
clc
derad=pi/180;%角度转弧度
radeg=180/pi;
twpi=2*pi;
kelm=8;%阵元数
dd=0.5;%阵元间距
d=0:dd:（kelm-1）*dd;
iwave=3;%信源数
theta=[10 20 30];%波达方向
snr=10;%信噪比
n=500;%快拍数
A=exp（-j*twpi*d.‘*sin（theta*derad））;%方向矢量8*3
S=randn（iwaven）;%3*正态随机分布矩阵3*500
X=A*S;%接收信号
X1=awgn（Xsnr‘measured‘）;%添加噪声
Rxx=X1*X1‘/n;%计算协方差矩阵
InvS=inv（Rxx）;
[EVD]=eig（Rxx）;%D为特征值，EV为特征向量
EVA=diag（D）‘;%取特征矩阵的对角线上元素并转置，结果为一个行向量
[EVAI]=sort（EVA）;%
F=EV;
EV=EV（:I）;%对应特征向量的位置
f=EV;
EVA=fliplr（EVA）;%左右翻转，改为从左至右元素从大到小排列，即特征值从大到小排列
EV=fliplr（EV）;%对应特征向量也进行翻转，特征值对应的特征向量
%构造music函数
for i=1:361
    angle（i）=（i-181）/2;
    theta_2=derad*angle（i）;
    a=exp（-j*twpi*d*sin（theta_2））.‘;
    N=iwave;
    En=EV（:N+1:kelm）;
    SP（i）=（a‘*a）/（a‘*En*En‘*a）;
end
%绘图
sp=SP;
SP=abs（SP）;