esprint.m

ULA1一个天线阵的方向图,可以用来估计信号的来波方向,face_direction_pic面阵的方向图 music,MN-music,capon,ESPRINT为四种进测量来波方向的算法 该程序包主

clear;close all;clc;
%=========参数初始化=========%
num=10;%阵元数
signal_number=3;%信号数
a1=20;%信号1方向
a2=80;%信号2方向
a3=40;%信号3方向
f=10e9;%载波频率
f1=20;%信号1载波频率
f2=10;%信号2载波频率
f3=5;%信号3载波频率
c=3e8;%光速
lmda=c/f;%波长
d=lmda/2;%阵元间距
% thita=0:1:180;
SNR1=20;%信噪比
SNR2=20;%信噪比
SNR3=20;%信噪比
fs=500;
Ts=[0:1/fs:1];
snapshot=length(Ts);
s1=10^(SNR1/20)*sin(2*pi*f1*Ts);%信号1包络
s2=10^(SNR2/20)*sin(2*pi*f2*Ts);%信号2包络
s3=10^(SNR3/20)*sin(2*pi*f3*Ts);%信号2包络
%子阵1
athita1=exp(j*2*pi*d/lmda*cos(a1*pi/180)*[0:num-2]).';%信号导向矢量
athita2=exp(j*2*pi*d/lmda*cos(a2*pi/180)*[0:num-2]).';%信号导向矢量
athita3=exp(j*2*pi*d/lmda*cos(a3*pi/180)*[0:num-2]).';%信号导向矢量

x=athita1*s1+athita2*s2+athita3*s3+sqrt(0.5)*(randn(num-1,snapshot)+j*randn(num-1,snapshot));%输入
% x=athita1*s1+athita2*s2;
%子阵2
athita1=exp(j*2*pi*d/lmda*cos(a1*pi/180)*[1:num-1]).';%信号导向矢量
athita2=exp(j*2*pi*d/lmda*cos(a2*pi/180)*[1:num-1]).';%信号导向矢量
athita3=exp(j*2*pi*d/lmda*cos(a3*pi/180)*[1:num-1]).';%信号导向矢量
y=athita1*s1+athita2*s2+athita3*s3+sqrt(0.5)*(randn(num-1,snapshot)+j*randn(num-1,snapshot));%输入
% y=athita1*s1+athita2*s2;
%==============ESPRINT==============%
Rxx=x*x'/snapshot;
Rxy=x*y'/snapshot;
[V,D]=eig(Rxx);
s=0;
for i=1:num-1-signal_number;
    s=sum(D(:,i))+s;
end
thigama=s/(num-1-signal_number);
Cxx=Rxx-thigama*eye(size(Rxx));
Cxy=Rxy-thigama*diag([ones(1,num-1-1)],-1);  %diag([1 1 1 1 1 1 1 1],-1)
[p,q]=eig(Cxx,Cxy);
%[p,q]=eig(Rxx,Rxy);
for i=1:signal_number;
    thita(i)=abs(acos(-j*(log(q(i,i)))*lmda/(2*pi*d))*180/pi);
end
thita