space_tempral.m

时空级联算法用于海面低空目标相干信号源的DOA估计

clear
clc
%%%%%%%%%%%%input%%%%%%%%%%%%%
N=10;%number of elements 
M=128;%snapshots
hr=20;%radar height
ht=400;%target height
d=0.5;%space in lambda unit
snr=30;%信噪比
p=0.9;%反射系数
f=500;%多普勒频率
deltat=10^-3;%采样时间
Rd=10^4;
%%%%%%%%%%%input%%%%%%%%%%%%%%%

fai=0.2*pi;
sintheta1=(ht-hr)/Rd;
theta1=asin(sintheta1);
theta2=-asin(2*hr/Rd+sintheta1);
arg=[theta1',theta2'];
sita=sin(arg); 
windows=hanning(N);%加窗，N点海明窗
w=ones(1,N);%权向量
wq=windows;

%方向图（加窗，瞄准方向0度）
for nn=0:round((80+80)/0.05)
    n=1:N;
    A0=exp(-i*(2*pi*d*(n-1)*sin((-80+nn*0.05)*pi/180)));
    F(nn+1)=0;
    F(nn+1)=abs(w*A0.');
end
F0=F/max(F);
figure(1)
theta0=-80:0.05:80;
plot(theta0,F0,'k')
xlabel('deg')
ylabel('G/dB')
grid on
print -dtiff DBF
%波束宽度
thetaB=50.8/(N*d)

for m=1:M
    for n=1:N
        A(n,1)=exp(-i*(2*pi*d*(n-1)*sita(1,1))); %方向矩阵
        A(n,2)=exp(-i*(2*pi*d*(n-1)*sita(1,2)));
    end
    b=exp(i*2*pi*(m-1)*deltat*f);
    s0=10.^(snr/20)*b;
    S=s0*[1;p*exp(i*fai)];
    x0=A*S+randn(N,1)+i*randn(N,1);
    y=wq'*x0;%输出数据向量
    y1(m)=y;
end 
y2=abs(fft(y1));%对波束形成数据作FFT变换
[xx,ss]=max(y2);
a0=(1/M)*(1/deltat);
aa=1:M;
ss0=aa*a0;
figure(2)
plot(ss0,y2,'k') 
xlabel('frequency/Hz')
print -dtiff Doppler
fd=ss*a0;

X=zeros(N,1);
for m=1:M
    fei(m)=exp(-i*2*pi*deltat*fd*(m-1));
    x=x0*fei(m);   
    X=X+x;
end

N1=-4;N2=4;tt=0.001;
for jj=0:round((N2-N1)/tt)
    theta01=N1+jj*tt;
    theta02=-asin(sin(theta01*pi/180)+2*hr/Rd);
    aa1=2*pi*d*sin(theta01*pi/180);
    aa2=2*pi*d*sin(theta02);
    n=1:N;
    e1=exp(-i*aa1*(n-1));
    e2=exp(-i*aa2*(n-1));
    An=[e1.',e2.'];
    PA=An*inv(An'*An)*An';
    TT=X'*PA*X;
    T(jj+1)=0;
    T(jj+1)=10*log10(TT);
end  

figure(3)    
pp=N1:tt:N2;
plot(pp,T,'k');
xlabel('\theta/deg')
ylabel('log-likelihood')
grid on
print -dtiff Elevation