test_ls_yuanzhen.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=50;%有用信号俯仰角
fai2=20;% 干扰信号俯仰角
a3=-40;

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);

  D=[0:1:M-1]';          %------- 阵元半径为整波长，均匀圆阵

for i=0:M-1
phi(:,i+1)=cos(a1*-2*pi*i/M);
end

for i=0:M-1
a(:,i+1) = exp(j*2*pi*sin(fai1)*phi(:,i+1));
end

 for i=0:M-1
phi(:,i+1)=cos(a2*-2*pi*i/M);
end

for i=0:M-1
b(:,i+1) = exp(j*2*pi*sin(fai2)*phi(:,i+1));
 end
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
% --------------- 方向图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);

sita=(-pi/2:pi/180:pi/2);%---------仰角
Asita=sita*180/pi;

psai=(-pi/2:pi/180:pi/2);%----------方位角
Bpsai=psai*180/pi;
for i=0:M-1
phi1(:,i+1)=cos(theta.'-2*pi*i/M);
end

for i=0:M-1
    for h=1:100
vv(:,h,i+1) = exp(j*2*pi*sin(fai(h))*phi1(:,i+1));
    end
end



for i=1:100
    for h=1:100
      aa1=reshape(vv(i,h,:),M,1);
      Bself(i,h)=10*log10(abs(w'*aa1).^2);
end
end
% 
% %---------------------------------------------------------------
% figure;
% set(gcf,'color','white');
% mesh(Bpsai,Asita,Bself);  %----------横坐标为方位角，纵坐标为仰角
% xlabel('X轴');
% Ylabel('Y轴');
% zlabel('dB');
figure;
set(gcf,'color','white');
mesh(theta*180/pi,fai*180/pi,Bself); 
hold on;%----------横坐标为方位角，纵坐标为仰角
xlabel('仰角 (度)');
Ylabel('方位角 (度)');
zlabel('增益 (dB)');

   plot3([a1,a1],[fai1,fai1],[-50,10],'k','linewidth',1.5)
plot3([a2,a2],[fai2,fai2],[-50,10],'linewidth',1.5)




% 
figure(2)
plot(fai*180/pi,20*log10(abs(G)))
figure(3)
plot(theta*180/pi,20*log10(abs(G.')))