% musicapp.m
% a simple application of the MUSIC algorithm
clear

N=200;
sigma=1;
t=(1:N)';
s=sin(5*2*pi/N*t)+sin(23*2*pi/N*t);	% a rank 4 signal
data=s+sigma*randn(N,1);

% estimate the autocorrelation matrix
M=19;
R=zeros(M+1,M+1);
for n=1:N-M
	R=R+data(n:n+M)*data(n:n+M)';
end
R=R/(N-M);

% calculate the rangespace of R
[U S V]=svd(R);
Neig=4;			% the 4 first left eigenvectors
Ssub=U(:,1:Neig);	% span the signal subspace.

% calculate the pseudospectrum
m=1;
steer=0:M;
for angle=-pi:pi/100:pi
	aflip=exp(-j*angle*steer);	% the transposed steering vector for argument angle
	dum=(norm(aflip*Ssub))^2/(M+1);	% norm(aflip)^2=M+1
	P(m)=1/(1-dum);			% this screws up if you remove the noise
	m=m+1;
end

% plot the data and the pseudospectrum.
figure(1), clf, hold on
plot(data)
title('Two sinusoids in white noise')
figure(2), clf, hold on
plot([-pi:pi/100:pi],P)
title('MUSIC pseudospectrum')

% There are a lot of experiments to be carried out.
%	Vary N
%	Vary sigma - too much noise and you may miss the signal subspace. 
%	Vary Neig, we do not always know the rank of the signal subspace
%	Vary M, it affects resolution