%Demo done during Session 17 illustrating
%use of ESPRIT to estimate angle of arrival
%of each source incident on uniform linear
%array of antennas
clear all
clf
set(0,'defaultaxesfontsize',22);
M=12;
Nb=32;
P=3;
%angles=[25 80 135]*(pi/180);
angdeg=[80 90 115];
angles=angdeg*(pi/180);
powers=[10 10 10];
dlambda=0.5;
%generate random bits of information
poles=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(:,k)=exp(j*mu*(0:M-1)).';
end
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)*A(:,1)*b;
X=X+powers(2)*A(:,2)*(randn(1,Nb)+j*randn(1,Nb));
omega=3*pi/4;
X=X+powers(3)*A(:,3)*exp(j*omega*(0:Nb-1));
%add some noise
X=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 frequencies
polar(0,1,'.')
hold on
plot(real(poles),imag(poles),'kx','MarkerSize',18,'Linewidth',4);
plot(real(Phivec),imag(Phivec),'ro','MarkerSize',18,'Linewidth',4);
hold off
legend('True "poles"','ESPRIT eigenvalues')
angests=sort(acos(angle(Phivec)/pi)*(180/pi))
angdeg
% 
hold off
pause
Pn=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.'));
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