clear?all
close?all
clc
derad=pi/180;%角度轉(zhuǎn)弧度
radeg=180/pi;
twpi=2*pi;
kelm=8;%陣元數(shù)
dd=0.5;%陣元間距
d=0:dd:（kelm-1）*dd;
iwave=3;%信源數(shù)
theta=[10?20?30];%波達方向
snr=10;%信噪比
n=500;%快拍數(shù)
A=exp（-j*twpi*d.‘*sin（theta*derad））;%方向矢量8*3
S=randn（iwaven）;%3*正態(tài)隨機分布矩陣3*500
X=A*S;%接收信號
X1=awgn（Xsnr‘measured‘）;%添加噪聲
Rxx=X1*X1‘/n;%計算協(xié)方差矩陣
InvS=inv（Rxx）;
[EVD]=eig（Rxx）;%D為特征值，EV為特征向量
EVA=diag（D）‘;%取特征矩陣的對角線上元素并轉(zhuǎn)置，結(jié)果為一個行向量
[EVAI]=sort（EVA）;%
F=EV;
EV=EV（:I）;%對應特征向量的位置
f=EV;
EVA=fliplr（EVA）;%左右翻轉(zhuǎn)，改為從左至右元素從大到小排列，即特征值從大到小排列
EV=fliplr（EV）;%對應特征向量也進行翻轉(zhuǎn)，特征值對應的特征向量
%構造music函數(shù)
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）;