%Demo done during Session 17 illustrating%use of ESPRIT to estimate angle of arrival%of each source incident on uniform linear%array of antennasclear allclfset(0,'defaultaxesfontsize',22);M=12;Nb=32;P=3;%angles=[25 80 135]*(pi/180);angdeg=[80 90 115];angles=angdeg*(pi/180);dlambda=0.5;%generate random bits of informationpoles=zeros(1,P);X=zeros(M,Nb);Rideal=zeros(M,M);for k=1:P,mu=pi*cos(angles(1,k));poles(1,k)=exp(j*mu);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 noiseX=X+0.1*(randn(M,Nb)+j*randn(M,Nb));Rxx=X*X'/Nb;%for n=1:Nb;%Rxx=Rxx+X(:,n)*X(:,n)';%nd%Rbb=Rxx(M:-1:1,:);  Rb=conj(Rbb(:,M:-1:1));%Rxxfb=0.5*(Rxx+Rb);[E,D,V]=svd(Rxx);diag(D)%ESPRIT algorithm:Es=E(:,1:P);Es1=Es(1:M-1,:); Es2=Es(2:M,:);Psi=Es1\Es2;[T,Phi]=eig(Psi);Phivec=diag(Phi);%plot eigenvalues from ESPRIT and compare with true frequenciespolar(0,1,'.')hold onplot(real(poles),imag(poles),'kx','MarkerSize',18,'Linewidth',4);plot(real(Phivec),imag(Phivec),'ro','MarkerSize',18,'Linewidth',4);hold offlegend('True "poles"','ESPRIT eigenvalues')pause% %compute roots of each noise eigenvector and%to show that each has a root at the angle of%each incident source.%nc=0;s='kbgrmcykb';freqomega=linspace(-pi,pi,1024);theta=acos(freqomega/pi)*(180/pi);for nk=P+1:M;enroots=roots(E(:,nk));nc=nc+1;%plot eigenvalues from ESPRIT and compare with true frequenciespolar(0,1,'.')hold onplot(real(poles),imag(poles),'kx','MarkerSize',18,'Linewidth',4);plot(real(enroots),imag(enroots),'ro','MarkerSize',18,'Linewidth',4);title('Roots of Noise Eigenvector'),...hold offpauseendPn=E(:,P+1:M)*E(:,P+1:M)';Rinv=inv(Rxx);%compute MUSIC spatial spectrum using all noise eigenvectors%plot and compare with Minimum Variance spatial 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