Estimation of signal parameters via rotational invariance techniques

In estimation theory, estimation of signal parameters via rotational invariant techniques (ESPRIT) is a technique to determine parameters of a mixture of sinusoids in a background noise. This technique is first proposed for frequency estimation, however, with the introduction of phased-array systems in daily use technology, it is also used for Angle of arrival estimations as well.
General description.
Let us define a signal vector as,
formula_1
where formula_2 represents the radial frequency of kth sinusoid. Let us construct a Vandermonde matrix for K number of sinusoids, such as
formula_3
Let us divide the matrix A into two sets, such as
formula_4
and
formula_5
where formula_6 is an identity matrix of size (m-1)-by-(m-1). It is clear that, formula_7 contains the first (m-1) rows of A, while formula_8 contains the last (m-1) rows of A.
Now, let us write the following relation,
formula_9
Here, H is a diagonal matrix, where its diagonal elements can be written in a vector, such as,
formula_10
In other words, the diagonal elements of H, are the complex exponentials with the radial frequencies of the set formula_11. Here, it is clear that H applies a rotation to the matrix formula_7. ESPRIT exploits similar rotations from the covariance matrix of the measured data.
To understand the algorithm itself, let us denote R as the covariance matrix of the measured data. By computing the eigenvalue decomposition of R (via algorithms like singular value decomposition), the following can be written,
formula_13
where E is a diagonal matrix that contains the eigenvalues of R, in a decreasing order. Here, by finding the eigenvalues that are higher than the variance of the noise, we can separate the orthonormal eigenvectors from U, that correspond to these eigenvalues. This can be noted as formula_14 where we kept only the first K columns. 
As similar before, we can make the following separation on S,
formula_15 and formula_16.
Moreover, there exists a relation between S and A such as formula_17, where the content of the matrix F is known, but irrelevant for the current subject. We can derive the following relations,
formula_18
(where we made use of formula_19 and formula_20).
It is clear that the matrix P contains rotational information with respect to the frequency contents, such that the rotation on the first set of orthonormal eigenvectors yield to the second set. Moreover, the eigenvalues of P are equal to the diagonal elements of H. 
Therefore, by solving the following equation for P,
formula_21 
we can estimate the frequency content. To achieve this, the above equation can be solved with pseudo inverse (via Least squares) method.
To do so, formula_22 can be written. 
Finally, by finding the angles of the eigenvalues of P, one can estimate the set formula_11.
Algorithm example.
A pseudo code is given below for the implementation of ESPRIT algorithm.
 function esprit("y", "model_order", "number_of_sources"):
 "m" = model_order
 "n" = number_of_sources
 create covariance matrix R, from the noisy measurements y. Size of R will be (m-by-m).
 compute the svd of R
 [U, E, V] = svd(R)
 obtain the orthonormal eigenvectors corresponding to the sources
 "S" = U(:, 1:n) 
 split the orthonormal eigenvectors in two
 S1 = S(1:m-1, :) and S2 = S(2:m, :)
 compute P via LS (MATLAB's backslash operator)
 "P" = S1\S2 
 find the angles of the eigenvalues of P
 "w" = angle(eig(P)) / (2*pi*elspacing)
 doa=asind(w) %return the doa angle by taking the arcsin in degrees 
 return 'doa"