test_ls_fangzhen0.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=-40;

M=4;
a11=exp(-j*pi*[0:M-1]'*(cos(a1*pi/180)*cos(fai1*pi/180)));
a22=exp(-j*pi*[0:M-1]'*(sin(a1*pi/180)*cos(fai1*pi/180)));
a=kron(a11,a22);

b11=exp(-j*pi*[0:M-1]'*(cos(a2*pi/180)*cos(fai2*pi/180)));
b22=exp(-j*pi*[0:M-1]'*(sin(a2*pi/180)*cos(fai2*pi/180)));
b=kron(b11,b22);
% 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(4*M,N+NN+1e4);
  tic
K=400;  %一个周期的采样点
K1=128;
K2=80;  % 参与LS运算的快拍数

gold_rate=6.4*1e6;
w0=eye(4*M,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:K+(n-1)*K))*(x(:,1+(n-1)*K:K+(n-1)*K))')*a;
% --------------- 方向图output ------------------
w44=reshape(w,4,4);              
Q=100;
theta=linspace(0,pi,Q);
fai=linspace(0,pi/2,Q);
A=exp(-j*pi.*[0:M-1]'*(kron(cos(theta),cos(fai))));
B=exp(-j*pi.*[0:M-1]'*(kron(sin(theta),cos(fai))));
% 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);


for m=1:100
    for k=1:100
        Grow=exp(-j*pi*[0:M-1]'*(cos(theta(m))*cos(fai(k))));
        Gcol=exp(-j*pi*[0:M-1]'*(sin(theta(m))*cos(fai(k))));
        G1=kron(Grow,Gcol);
        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],'linewidth',1.5)
plot3([a2,a2],[fai2,fai2],[-50,10],'linewidth',1.5)
%     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.')))
y=w'*x;
figure,plot(abs(fftshift((fft(y)))))