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