ESPRIT算法及估计性能误差分析

ESPRIT算法程序和算法的估计性能误差分析，DOA估计

%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;%信元数
sensor_number=8;%原阵元数
m=7;%子阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=pi/4;%信号频率
l=2*pi*3e8/w;%信号波长 
d=0.5*l;%阵元间距

source_doa=50;%信号的入射角度

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]）;%仿真信号

%x=awgn（ssnr）;
x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪声后的阵列接收信号

x1=x（1:m:）;%子阵1接受的数据矢量
x2=x（2:（m+1）:）;%子阵2接受的数据矢量

%对两个子阵的模型进行合并
X=[x1;x2];
R=X*X‘/snapshot_number;
%对R进行奇异值分解
[USV]=svd（R）;
Us=U（:1:source_number）;
disp（Us）;
Us1=Us（1:m:）;
Us2=Us（（m+1）:2*m:）;

%按照公式得到旋转不变矩阵M
M=pinv（Us1）*Us2;
disp（‘M‘）;
disp（M）;
%对得到的旋转不变矩阵进行特征分解
[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）;

%求解均方根误差和标准偏差
%E_doa=sum（aaa（1:））/300;%做300次试验的均值
E_doa=mean（aaa）;
disp（E_doa）;

RMSE_doa=sqrt（sum（（aaa（1:）-source_doa）.^2）/300）;%做300次试验的均方根误差
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]）;%仿真信号
%x=awgn（ssnr）;
x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪声后的阵列接收信号

x1=x（1:m:）;%子阵1接受的数据矢量
x2=x（2:（m+1）:）;%子阵2接受的数据矢量

%对两个子阵的模型进行合并
X=[x1;x2];
R=X*X‘/snapshot_number;
%对R进行奇异值分解
[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”进行特征分解，得到矩阵E
[ESaVa]=svd（Us12‘*Us12）;
disp（‘E‘）;
disp（E）;
disp（Sa）;
%将E分解为四个小矩阵
E11=E（11）;
E12=E（12）;
E21=E（21）;
E22=E（22）;
%按照公式得到旋转不变矩阵M
M=-（E12*（inv（E22）））;
disp（‘M‘）;
disp（M）;
%对得到的旋转不变矩阵进行特征分解
[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）;

%求解均方根误差和标准偏差
%E_doa=sum（aaa（1:））/300;%做300次试验的均值
E_doa=mean（aaa）;
disp（E_doa）;

RMSE_doa=sqrt（sum（（aaa（1:）-source_doa）.^2）/300）;%做300次试验的均方根误差
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]）;%仿真信号
% %x=awgn（ssnr）;
% x=A*s+（1/sqrt（2））*（randn（sensor_numberN_x）+j*randn（sensor_numberN_x））;%加了高斯白噪声后的阵列接收信号
% 
% R=x*x‘/snapshot_number;
% disp（R）;
% %对接收到的数据协

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           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算法估计性能分析.pdf