%LS_ESPRIT?ALOGRITHM
%DOA?ESTIMATION?BY?LS_ESPRIT?ALOGRITHM
clear?all;
%close?all;
clc;

bbb=zeros（111）;
for?kk=1:11
snr=[-10?-8?-6?-4?-2?0?2?4?6?8?10];%信噪比（dB）

aaa=zeros（1300）;
for?k=1:300

source_number=1;%信元數(shù)
sensor_number=8;%原陣元數(shù)
m=7;%子陣元數(shù)
N_x=1024;?%信號(hào)長度
snapshot_number=N_x;%快拍數(shù)
w=pi/4;%信號(hào)頻率
l=2*pi*3e8/w;%信號(hào)波長?
d=0.5*l;%陣元間距

source_doa=50;%信號(hào)的入射角度

A=[exp（-j*（0:sensor_number-1）*d*2*pi*sin（source_doa*pi/180）/l）].‘;
s=10.^（（snr（kk）/2）/10）*exp（j*w*[0:N_x-1]）;%仿真信號(hào)

%x=awgn（ssnr）;
x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪聲后的陣列接收信號(hào)

x1=x（1:m:）;%子陣1接受的數(shù)據(jù)矢量
x2=x（2:（m+1）:）;%子陣2接受的數(shù)據(jù)矢量

%對兩個(gè)子陣的模型進(jìn)行合并
X=[x1;x2];
R=X*X‘/snapshot_number;
%對R進(jìn)行奇異值分解
[USV]=svd（R）;
Us=U（:1:source_number）;
disp（Us）;
Us1=Us（1:m:）;
Us2=Us（（m+1）:2*m:）;

%按照公式得到旋轉(zhuǎn)不變矩陣M
M=pinv（Us1）*Us2;
disp（‘M‘）;
disp（M）;
%對得到的旋轉(zhuǎn)不變矩陣進(jìn)行特征分解
[VmDm]=eig（M）;
disp（Dm）;
estimated_doa=-asin（angle（Dm（11））/pi）*180/pi;
disp（estimated_doa）;

aaa（:k）=estimated_doa;

disp（‘estimated_doa‘）;
disp（estimated_doa）;
end
disp（aaa）;

%求解均方根誤差和標(biāo)準(zhǔn)偏差
%E_doa=sum（aaa（1:））/300;%做300次試驗(yàn)的均值
E_doa=mean（aaa）;
disp（E_doa）;

RMSE_doa=sqrt（sum（（aaa（1:）-source_doa）.^2）/300）;%做300次試驗(yàn)的均方根誤差
disp（‘RMSE_doa‘）;
disp（RMSE_doa）;

bbb（:kk）=RMSE_doa;
end
disp（bbb）;
plot（snrbbb（1:）‘k*-‘）;
save?LS_ESPRIT_snr_rmse.mat;
%TLS_ESPRIT
bbb=zeros（111）;
for?kk=1:11
snr=[-10?-8?-6?-4?-2?0?2?4?6?8?10];%信噪比（dB）

aaa=zeros（1300）;
for?k=1:300
????
s=sqrt（10.^（snr（kk）/10））*exp（j*w*[0:N_x-1]）;%仿真信號(hào)
%x=awgn（ssnr）;
x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪聲后的陣列接收信號(hào)

x1=x（1:m:）;%子陣1接受的數(shù)據(jù)矢量
x2=x（2:（m+1）:）;%子陣2接受的數(shù)據(jù)矢量

%對兩個(gè)子陣的模型進(jìn)行合并
X=[x1;x2];
R=X*X‘/snapshot_number;
%對R進(jìn)行奇異值分解
[USV]=svd（R）;
Us=U（:1:source_number）;
disp（Us）;
Us1=Us（1:m:）;
Us2=Us（（m+1）:2*m:）;
Us12=[Us1?Us2];
%形成矩陣Us12
Us12=[Us1Us2];
%對“Us12‘*Us12”進(jìn)行特征分解，得到矩陣E
[ESaVa]=svd（Us12‘*Us12）;
disp（‘E‘）;
disp（E）;
disp（Sa）;
%將E分解為四個(gè)小矩陣
E11=E（11）;
E12=E（12）;
E21=E（21）;
E22=E（22）;
%按照公式得到旋轉(zhuǎn)不變矩陣M
M=-（E12*（inv（E22）））;
disp（‘M‘）;
disp（M）;
%對得到的旋轉(zhuǎn)不變矩陣進(jìn)行特征分解
[VmDm]=eig（M）;
disp（Dm）;
estimated_doa=-asin（angle（Dm（11））/pi）*180/pi;
disp（estimated_doa）;

aaa（:k）=estimated_doa;

disp（‘estimated_doa‘）;
disp（estimated_doa）;
end
disp（aaa）;

%求解均方根誤差和標(biāo)準(zhǔn)偏差
%E_doa=sum（aaa（1:））/300;%做300次試驗(yàn)的均值
E_doa=mean（aaa）;
disp（E_doa）;

RMSE_doa=sqrt（sum（（aaa（1:）-source_doa）.^2）/300）;%做300次試驗(yàn)的均方根誤差
disp（‘RMSE_doa‘）;
disp（RMSE_doa）;

bbb（:kk）=RMSE_doa;
end
disp（bbb）;
hold?on
plot（snrbbb（1:）‘rs-‘）;
save?TLS_ESPRIT_snr_rmse.mat;
%?%TAM
%?bbb=zeros（111）;
%?for?kk=1:11
%?snr=[-10?-8?-6?-4?-2?0?2?4?6?8?10];%信噪比（dB）
%?
%?aaa=zeros（1300）;
%?for?k=1:300
%?????
%?s=sqrt（10.^（snr（kk）/10））*exp（j*w*[0:N_x-1]）;%仿真信號(hào)
%?%x=awgn（ssnr）;
%?x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪聲后的陣列接收信號(hào)
%?
%?R=x*x‘/snapshot_number;
%?disp（R）;
%?%對接收到的數(shù)據(jù)協(xié)

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----
?????目錄???????????0??2014-12-08?21:14??ESPRIT\
?????文件????????9502??2014-12-08?21:14??ESPRIT\esprit_rmse_snr.m
?????文件????????4715??2010-01-03?12:23??ESPRIT\esprit_rmse_zhenyuan.m
?????文件??????497940??2014-12-08?21:14??ESPRIT\LS_ESPRIT_snr_rmse.mat
?????文件?????1601958??2010-01-03?12:23??ESPRIT\RB_esprit_snr_rmse.mat
?????文件?????1346930??2010-01-03?12:23??ESPRIT\resprit_snr_rmse.mat
?????文件??????497063??2010-01-03?12:24??ESPRIT\TAM_snr_rmse.mat
?????文件??????498611??2010-01-03?12:24??ESPRIT\TLS_ESPRIT_snr_rmse.mat
?????文件???????12376??2010-01-03?12:24??ESPRIT\untitled.fig
?????文件?????1887502??2014-12-08?21:14??ESPRIT算法估計(jì)性能分析.pdf