%Demo done during Session 12 illustrating
%adaptive spatial filtering for digital communications.
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);
%generate random bits of information
X=zeros(M,Nb);
Rideal=zeros(M,M);
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+powers(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);
rdx=zeros(M,1);
w=1; alpha=1/w;
for n=1:Nb;
Rxx=Rxx+X(:,n)*X(:,n)';
rdx=rdx+conj(b(n))*X(:,n);
end
%initial filter coefficient vectors for LMS and RLS
RI=inv(Rideal+.001*eye(M));
%RI=inv(Rxx);
%hest=RI*rdx;
hest=RI*a;
%hest=a;
%Create Blocking Matrix:
a=a/M;
B=null(a');
%
%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);
%iterating over successive data blocks of length M
for n=1:Nb;
%ideal case first
y(1,n)=hest'*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
%create h from w for LMS
hlms=a-B*wlms;
%
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