%Maximum SIR beamforming
d=.5;
N= 5;
sig2=.001;        % noise variance
theta=-pi/2:.01:pi/2;
ang=theta*180/pi;
th0=30*pi/180;      % receive angle
th1=-20*pi/180;          % first interferer angle
th2=40*pi/180;           % second interferer angle

s=1;                %signal amplitude

a=[];
n1=[];
n2=[];
for n=1:N
   a=[a; exp(1j*pi*(n-(N+1)/2)*d*sin(th0))];       %desired signal array vector
    n1=[n1; exp(1j*pi*(n-(N+1)/2)*d*sin(th1))];    % interferer 1 array vector
  n2=[n2; exp(1j*pi*(n-(N+1)/2)*d*sin(th2))];      % interferer 2 array vector
end    
Rss=a*a.';              % signal correlation matrix  
R11=n1*n1';             % interferer 1 correlation matrix  
R22=n2*n2';             % interferer 2 correlation matrix  
A=[n1 n2];
Rnn=sig2*eye(N);        % noise correlation matrix       
Rii=A*A';               % total interferer correlation matrix
Ruu=Rii+Rnn;            % total undesired signal correlation matrix
Rxx=Ruu+Rss;            % total array correlation matrix

w=inv(Rxx)*a;           % calculated weight vector

for j=1:length(theta)
    th=theta(j);
aa=[];
 for n=1:N
     aa=[aa; exp(1j*pi*(n-(N+1)/2)*d*sin(th))];  
 end
 %aa=[exp(1j*2*pi*d*sin(th)) exp(1j*pi*d*sin(th)) 1 exp(-1j*pi*d*sin(th)) exp(-1j*2*pi*d*sin(th))].';
    y(j)=w'*aa;
end
figure;
plot(ang,abs(y)/max(abs(y)),'k')
axis([-90 90 0 1])
xlabel('\theta')
ylabel('|AF(\theta)|')
num2str(Rxx,2)