% Robust Beamforming Using Semidefinite Programming
% copyright Sergiy A. Vorobyov, Alex B. Gershman; May 22, 2001
% email: svor@mail.ece.mcmaster.ca

clear;
clc;
n_sampl=30;
n=10;  % sensors
dl=0.5; % d/lambda
sqrt05 = sqrt(0.5);
tet= [30, 50];
tet_s=3;
tet_s_r=3;
m  = length(tet);
p  = 1000;
SNR_dB = [-20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 ...
        -6 -5 -4 -3 -2 -1 0 1 2 3 4 5];
%SNR_dB = [-20 -19.5 -19 -18.5 -18 -17.5 -17 -16.5 -16 -15.5 -15 -14.5 -14 -13.5 -13 -12.5...
%        -12 -11.5 -11 -10.5 -10 -9.5 -9 -8.5 -8 -7.5 -7 -6.5 -6 -5.5 -5];
for i=1:length(SNR_dB), SNR(i) = 10^(SNR_dB(i)/10); end;
Am = exp(1i*2.0*[0:n-1]'*pi*dl*sin(tet*pi/180.0)); % interf. steering matr.

SINRs=zeros(1,length(SNR)); SINRr=zeros(1,length(SNR)); 
SINRrob=zeros(1,length(SNR)); SINRdl=zeros(1,length(SNR));

for iii=1:100;
    
    as = exp(1i*2.0*[0:n-1]'*pi*dl*sin(tet_s*pi/180.0));
    var=4;
    b1 = exp(1i*2.0*[0:n-1]'*pi*dl*sin((tet_s+rand*sqrt(3*var))*pi/180.0));
    b2 = exp(1i*2.0*[0:n-1]'*pi*dl*sin((tet_s+rand*sqrt(3*var))*pi/180.0));
    b3 = exp(1i*2.0*[0:n-1]'*pi*dl*sin((tet_s+rand*sqrt(3*var))*pi/180.0));
    b4 = exp(1i*2.0*[0:n-1]'*pi*dl*sin((tet_s+rand*sqrt(3*var))*pi/180.0));
    psi1 = 2*rand*pi; psi2 = 2*rand*pi; psi3 = 2*rand*pi; psi4 = 2*rand*pi;
    as_r = as + exp(1i*psi1)*b1 + exp(1i*psi2)*b2 + exp(1i*psi3)*b3 + exp(1i*psi4)*b4;
    as_r = sqrt(10)*as_r/norm(as_r)';
    
    Ncount=0;
    for ps=SNR;
        Ncount=Ncount+1;
nk=n_sampl;

Rs = zeros(n);
Rss = zeros(n);

for l=1:nk
     s = (randn(m,1) + 1i*randn(m,1)).*sqrt(p)'*sqrt05; % waveform
     ssig=(randn(1,1) + 1i*randn(1,1))*sqrt(ps)*sqrt05;  % signal sig^2=1
     z = (randn(n,1) + 1i*randn(n,1))*sqrt05;           % noise sig^2=1
     x = Am*s + z;     % snapshot without signal (radar case)
     xs = Am*s + ssig*as_r + z; % snapshot with signal (passive location case)
     Rs = Rs + x*x';       
     Rss = Rss + xs*xs';
end

     Rs=Rs/nk;  % sample cov. matrix
     Rss=Rss/nk;  % sample cov. matrix


w=inv(Rs)*as;            % weight vector (SMI)
ws=inv(Rss)*as;

% diagonal loading
Rssdl = Rss + eye(n,n)*0.5;
wsdl=inv(Rssdl)*as;

% robust beamformer (projection to signal and interference subspace (S&IS))
[RR,D]=eig(Rss);   % eigenvalue decomposition of covariance matrix
RR_subs=RR(:,n-2:n); 
as_new=[RR_subs,zeros(n,n-3)]*as;  % projection of waveform to S&IS
ws_r=inv(Rss)*as_new;   % beamformer

R=p*Am*Am'+eye(n,n);        % exact cov. matrix

% interior-point optimization based robust beamformer
wrob = robustbeam(Rss,as,3,0);

SINRs1(Ncount)=ps*(abs(ws'*as_r)*abs(ws'*as_r))/abs(ws'*R*ws);  % SINR
SINRdl1(Ncount)=ps*(abs(wsdl'*as_r)*abs(wsdl'*as_r))/abs(wsdl'*R*wsdl);
SINRr1(Ncount)=ps*(abs(ws_r'*as_r)*abs(ws_r'*as_r))/abs(ws_r'*R*ws_r);
SINRrob1(Ncount)=ps*(abs(wrob'*as_r)*abs(wrob'*as_r))/abs(wrob'*R*wrob); % SINR for robust beamformer
SINRopt(Ncount)=ps*as_r'*inv(R)*as_r;

 end
SINRs     = ((iii-1)/iii)*SINRs+(1/iii)*SINRs1;
SINRdl    = ((iii-1)/iii)*SINRdl+(1/iii)*SINRdl1;
SINRr     = ((iii-1)/iii)*SINRr+(1/iii)*SINRr1;
SINRrob   = ((iii-1)/iii)*SINRrob+(1/iii)*SINRrob1;
  
end;     
figure(6)
plot(SNR_dB,10*log10(SINRopt),'m',SNR_dB,10*log10(SINRrob),'-*',...
    SNR_dB,10*log10(SINRs),'-d',SNR_dB,10*log10(SINRdl),'-^',...
    SNR_dB,10*log10(SINRr),'-o'), grid
legend('OPTIMAL SINR','PROPOSED ROBUST BEAMFORMER','SMI BEAMFORMER',...
        'LSMI BEAMFORMER','EIGENSPACE-BASED BEAMFORMER');
Xlabel('SNR (DB)')
Ylabel('OUTPUT SINR (DB)')