%Demo done during Session 14 illustrating
%adaptive spatial filtering for digital communications.
%First adaptation is done based on angle-or-arrival only
%then we switch to decision-directed mode
clear all
clf
set(0,'defaultaxesfontsize',22);
M=12;
Nb=256;
P=3;
angles=[45 130 95]*(pi/180);
powers=[10 10 1];
%angles=[80 90 115]*(pi/180);
dlambda=0.5;
A=zeros(M,P);
X=zeros(M,Nb);
Rideal=zeros(M,M);
%generate random bits of information
%and create data matrix X
for k=1:P,
mu=pi*cos(angles(1,k));
a=exp(j*mu*(0:M-1)).';
Rideal=Rideal+a*a';
br=ones(1,Nb);
temp=rand(1,Nb);
br(find(temp<.5))=-1;
bi=ones(1,Nb);
temp=rand(1,Nb);
bi(find(temp<.5))=-1;
b=br+j*bi;
X=X+power(1,k)*a*b;
end
%add some noise
X=X+0.1*(randn(M,Nb)+j*randn(M,Nb));
%estimation of optimal filter coefficient vector
%based on block time-averaged estimates of Rxx and rdx
Rxx=zeros(M,M); 
w=1; alpha=1/w;
for n=1:Nb;
Rxx=Rxx+X(:,n)*X(:,n)';
end
%initial filter coefficient vectors for LMS and RLS
RI=inv(Rideal+.01*eye(M));
hest=RI*a;
%hest=a;
%Create Blocking Matrix:
a=a/M;
B=null(a');
U=[a B];
%
%Initialize LMS
wlms=zeros(M-1,1);
wlms(M/2,1)=1;
wrls=wlms;
Rxx=Rxx/Nb;
%compute step size for LMS
Ruu=B'*Rxx*B;
mu=1/(4*trace(Ruu));
Rinv=0.1*eye(M-1,M-1);
%
%run LMS and RLS for first Ns symbols using known
% direction only --no training sequence
%
Rxxs=zeros(M,M);
Ns=64;
for n=1:Ns;
%ideal case first
y(1,n)=hest'*X(:,n);
Rxxs=Rxxs+X(:,n)*X(:,n)';
%
s(n)=a'*X(:,n);
u=B'*X(:,n);
%LMS adaptation
ylms(1,n)=s(n)-wlms'*u;
wlmsold=wlms;
wlms=wlmsold+mu*conj(ylms(1,n))*u;
%RLS adaptation
wrlsold=wrls;
errrls(1,n)=s(n)-wrlsold'*u;
Rinvold=Rinv;
f=Rinvold*u;
murls=f'*u;
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f');
wrls=wrlsold+conj(errrls(1,n))*Kvec;
hrls=a+B*wrls;
yrls(1,n)=hrls'*X(:,n);
end
hlms=a+B*wlms;
%
%switch over to Decision-Directed Mode
%
c1=a'*Rxxs*a;  
bc=B'*Rxxs*a;  
c2=c1-bc'*Rinv*bc;
Rxxinv=[1 -bc'*Rinv ; -Rinv*bc Rinv*bc*bc'*Rinv+c2*Rinv]/c2;
%
%Rinv=inv(Rxxs);
Rinv=U*Rxxinv*U';
norm(Rinv-inv(Rxxs))
%
%Rinv=inv(Rxxs/Ns);
for n=Ns+1:Nb;
%ideal case first
y(1,n)=hest'*X(:,n);
%
%LMS adaptation
ytlms(1,n)=hlms'*X(:,n);
if real(ytlms(1,n)) > 0, btr=1; end
if real(ytlms(1,n)) < 0, btr=-1; end
if imag(ytlms(1,n)) > 0, bti=1; end
if imag(ytlms(1,n)) < 0, bti=-1; end
bt=btr+j*bti;
errlms=bt-ytlms(1,n);
hlmsold=hlms;
hlms=hlmsold+mu*conj(errlms)*X(:,n);
ylms(1,n)=hlms'*X(:,n);
%RLS adaptation
ytrls(1,n)=hrls'*X(:,n);
if real(ytrls(1,n)) > 0, btr=1; end
if real(ytrls(1,n)) < 0, btr=-1; end
if imag(ytrls(1,n)) > 0, bti=1; end
if imag(ytrls(1,n)) < 0, bti=-1; end
bt=btr+j*bti;
hrlsold=hrls;
errrls=bt-hrlsold'*X(:,n);
Rinvold=Rinv;
f=Rinvold*X(:,n);
murls=f'*X(:,n);
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f');
hrls=hrlsold+conj(errrls)*Kvec;
yrls(1,n)=hrls'*X(:,n);
end
%
ell=0;
for theta=0:.1:180,
ell=ell+1;
atheta=exp(j*pi*cos(theta*(pi/180))*(0:M-1));
Shest(ell)=20*log10(abs(hest'*atheta.'));
Shlms(ell)=20*log10(abs(hlms'*atheta.'));
Shrls(ell)=20*log10(abs(hrls'*atheta.'));
end
Shest=Shest-max(Shest);
Shlms=Shlms-max(Shlms);
Shrls=Shrls-max(Shrls);
theta=0:.1:180;
plot(theta,Shest,'k','Linewidth',4);
axis([0 180 -40 0]);
hold on
plot(theta,Shlms,'b','Linewidth',4);
plot(theta,Shrls,'r','Linewidth',4);
angdeg=angles*(180/pi);
for k=1:P,
plot([angdeg(1,k) angdeg(1,k)], [-40 0], 'm','Linewidth',4);
end
legend('MV with AOA','LMS Estimate', 'RLS Estimate');
title('Adapted Beam Patterns');
xlabel('Angle (Degrees)'); 
ylabel('Beam Magnitude (dB)');
hold off
pause
%plot data sequence obtained after beamforming
plot(real(b),imag(b),'mx','MarkerSize',20,'Linewidth',4);
axis([-2 2 -2 2]);
hold on
Ns=64;
plot(real(y),imag(y),'kx','MarkerSize',20,'Linewidth',4);
plot(real(ylms(1,Ns:Nb)),imag(ylms(1,Ns:Nb)),'b+','MarkerSize',20,'Linewidth',4);
plot(real(yrls(1,Ns:Nb)),imag(yrls(1,Ns:Nb)),'ro','MarkerSize',20,'Linewidth',4);
legend('True Bits','MV with AOA','LMS Estimate', 'RLS Estimate');
title('LMS and RLS Output Constellation');
hold off
pause
%spatial spectrum estimates
Rbb=Rxx(M:-1:1,:);  Rb=conj(Rbb(:,M:-1:1));
Rfb=0.5*(Rxx+Rb);
Rinv=inv(Rfb);
%Eigenanalysis
[E,D,V]=svd(Rfb);
Pn=E(:,P+1:M)*E(:,P+1:M)';
%compute estimate of spectrum, plot and compare
%with true spectrum
%
ell=0;
for theta=0:.1:180,
ell=ell+1;
atheta=exp(j*pi*cos(theta*(pi/180))*(0:M-1));
SMinVar(ell)=-10*log10(real(conj(atheta)*Rinv*atheta.'));
SMUSIC(ell)=-10*log10(real(conj(atheta)*Pn*atheta.'));
end
theta=0:.1:180;
plot(theta,(SMUSIC-max(SMUSIC)),'r','Linewidth',4);
axis([0 180 -60 0]);
hold on
plot(theta,(SMinVar-max(SMinVar)),'b','Linewidth',4);
for k=1:P,
plot([angdeg(1,k) angdeg(1,k)], [-60 0], 'm','Linewidth',4);
end
legend('MUSIC','Min Variance');
title('Spatial Spectrum Estimates');
xlabel('Arrival Angle (Degrees)'); 
ylabel('Spectral Magnitude (dB)'); 
hold off