% espritapp.m
% a simple application of the ESPRIT 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, and the signal subspace
[U S V]=svd(R);
Neig=4;			% the 4 first left eigenvectors
Ssub=U(:,1:Neig);	% span the signal subspace.

% now for the ESPRIT version
S1=Ssub(1:M,:);
S2=Ssub(2:M+1,:);
poles=eig(S1\S2);

% plot the poles
figure(1), clf, axis equal, hold on
plot(real(poles),imag(poles),'rp')
t=0:1000;
t=2*pi/1000*t;
plot(cos(t),sin(t),'b',[-1.2 1.2], [0 0],'k',[0 0],[-1.2 1.2],'k')
title('Detected frequencies')
% Check the angles, should be 5 (23) /200 *2 *pi

% 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 