%Demo done during Session 12 illustrating%adaptive spatial filtering for digital communications.clear allclfset(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 informationX=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 noiseX=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 rdxRxx=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 RLSRI=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 LMSwlms=zeros(M-1,1);wlms(M/2,1)=1;wrls=wlms;Rxx=Rxx/Nb;%compute step size for LMSRuu=B'*Rxx*B;mu=1/(4*trace(Ruu));Rinv=0.1*eye(M-1,M-1);%iterating over successive data blocks of length Mfor n=1:Nb;%ideal case firsty(1,n)=hest'*X(:,n);%s(n)=a'*X(:,n);u=B'*X(:,n);%LMS adaptationylms(1,n)=s(n)-wlms'*u;wlmsold=wlms;wlms=wlmsold+mu*conj(ylms(1,n))*u;%RLS adaptationwrlsold=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 LMShlms=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.'));endShest=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 onplot(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);endlegend('MV with AOA','LMS Estimate', 'RLS Estimate');title('Adapted Beam Patterns');xlabel('Angle (Degrees)'); ylabel('Beam Magnitude (dB)');hold offpause%plot data sequence obtained after beamformingplot(real(b),imag(b),'mx','MarkerSize',20,'Linewidth',4);axis([-2 2 -2 2]);hold onNs=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 offpause%spatial spectrum estimatesRbb=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.'));endtheta=0:.1:180;plot(theta,(SMUSIC-max(SMUSIC)),'r','Linewidth',4);axis([0 180 -60 0]);hold onplot(theta,(SMinVar-max(SMinVar)),'b','Linewidth',4);for k=1:P,plot([angdeg(1,k) angdeg(1,k)], [-60 0], 'm','Linewidth',4);endlegend('MUSIC','Min Variance');title('Spatial Spectrum Estimates');xlabel('Arrival Angle (Degrees)'); ylabel('Spectral Magnitude (dB)'); hold off