uesprit.m

用于信号处理中的检测DOA估计

function doa = uesprit( x,d,N )

X=x.';
J1=[eye(N-1) zeros(N-1,1)];   % (N-1)*N
J2=[zeros(N-1,1) eye(N-1)];

Qn=unitary(N).';
Qn_1=unitary(N-1).';

Y=Qn'*X;

[U S V]=svd([real(Y) imag(Y)]);
Es=U(:,1:d);

K1=real(Qn_1'*J2*Qn);
K2=imag(Qn_1'*J2*Qn);

Lu1=K1*Es;
Lu2=K2*Es;
fai=pinv(Lu1)*Lu2;
omeig=eig(fai);

mue=2*atan(omeig);
u=mue/pi;
if(u<0)
    u=-u;
elseif(u>0)
    u=1-u;
end

theta=asin(u)*180/pi;
doa=theta;


clear all;
clc;
%所需信号来自ang1度方向，而干扰信号来自ang2度方向，频率为freq=30MHZ,也即波长为wavelength=10米。阵元之间的间距为dist=wavelength/2米，
%假设阵列有30个阵元组成，记右边第一个阵元接收的信号为Sd，则第二个阵元接收的信号为Sd*EXP(-i*2*pi*d*sin(ang)/wavelength),
%以此类推下去。
ratio_d_and_w=0.5;
%N_singal is numbers of sampling singal;
N_signal=10000;% 样本数

ang1=(30*pi/180);%所需信号的方向
SNR=5;%信噪比
ASd=sqrt(10.^(SNR/10));  
ang2=(65*pi/180);%干扰信号的方向
INR=45;%干噪比
ASi=sqrt(10.^(INR/10));


N_array=30;%阵列数

Sd=ASd*(randn(1,N_signal)+i*randn(1,N_signal));%Sd为所需信号1*10000

Si=ASi*(randn(1,N_signal)+i*randn(1,N_signal));%Si为干扰信号1*10000

Ni=randn(N_signal,N_array)+i*randn(N_signal,N_array);%Ni内噪声10000*30

Desired_Array=zeros(N_signal,N_array);%10000*30
Interferential_Array=zeros(N_signal,N_array);%10000*30
%f=test(ang2,N_array,ratio_d_and_w);%30*1

for LL=1:N_signal
    Interferential_Array(LL,:)=Si(LL)*test(-ang2,N_array,ratio_d_and_w).';
    Desired_Array(LL,:)=Sd(LL)*test(-ang1,N_array,ratio_d_and_w).';
end
%a_a=Desired_Array(1,:)
X=zeros(N_signal,N_array);

%X=Ni;
X=Desired_Array + Interferential_Array +Ni;

Rx=cov(X);
%AA=size(Rx);
S=test(ang1,N_array,ratio_d_and_w);
Wopt=inv(Rx)*S/(S'*inv(Rx)*S);%最优权矢量

%从-90度到90度，间隔为0.1度；信号的方向
x_axis=-90:0.5:90;
ang=(x_axis)*pi/180;
Yy=zeros(1,length(ang));
for L=1:length(ang)
    Z2=test(ang(L),N_array,ratio_d_and_w);
    Yy(L)=(abs(Wopt'*Z2))^2;
end
subplot(2,1,1)
plot(x_axis,Yy);
xlabel('角度')
ylabel('Capon线性波束图')
subplot(2,1,2)
plot(x_axis,N_array*log10(Yy));
xlabel('角度')
ylabel('Capon对数波束图')
axis([-100 100 -300 0]);
grid on;
title('distance/wavelength=0.5时的方向图 SNR=5;INR=35')