test_ls_xianzhen0.m

分别仿真16阵元的十字阵列、线阵、面阵和圆阵的LS算法的方向图

% --------- code rate=6.4e6 (20080707)-wanglei---------
% LS算法可以采用2～3倍的阵元个数的快拍数参与运算
% 收敛结果也比较好。降低运算量。
clc
clear all
close all
fs=20e6;  % sampling freq
fc=0;         % 
N=1200;
NN=700;
load s1.mat   % 训练序列
load s2.mat   % interfere signal,AM
load s4.mat
% load s6.mat
sss=[s4(1:NN) s1(1:N) s4(1:1e4)];; % 发射数据  s1训练序列，s4有用信息
s3=hilbert(s2);  % change to complex
%s3、s6干扰信号
sir=-10;
pii=10^(-sir/20);
s=[sss;pii*s3(1:N+NN+1e4)]; ;  %s3干扰信号

a1=90; %有用信号方向
a2=50;  % 干扰信号方向
fai1=30;%有用信号俯仰角
fai2=50;% 干扰信号俯仰角
a3=-20;

M=16;
% a11=exp(-j*pi*[-M:M]'*(cos(a1*pi/180)*cos(fai1*pi/180)));
% a22=exp(-j*pi*[-M:M]'*(sin(a1*pi/180)*cos(fai1*pi/180)));
% a=kron(a11,a22);
a=[exp(-j*pi*[0:M-1]'*(cos(a1*pi/180)*cos(fai1*pi/180)))] ;
b=[exp(-j*pi*[0:M-1]'*(cos(a2*pi/180)*cos(fai2*pi/180)))] ;
% c=exp(-j*pi*[0:M-1]'*sin((a1+5)*pi/180));
A=[a,b];

snr=10;      % snr: in dB
randn('state',000)
x=A*s+10^(-snr/20)*randn(16,N+NN+1e4);
  tic
K=400;  %一个周期的采样点
K1=128;
K2=80;  % 参与LS运算的快拍数

gold_rate=6.4*1e6;
w0=eye(16,1);  
t1=0; % 采样开始时刻, 20 samples shift
T1=5e-3; % 采样持续时间
m=0;
%--------xiangguan LS ------------
for n=1:6
    y=w0'*x(:,1+(n-1)*K:K+(n-1)*K);
   [yy,ptinit1]=xiangguan(y,fs,gold_rate,1,K,8,K1);
%    ptinit1=ptinit1+3;
   aa(n)=max(abs(yy));
   if aa(n)>4
       delay1(n)=ptinit1/fs; 
       sr=signalgen(1,t1,fs,T1,delay1(n),fc,gold_rate,K1);
       break
   end
end    % end of xiangguan
r=sr(1:K2);
tic
w=inv((x(:,1+(n-1)*K:K2+(n-1)*K))*(x(:,1+(n-1)*K:K2+(n-1)*K))')*x(:,1+(n-1)*K:K2+(n-1)*K)*r';
%  w=inv((x(:,1+(n-1)*K:K+(n-1)*K))*(x(:,1+(n-1)*K:K+(n-1)*K))')*x(:,1+(n-1)*K:K+(n-1)*K)*r';
toc
% w=inv((x(:,1+(n-1)*K:K2+(n-1)*K))*(x(:,1+(n-1)*K:K2+(n-1)*K))')*a;
% --------------- 方向图output ------------------
% for i=1:4
%     Grow(i,:)=conj(w44(i,:))*A;
%     Gcol(i,:)=w44(:,i)'*B;
% end
%     GROW=sum(Grow);
%     GCOL=sum(Gcol);
%     G=GROW.*GCOL;
%     G=reshape(G,100,100);

Q=100;
theta=linspace(0,pi,Q);
fai=linspace(0,pi/2,Q);
% for m=1:100
%     for k=1:100
%         G1=[exp(-j*pi*[0:M-1]'*cos(theta(m))*cos(fai(k)))] ;
%         G(m,k)=w'*G1;
%     end
% end

for m=1:100
    for k=1:100
        G1=[exp(-j*pi*[0:M-1]'*cos(theta(m))*cos(fai(k)))] ;
        G(m,k)=w'*G1;
    end
end


    mesh(theta*180/pi,fai*180/pi,20*log10(abs(G)))
    hold on
   plot3([a1,a1],[fai1,fai1],[-50,10],'k','linewidth',1.5)
plot3([a2,a2],[fai2,fai2],[-50,10],'linewidth',1.5)
% mesh(abs(fft2(w,100,100)))
%     a1=-10; %有用信号方向
% a2=30;  % 干扰信号方向
% fai1=20;%有用信号俯仰角
% fai2=40;% 干扰信号俯仰角
% mesh(abs(fft2(w,100,100)))
% hold on
% plot3([a1,a1],[fai1,fai1],[-50,0],'linewidth',1.5)
% plot3([a2,a2],[fai2,fai2],[-50,0],'linewidth',1.5)
% P1=20*log10(abs(G(:,20))./max(abs(G(:,20))));
% plot(P1)
% figure
% g1=plot(phi*180/pi,P1,'b-','linewidth',1.5);
% hold on
% plot([a1,a1],[-50,0],'b--','linewidth',1.5);
% plot([a2,a2],[-50,0],'b--','linewidth',1.5);
% % axis([-90,90,-40,0]);
% xlabel('\fontsize{12}DOA');
% ylabel('\fontsize{12}Power');
% title('LS(带训练序列)');
% 
figure(2)
plot(fai*180/pi,20*log10(abs(G)))
figure(3)
plot(theta*180/pi,20*log10(abs(G.')))