esprit_1.m

一个自己写的ESPRIT的程序

clc;
clear all;
N=8;%阵元数目
N1=N-1;%子阵阵元数
p=5; %信源数 
M=1024; %快拍数
theta=[10 20 30 40 50];%波束的到达方向
lamda=0.50;    %波长
d=lamda/2;%阵元间距
SNR=25;
z=[0,0,0,0,0,0,0;1,0,0,0,0,0,0;0,1,0,0,0,0,0;0,0,1,0,0,0,0;0,0,0,1,0,0,0;0,0,0,0,1,0,0;0,0,0,0,0,1,0];
u(1:8)=[1.00,1.2,0.85,1.12,1.2,0.9,1.2,0.86];
f(1:8)=pi*[0,2.32,-5.28,8.91,-4.69,6.49,-2.11,1.19]/180;
tt=[u(1)*exp(j*f(1)),u(2)*exp(j*f(2)),u(3)*exp(j*f(3)),u(4)*exp(j*f(4)),u(5)*exp(j*f(5)),u(6)*exp(j*f(6)),u(7)*exp(j*f(7)),u(8)*exp(j*f(8))]
tt=diag(tt);
%%%%%%%%%%%%%  空间信号采样 %%%%%%%%%%%%
fs=1024; %采样频率
n=[0:N-1]';
A=exp(j*2*pi*n*sin(theta/180*pi)*d/lamda);%阵列流型
t=(0:M-1)/fs;
b=10*[1:p]';    %线性调频系数
S=exp(j*2*pi*(1/2*b*t.^2+randn(p,1)*ones(1,M))); % 为了便于仿真，不考虑载频
sigma_2=10^(SNR/10);
X=sqrt(sigma_2)*A*S+(0.707*randn(N,M)+j*0.707*randn(N,M)); % 加白噪声 SNR
XX=sqrt(sigma_2)*tt*A*S+(0.707*randn(N,M)+j*0.707*randn(N,M)); % 加白噪声 SNR
X1=X(1:N-1,:);%子阵1接收的数据
X2=X(2:N,:);%子阵2接收的数据
X11=XX(1:N-1,:);%子阵1接收的数据
X22=XX(2:N,:);%子阵2接收的数据
%%%%%%%%%%%%%%%%    ESPRIT 算法       %%%%%%%%%%%%%%%%%%%%
R11=X1*X1'/M;      %估计子阵1的相关矩阵
R12=X1*X2'/M;       %估计两子阵的互相关矩阵
D1= eig(R11);        %对相关矩阵R11作特征值分解
u2=mean(D1(1:N1-p));  %估计噪声的方差
C11=R11-u2*eye(N1,N1); %计算子阵1 的协方差矩阵
C12=R12-u2*z; %计算子阵1 的协方差矩阵
D=eig(C11,C12);         % 对C11和C12作广义特征值分解
%theta_1 = asin((angle(D)/pi))*180/pi %计算来波方向
D_abs=abs(D);
D_big=zeros(1,p);
%%%%%%%%%%%%%%%%    误差ESPRIT 算法       %%%%%%%%%%%%%%%%%%%%
R111=X11*X11'/M;      %估计子阵1的相关矩阵
R122=X11*X22'/M;       %估计两子阵的互相关矩阵
D11= eig(R111);        %对相关矩阵R11作特征值分解
u22=mean(D11(1:N1-p)),  %估计噪声的方差
C111=R111-u22*eye(N1,N1); %计算子阵1 的协方差矩阵
C122=R122-u22*z; %计算子阵1 的协方差矩阵
DD=eig(C111,C122);         % 对C11和C12作广义特征值分解
D1_abs=abs(DD);
D1_big=zeros(1,p);
for i=1:p     %找到P个大特征值
    [m1,k]=max(D);
    D_big(i)=D(k);
    D(k)=0;
    
%     D1_big(i)=DD(k);

 end
theta_ = asin((angle(D_big)/pi))*180/pi %计算来波方向
theta_1= asin((angle(D1_big)/pi))*180/pi %计算来波方向