% Multiple beamforming. 
clear;clc;close all;
starttime=now;
%%%%%%%%%%%%%%%%%%The initial condition%%%%%%%%%%%%%%%%%%%%
N=8;                                            % the array elments no.
d=0.5;                                          % in wavelength unit.
Var=1;                                          % white noise variance.
p=1;                                            % no. of desired signal
SNR=[10 40]';                                   % source SNR in dB.
Angs1=[0.2 -0.4]';                              % the incident angle vector in radians.
L=length(Angs1);                                % sources no.
snapshot=200;                                   % Sampling data number.
sir_data=200;                                   % Iteration numbers of LMS algorithm
experimentno=20;                                %  Averaging experimentno.
g=ones(p,1);
ave_sir=zeros(1,sir_data);
sir_end=0;
step=10^(-7);                                   % Stepsize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j=sqrt(-1);
n=[0:1:N-1]';
pj=zeros(L,1);
for i=1:L
    pj(i)=Var*10^(SNR(i)/10);
end
amp=sqrt(pj);
%%%%%%%% Far-Field %%%%%%%%
angf=2*pi*d*Angs1;
aaf=n*angf';
D=exp(j*aaf);                                   % the Direction Matrix.                 
%%%%%%%%%%%%%%%% generate the autocorrelation matrix %%%%%
ny=zeros(N,1);
for ite=1:experimentno
    ite
    x=zeros(N,snapshot);
    %%%%% Generate the BPSK source %%%%%
    %     sam_t=ones(L,1);                                % the normalized sampling interval.
    %     pha=zeros(L,1);                                 % sources initial phase.
    %     chg_pha=ones(L,1);                              % Data length to exchage BPSK.
    %     factor=ones(L,1);
    %     for k=1:snapshot
    %         for kk=1:L
    %             if rem(k,chg_pha(kk)) == 0
    %                 b_level=round(rand);
    %                 factor(kk)=round(b_level-0.5)*factor(kk);
    %             end
    %         end  
    %         bpha=pi/2*factor;
    %         source(:,k)=amp.*exp(j*(2*pi*(k-1)*sam_t+bpha));
    %     end
    %%%%%%% Genrate project's source %%%%%%
    source=kron(amp,ones(1,snapshot));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x=D*source;
    samn=randn(N,snapshot)+j*randn(N,snapshot);
    samn=samn/sqrt(2)*Var;
    y=x+samn;
    Rxx=y*y'/snapshot;
    %%%%%%%%% MVDR %%%%%%%%%
    w=inv(Rxx)*D(:,1:p)*inv(D(:,1:p)'*inv(Rxx)*D(:,1:p))*g;         % MVDR weighting vector
    mvdr(ite,:)=pattern(n*d,w,ny);
    sir_end=sir_end+Sir(D,pj,w,L,p);
    %%%%%%% MVDR-LMS %%%%%%%
    w_int=D(:,1:p)*inv(D(:,1:p)'*D(:,1:p))*g;
    sir_int(1)=Sir(D,pj,w_int,L,p);
%     wq=w_int;
%     P=eye(N)-D(:,1:p)*inv(D(:,1:p)'*D(:,1:p))*D(:,1:p)';
    for data=2:sir_data
        error=source(1,data)-w_int'*y(:,data);
        w_int=w_int+step*conj(error)*y(data);                               % Haykin
%         w_int=P*(w_int-satep*y(:,data)*conj(w_int'*y(:,data)))+wq;           % Frost
        sir_int(data)=Sir(D,pj,w_int,L,p);
    end
    mvdr_lms(ite,:)=pattern(n*d,w_int,ny);
    ave_sir=ave_sir+sir_int;
end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mvdr=abs(mvdr).^2;
pat=Normalization(mvdr,experimentno);

mvdr_lms=abs(mvdr_lms).^2;
pat2=Normalization(mvdr_lms,experimentno);

ave_sir=10*log10(ave_sir/experimentno);
sir_end=10*log10(sir_end/experimentno);

s1=mat2str(SNR');
s2=mat2str(Angs1');

figure
k=sin([-90:0.5:90]'*pi/180);
k1=[1:1:sir_data];
plot(k,pat,'.-g',k,pat2,'.-b');
legend('MVDR','MVDR-LMS');
xlabel('Angles');
ylabel('Power Gain in dB')
title(sprintf('Array elements# =%g, SNR=%sdB from degree %s, experiment# =%g, snapshot=%g\nOTIR of MVDR=%gdB, OTIR of MVDR-LMS after %g iteration=%gdB',N,s1,s2,experimentno,snapshot,sir_end,sir_data,sir_int(end)))
for u=1:L
    line(Angs1(u),[-70:0.1:0])
end
line([-1:0.002:1],0);
axis([-1 1 -70 10]);

totaltime=now-starttime;
tt=datestr(totaltime,13)