%%%%%%%%%%%%%%%%%%%
%% Conjugate Gradient %%
%%%%%%%%%%%%%%%%%%%

%----- Givens -----%
K=20;    % total number of data samples
sig2=.001;
d = .5;  % element spacing in terms of wavelength d = lambda/2
N=8;
%randn('seed')
n=sqrt(sig2)*randn(N,K);
%thetaS = input('   What is the desired users AOA (in degrees)?   ');
thetaS=45;
%thetaI = input('   What is the interferers AOA(in degrees)?   ');
theta1 = -30;
%----- Desired Signal & Interferer -----%
TT=1E-3;
tt=[1:K]*TT/K;
S=cos(pi*tt/TT);
thetaS = thetaS*pi/180;                  % desired user AOA
I1=randn(1,K);
theta1 = theta1*pi/180;  % interferer AOA
I2=randn(1,K);
theta2=0;
theta2=theta2*pi/180;
vS = []; v1 = []; v2=[];
ij=1:N;
vS=exp(1j*(ij-1)*2*pi*d*sin(thetaS));
v1=exp(1j*(ij-1)*2*pi*d*sin(theta1));
v2=exp(1j*(ij-1)*2*pi*d*sin(theta2));
x=(vS.'*S+v1.'*I1+v2.'*I2)+n;	%array input with noise
w=.1*ones(N,1)*(1+1j);			% initialize weights
A=x.';					% define matrix of array values for K time samples
r=S.'-A*w;				% calculate residual vector
rold=r;
D=A'*r;
rit=zeros(1,N);
for jj=1:K
   nu=A'*r;
   den=A*D;
   mu=real(nu'*nu)/real(den'*den);
   w=w-mu*D;
   r=r+mu*A*D;
   den=A'*rold;
   nu=A'*r;
   al=real(nu'*nu)/real(den'*den)
rold=r;
D=A'*r-al*D;
rit(jj)=r'*r;
rit(jj)=norm(r);
dis(jj)=norm(D);
end
figure;plot(1:K,rit,'k')
xlabel('Iteration no.')
Ylabel('||r||')

w=w/abs(w(1));

 %----- Plot Results -----%
theta = -pi/2:.01:pi/2;
AF = zeros(1,length(theta));

% Determine the array factor for linear array with new weights
for ii = 1:N
   ww=w(ii);
    AF = AF + ww*exp(1j*(ii-1)*2*pi*d*sin(theta));
end

figure;
plot(theta*180/pi,abs(AF)/max(abs(AF)),'k')
xlabel('AOA (deg)')
ylabel('|AF|')
axis([-90 90 0 1.1])
set(gca,'xtick',[-90 -60 -30 0 30 60 90])
grid on

% solve using the MATLAB conjugate gradient canned program LSQR(A,d)

www=lsqr(A,S');

ii=1:N;
a=exp(1j*(ii-1)'*2*pi*d*sin(theta));
AF=www.'*a;
figure;
plot(theta*180/pi,abs(AF)/max(abs(AF)),'k')
xlabel('AOA (deg)')
ylabel('|AF|')
axis([-90 90 0 1.1])
set(gca,'xtick',[-90 -60 -30 0 30 60 90])
grid on