Digest of ESPRIT Algorithm

For most folks, ESPRIT means a fashion branch (http://www.esprit.eu/). But for nerds, ESPRIT stands for Estimation of Signal Parameter via Rotational Invariance Technique. It is a signal processing algorithm used to process things like signal coming from array. ESPRIT is often mentioned together with MUSIC algorithm. Again, MUSIC is not the real music played by Mozart but stands for MUltiple SIgnal Classification.

Recently I spent some time to understand the details of ESPRIT algorithm and here is what I have learnt. The biggest difference between ESPRIT and MUSIC is that ESPRIT has dual receivers. The diagram below is from a tutorial found in Internet (http://www.girdsystems.com/pdf/GIRD_Systems_Intro_to_MUSIC_ESPRIT.pdf). As shown here, to make ESPRIT algorithm possible, there are many pairs of sensors. The angle of arrive signal is obtained by comparing these two received signals.

The best way for understanding such algorithm is to run a Matlab script. The website of a book named "Spectral Analysis of Signals" provide a well written example of Matlab code for ESPRIT. The code is as below:

function w=esprit(y,n,m)
%
% The ESPRIT method for frequency estimation.
%
%  w=esprit(y,n,m);
%
%      y  ->  the data vector
%      n  ->  the model order
%      m  ->  the order of the covariance matrix in (4.5.14)
%      w  <-  the frequency estimates
%

% Copyright 1996 by R. Moses

y=y(:);
N=length(y);                       % data length

% compute the sample covariance matrix
R=zeros(m,m);
for i = m : N,
   R=R+y(i:-1:i-m+1)*y(i:-1:i-m+1)'/N;
end

% to use the forward-backward approach, uncomment the next line
% R=(R+fliplr(eye(m))*R.'*fliplr(eye(m)))/2;

% get the eigendecomposition of R; use svd because it sorts eigenvalues
[U,D,V]=svd(R);
S=U(:,1:n);

phi = S(1:m-1,:)\S(2:m,:);

w=-angle(eig(phi));
return

R in the code is a m x m matrix. The algorithm divides S(1:m, :) to two parts of S(1:m-1, :) and S(2:m, :). It is like dividing a large array to two smaller arrays with overlap of m-2 elements. I think this can be further fine tuned. For example, S(1:m, :) can be divided to S(1:m-2, :) and S(3:m, :) with one less overlapping element. "S(1:m-1,:)\S(2:m,:)" is to find least square solution of S(1:m-1,:)*X = S(2:m,:). By this method, it finds the delta between received signal of these two slightly different arrays.