%Demo done during Session 12 illustrating
%adaptive spatial filtering for digital communications.
clear all
clf
set(0,'defaultaxesfontsize',22);
M=12;
Nb=32;
P=3;
angles=[25 80 135]*(pi/180);
%angles=[80 90 115]*(pi/180);
dlambda=0.5;
%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+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;
%Initialize LMS
hlms=zeros(M,1);
hrls=hlms;
Rxx=Rxx/Nb;
mu=2/(2*trace(Rxx));
Rinv=0.1*eye(M,M);
%iterating over successive data blocks of length M
for n=1:Nb;
y(1,n)=hest'*X(:,n);
%LMS adaptation
ylms(1,n)=hlms'*X(:,n);
errlms=b(n)-ylms(1,n);
hlmsold=hlms;
hlms=hlmsold+mu*conj(errlms)*X(:,n);
%RLS adaptation
hrlsold=hrls;
errrls=b(n)-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
plot(real(y),imag(y),'kx','MarkerSize',20,'Linewidth',4);
plot(real(ylms),imag(ylms),'b+','MarkerSize',20,'Linewidth',4);
plot(real(yrls),imag(yrls),'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