musicdemo.m

matlab code for direction finding algorithm of MUSIC

%Demo done during Session 16 illustrating
%that the DSFT of each noise eigenvector
%exhibits a null at the spatial frequency
%associated with each source
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));
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)
% 
%compute DSFT of each noise eigenvector and
%to show that each has a null 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;
eigenspec=fftshift(fft(E(:,nk),1024));
fmag=20*log10(abs(eigenspec));
fmag=fmag-max(fmag);
nc=nc+1;
plot(theta,fmag,s(nc),'LineWidth',4);
axis([0 180 -60 0]);
hold on
angdeg=angles*(180/pi);
for p=1:P;
plot([angdeg(1,p) angdeg(1,p)],[-60 0],'k','LineWidth',3);
end
title('DSFT of Noise Eigenvector'),...
xlabel('Angle (Degrees)'), ylabel('Magnitude'), grid
pause
end
hold off
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