%---------------------------------------------------------
%--
%---------------------------------------------------------
clc
clear

Imag=sqrt(-1);
array_num=10;                                                               %sensor number
source_num=3;                                                               %source number
lmda=0.03;                                                                  %unit is:meter
d=lmda/2;                                                                   %sensor space
cx=[0:array_num-1]'*d;                                                      %sensor Xcord
cy=[0:array_num-1]'*0;                                                      %sensor Ycord
position=[cx cy];

doa=[0 -20 40]*pi/180;                                 %treat the last two source as a jammer                           
A=exp(Imag*(2*pi/lmda)*position*[sin(doa);cos(doa)]);   %array steer matrix

snap_num=256;                                           %snapshot number
ChirpRate=[50 60 150]';                                 %chirp rate for source signal  
Fs=1.0e+3;                                              %smaple frequency
SampleTime=[0:1/Fs:1];
SampleNum=length(SampleTime);
s=exp(Imag*pi*ChirpRate*(SampleTime.^2));               %the baseband source signal

PhiAmpErrorConsider=0;                                   
if PhiAmpErrorConsider
    DeltaAm=0.1;                                                           
    DeltaPh=5;                                          %unit is deg                                                        
    Ame=[1 1+(rand(1,array_num-1)-0.5)*DeltaAm*sqrt(12)];                  
    Phe=[0 (rand(1,array_num-1)-0.5)*DeltaPh*sqrt(12)];                     
else
    Ame=ones(1,array_num);
    Phe=zeros(1,array_num);
end
G=diag(Ame.*exp(Imag*Phe*pi/180));                       

NoiseConsider=1;                                                        
if NoiseConsider
        Noise=sqrt(0.5)*(randn(array_num,snap_num)+Imag*randn(array_num,snap_num));
else
        Noise=zeros(array_num,snap_num);
end

EnvelopeConsider=1;                                                       
if EnvelopeConsider                                                         
     TimeWidth=[32 32 32];                                                %Time Width.unit is sample point
    for i=1:source_num
       Const=ceil(SampleNum/TimeWidth(i));
       envelope(i,:)=reshape(ones(TimeWidth(i),1)*(2*randint(1,Const)-1),1,TimeWidth(i)*Const);
    end
else 
    envelope=ones(source_num,SampleNum);
end
%-------------------------------------
%---get the receive data------
%-------------------------------------
SNR=[20 20 40];                                         %Signal Noise Rate.dB
Am_n=mean2(abs(Noise));
Am=Am_n*10.^(SNR/20);
S=envelope(:,1:SampleNum).*(diag(Am)*s);  
Data=G*A*S(:,1:snap_num)+Noise;                         %sensor array received data
Data_sig=G*A(:,1)*S(1,1:snap_num);                      
Data_jamn=G*A(:,source_num-1:source_num)*S(source_num-1:source_num,1:snap_num)+Noise;  %jammer and noise data
Rxx=Data*Data'/snap_num;  
Rx_sig=Data_sig*Data_sig'/snap_num;
Rx_jamn=Data_jamn*Data_jamn'/snap_num;
Rxs_vec=Data*S(1,1:snap_num)'/snap_num;
InvRxx=inv(Rxx);

Wopt_s=G*A(:,1);             %static weight. Wopt_s=ste_vec(doa(sig))

[V,D]=eig(Rx_sig,Rx_jamn); [D,Indice]=sort(abs(diag(D))); V=V(:,Indice);
Wopt_SNR=V(:,length(D));     %the max SNR weight
Wopt_MSE=InvRxx*Rxs_vec;
% Wopt_LCMV=InvRxx*Wopt_s./(Wopt_s'*InvRxx*Wopt_s);
Wopt_LCMV=inv(Rx_jamn)*Wopt_s./(Wopt_s'*InvRxx*Wopt_s);

Wopt_mat=[Wopt_s,Wopt_SNR,Wopt_MSE,Wopt_LCMV];
[nrow,ncol]=size(Wopt_mat);
%-------------------------------------
%----get the beamformer response----
%-------------------------------------
space_seek=0.05;
doa_seek=[-90:space_seek:90]*pi/180;
seek_num=length(doa_seek);
PdB=zeros(ncol,seek_num);
for k=1:ncol
    P_vec=[];
    for i=1:seek_num
       st_vec=exp(Imag*(2*pi/lmda)*position*[sin(doa_seek(i));cos(doa_seek(i))]);  %search steer vector  
       P_vec=[P_vec,st_vec'*Wopt_mat(:,k)*Wopt_mat(:,k)'*st_vec];                  %get the array response
    end
    P_vec=abs(P_vec);
    PdB(k,:)=10*log10(P_vec./max(P_vec)+1.0e-6);
end                                                             

figure;
hold on;
for i=1:source_num
    plot(doa(i)*180/pi,[-80:0.1:0],'-r');
end 
CurveNo=[1:ncol];
fid(1)=plot(doa_seek*180/pi,PdB(CurveNo(1),:),'-g');
fid(2)=plot(doa_seek*180/pi,PdB(CurveNo(2),:),'-.r');  
fid(3)=plot(doa_seek*180/pi,PdB(CurveNo(3),:),'-k'); 
fid(4)=plot(doa_seek*180/pi,PdB(CurveNo(4),:),':b'); 
axis([-90 90 -80 1]);
grid on
legend(fid,'Wopt','WoptSNR','WoptMSE','WoptLCMV');