View source for Estimation of signal parameters via rotational invariance techniques

{{Short description|Signal processing method}}
{{multiple issues|
{{Expert needed|statistics
| date = November 2023
| reason = article reads like an overly detailed proof, and lacks encyclopedic tone. It needs rewriting as per [[MOS:MATH#TONE]]
}}
{{Textbook|date=November 2023}}
{{Copy edit|date=October 2023}}}}
[[File:ESPRIT 2D.gif|thumb|350px|Example of separation into subarrays (2D ESPRIT)]]

'''Estimation of signal parameters via rotational invariant techniques (ESPRIT),''' is a technique to determine the parameters of a mixture of [[Sine wave|sinusoid]]s in background noise. This technique was first proposed for frequency estimation.<ref>{{Citation | last1=Paulraj | first1=A. | last2=Roy | first2=R. | last3=Kailath | first3=T. | title=Nineteenth Asilomar Conference on Circuits, Systems and Computers | isbn=978-0-8186-0729-5 | year=1985 | chapter=Estimation Of Signal Parameters Via Rotational Invariance Techniques - Esprit | doi=10.1109/ACSSC.1985.671426 | pages=83–89| s2cid=2293566 }}</ref> However, with the introduction of [[phased-array]] systems in everyday technology, it is also used for [[angle of arrival]] estimations.<ref>Volodymyr Vasylyshyn. The direction of arrival estimation using ESPRIT with sparse arrays.// Proc. 2009 European Radar Conference (EuRAD). – 30 Sept.-2 Oct. 2009. - Pp. 246 - 249. - [https://ieeexplore.ieee.org/abstract/document/5307000]</ref>

== One-dimensional ESPRIT ==
At instance <math display="inline">t
</math>, the <math display="inline">M
</math> (complex-valued) output signals (measurements) <math display="inline">y_m[t]</math>, <math display="inline">m = 1, \ldots, M
</math>, of the system are related to the <math display="inline">K
</math> (complex-valued) input signals <math display="inline">x_k[t] 
</math>, <math display="inline">k = 1, \ldots, K
</math>, as<math display="block">y_m[t] = \sum_{k=1}^K a_{m,k} x_k[t] + n_m[t],</math>where <math display="inline">n_m[t] 
</math> denotes the noise added by the system. The one-dimensional form of ESPRIT can be applied if the weights have the form <math display="inline">a_{m,k} = e^{-j(m-1)\omega_k}
</math>, whose phases are integer multiples of some [[Angular frequency|radial frequency]] <math>\omega_k
</math>. This frequency only depends on the index of the system's input, i.e. <math display="inline">k
</math>. The goal of ESPRIT is to estimate <math>\omega_k
</math>'s, given the outputs <math display="inline"> y_m[t] 
</math> and the number of input signals, <math display="inline">K
</math>. Since the radial frequencies are the actual objectives, <math display="inline">a_{m,k}
</math> is denoted as <math display="inline">a_m(\omega_k)
</math>.   

Collating the weights <math display="inline">a_m(\omega_k)
</math> as <math>\mathbf a(\omega_k) = [\,1 \,\ e^{-j\omega_k} \,\ e^{-j2\omega_k} \,\ ... \,\ e^{-j(M-1)\omega_k}\,]^\top</math> and the <math display="inline">M
</math> output signals at instance <math display="inline">t
</math> as <math>\mathbf y[t] = [\,y_1[t]\,\ y_2[t] \,\ ... \,\ y_M[t]\,]^\top</math>, <math display="block">\mathbf y[t] = \sum_{k=1}^K \mathbf a(\omega_k) x_k[t] + \mathbf{n}[t],</math>where <math>\mathbf n[t] = [\,n_1[t]\,\ n_2[t] \,\ ... \,\ n_M[t]\,]^\top</math>. Further, when the weight vectors <math>\mathbf a(\omega_1), \ \mathbf a(\omega_2),\ ..., \  \mathbf a(\omega_K)
</math> are put into a [[Vandermonde matrix]] <math>\mathbf A = [\,\mathbf a(\omega_1) \,\ \mathbf a(\omega_2) \,\ ... \,\ \mathbf a(\omega_K)\,]</math>, and the <math>K</math> inputs at instance <math display="inline">t
</math> into a vector <math>\mathbf x[t] = [\,x_1[t] \,\ ... \,\ x_K[t]\,]^\top</math>, we can write<math display="block"> \mathbf y[t] = \mathbf A\, \mathbf x[t]+\mathbf n[t].</math>With several measurements at instances <math display="inline">t = 1, 2, \dots, T
</math> and the notations <math display="inline">\mathbf{Y} = [\,\mathbf{y}[1] \,\ \mathbf{y}[2] \,\ \dots \,\ \mathbf{y}[ T ]\,]</math>, <math display="inline">\mathbf{X} = [\,\mathbf{x}[1] \,\ \mathbf{x}[2] \,\ \dots \,\ \mathbf{x}[ T ]\,]</math> and <math display="inline">\mathbf{N} = [\, \mathbf{n}[1] \,\ \mathbf{n}[2] \,\ \dots \,\ \mathbf{n}[ T ]\,]</math>, the model equation becomes<math display="block">\mathbf{Y} = \mathbf{A} \mathbf{X} + \mathbf{N}.</math>

=== Dividing into virtual sub-arrays ===

[[File:Max overlapping subarrays.png|thumb|350px|Maximum overlapping of two sub-arrays (''N'' denotes number of sensors in the array, ''m'' is the number of sensors in each sub-array, and <math>J_1</math> and <math>J_2</math> are selection matrices)]]

The weight vector <math display="inline">\mathbf a(\omega_k)
</math> has the property that adjacent entries are related.<math display="block">[\mathbf{a}(\omega_k)]_{m+1} = e^{-j\omega_k} [\mathbf{a}(\omega_k)]_{m}
</math>For the whole vector <math>\mathbf a(\omega_k)
</math>, the equation introduces two selection matrices <math>\mathbf{J}_1
</math>and <math>\mathbf{J}_2
</math>: <math>\mathbf J_1 = [\mathbf I_{M-1} \quad \mathbf 0]</math> and <math>\mathbf J_2 = [\mathbf 0 \quad \mathbf I_{M-1}]</math>. Here, <math> \mathbf I_{M-1} </math> is an identity matrix of size <math display="inline">(M-1)</math> and <math>\mathbf 0</math> is a vector of zeros.  

The vector <math>\mathbf J_1 \mathbf a(\omega_k)
</math> <math>[\mathbf J_2 \mathbf a(\omega_k)]
</math> contains all elements of <math>\mathbf a(\omega_k)
</math> except the last [first] one. Thus, <math>\mathbf J_2 \mathbf a(\omega_k) = e^{-j\omega_k} \mathbf J_1 \mathbf a(\omega_k)</math> and<math display="block">\mathbf J_2 \mathbf A = \mathbf J_1 \mathbf A \mathbf H
,\quad\text{where}\quad
{\mathbf {H}} := \begin{bmatrix} e^{-j\omega_1} &  \\  & e^{-j\omega_2} \\ & & \ddots \\ & & & e^{-j\omega_K} \end{bmatrix}
.</math>The above relation is the first major observation required for ESPRIT. The second major observation concerns the signal subspace that can be computed from the output signals. 

=== Signal subspace ===
The [[singular value decomposition]] (SVD) of <math display="inline">\mathbf Y
</math> is given as<math display="block">\mathbf{Y} = \mathbf U \mathbf \Sigma \mathbf V^\dagger</math>where <math display="inline">\mathbf U \in\mathbb C^{M\times M}
</math> and <math display="inline">\mathbf V \in \mathbb C^{T\times T}
</math> are unitary matrices and <math display="inline">\mathbf \Sigma
</math> is a diagonal matrix of size <math display="inline">M \times T</math>, that holds the singular values from the largest (top left) in descending order. The operator <math display="inline">\dagger
</math> denotes the complex-conjugate transpose (Hermitian transpose). 

Let us assume that <math display="inline">T \geq M</math>. Notice that we have <math display="inline">K</math> input signals. If there was no noise, there would only be <math display="inline">K</math> non-zero singular values. We assume that the <math display="inline">K</math> largest singular values stem from these input signals and other singular values are presumed to stem from noise. The matrices in the SVD of <math display="inline">\mathbf Y</math> can be partitioned into submatrices, where some submatrices correspond to the signal subspace and some correspond to the noise subspace.<math display="block"> \begin{aligned}
\mathbf{U} = 
\begin{bmatrix}
\mathbf{U}_\mathrm{S} &
\mathbf{U}_\mathrm{N}
\end{bmatrix},
& & 
\mathbf{\Sigma} = 
\begin{bmatrix}
\mathbf{\Sigma}_\mathrm{S} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{\Sigma}_\mathrm{N} & \mathbf{0}
\end{bmatrix},
& & 
\mathbf{V} = 
\begin{bmatrix}
\mathbf{V}_\mathrm{S} &
\mathbf{V}_\mathrm{N} &
\mathbf{V}_\mathrm{0}
\end{bmatrix}
\end{aligned},</math>where <math display="inline">\mathbf{U}_\mathrm{S} \in \mathbb{C}^{M \times K}
</math> and <math display="inline">\mathbf{V}_\mathrm{S} \in \mathbb{C}^{N \times K}
</math> contain the first <math display="inline">K</math> columns of <math display="inline">\mathbf U
</math> and <math display="inline">\mathbf V
</math>, respectively and <math display="inline">\mathbf{\Sigma}_\mathrm{S} \in \mathbb{C}^{K \times K}
</math>is a diagonal matrix comprising the <math display="inline">K</math> largest singular values. 

Thus, The SVD can be written as<math display="block"> \mathbf{Y} =  \mathbf U_\mathrm{S} \mathbf \Sigma_\mathrm{S} \mathbf V_\mathrm{S}^\dagger 
+ \mathbf U_\mathrm{N} \mathbf \Sigma_\mathrm{N} \mathbf V_\mathrm{N}^\dagger 
,</math>where <math display="inline">\mathbf{U}_\mathrm{S}
</math>, ⁣<math display="inline">\mathbf{V}_\mathrm{S}
</math>, and <math display="inline">\mathbf{\Sigma}_\mathrm{S}
</math> represent the contribution of the input signal <math display="inline">x_k[t]
</math> to <math display="inline">\mathbf Y
</math>. We term <math display="inline">\mathbf{U}_\mathrm{S}
</math> the signal subspace. In contrast, <math display="inline">\mathbf{U}_\mathrm{N}
</math>, <math display="inline">\mathbf{V}_\mathrm{ N }
</math>, and <math display="inline">\mathbf{\Sigma}_\mathrm{N}
</math> represent the contribution of noise <math display="inline">n_m[ t ]
</math> to <math display="inline">\mathbf Y
</math>. 

Hence, from the system model, we can write <math display="inline">\mathbf{A} \mathbf{X} = 
\mathbf U_\mathrm{S} \mathbf \Sigma_\mathrm{S} \mathbf V_\mathrm{S}^\dagger
</math> and <math display="inline">\mathbf{N} = 
\mathbf U_\mathrm{N} \mathbf\Sigma_\mathrm{N} \mathbf V_\mathrm{N}^\dagger
</math>. Also, from the former, we can write<math display="block"> \mathbf U_\mathrm{S}
= \mathbf{A} \mathbf{F},
</math>where <math> {\mathbf F } = \mathbf{X} \mathbf V_\mathrm{S} \mathbf \Sigma_\mathrm{S}^{-1}
</math>. In the sequel, it is only important that there exists such an invertible matrix <math display="inline">\mathbf F
</math> and its actual content will not be important.

'''Note:''' The signal subspace can also be extracted from the [[Spectral decomposition (Matrix)|spectral decomposition]] of the [[Autocorrelation matrix|auto-correlation matrix]] of the measurements, which is estimated as<math display="block">\mathbf{R}_\mathrm{YY} 
=  
\frac{1}{T}
\sum_{t=1}^T \mathbf{y}[t] \mathbf{y}[t] ^\dagger
=\frac{1}{T}
\mathbf{Y} \mathbf{Y}^\dagger
= 
\frac{1}{T}
\mathbf U 
{\mathbf \Sigma \mathbf \Sigma^\dagger}
\mathbf U^\dagger
= 
\frac{1}{T}
\mathbf U_\mathrm{S} \mathbf \Sigma_\mathrm{S}^2 \mathbf U_\mathrm{S}^\dagger
+
\frac{1}{T}
\mathbf U_\mathrm{N} \mathbf \Sigma_\mathrm{N}^2 \mathbf U_\mathrm{N}^\dagger.</math>

=== Estimation of radial frequencies ===
We have established two expressions so far: <math display="inline">\mathbf J_2 \mathbf A = \mathbf J_1 \mathbf A \mathbf H
</math> and <math display="inline">\mathbf U_\mathrm{S} = \mathbf{A} \mathbf{F}
</math>. Now, <math display="block"> \begin{aligned}
\mathbf J_2 \mathbf A 
= \mathbf J_1 \mathbf A \mathbf H
\implies
\mathbf J_2 (\mathbf U_\mathrm{S} \mathbf{F}^{-1})
= \mathbf J_1 (\mathbf U_\mathrm{S} \mathbf{F}^{-1}) \mathbf H
\implies
\mathbf S_2 = \mathbf S_1 \mathbf{P},
\end{aligned}</math>where <math> \mathbf{S}_1 = \mathbf J_1 \mathbf U_\mathrm{S} </math> and <math> \mathbf{S}_2 = \mathbf J_2 \mathbf U_\mathrm{S}</math> denote the truncated signal sub spaces, and <math display="block"> \mathbf{P} = \mathbf{F}^{-1} \mathbf H \mathbf{F}.</math>The above equation has the form of an [[eigenvalue decomposition]], and the phases of the eigenvalues in the diagonal matrix <math> \mathbf H</math> are used to estimate the radial frequencies.  

Thus, after solving for <math>\mathbf P</math> in the relation <math> \mathbf S_2 = \mathbf S_1  \mathbf{P}</math>, we would find the eigenvalues <math display="inline">\lambda_1,\ \lambda_2,\ \ldots,\ \lambda_K
</math> of <math> \mathbf P</math>, where <math>\lambda_k = \alpha_k \mathrm{e}^{j \omega_k}</math>, and the radial frequencies <math>\omega_1,\ \omega_2,\ \ldots,\ \omega_K
</math> are estimated as the phases (argument) of the eigenvalues. 

'''Remark:''' In general, <math> \mathbf S_1</math> is not invertible. One can use the [[least squares]] estimate <math> {\mathbf P} = (\mathbf S_1^\dagger {\mathbf S_1})^{-1} \mathbf S_1^\dagger {\mathbf S_2}</math>. An alternative would be the [[total least squares]] estimate.

=== Algorithm summary ===
'''Input:''' Measurements <math display="inline">\mathbf{Y} := [\,\mathbf{y}[1] \,\ \mathbf{y}[2] \,\ \dots \,\ \mathbf{y}[ T ]\,]
</math>, the number of input signals <math display="inline">K
</math> (estimate if not already known).
# Compute the singular value decomposition (SVD) of <math display="inline">\mathbf{Y} = \mathbf U \mathbf \Sigma \mathbf V^\dagger</math> and extract the signal subspace <math display="inline">\mathbf{U}_\mathrm{S} \in \mathbb{C}^{M \times K}
</math> as the first <math display="inline">K</math> columns of <math display="inline">\mathbf U </math>.
# Compute <math display="inline">\mathbf{S}_1 = \mathbf J_1 \mathbf U_\mathrm{S} </math> and <math display="inline">\mathbf{S}_2 = \mathbf J_2 \mathbf U_\mathrm{S}
</math>, where <math> \mathbf J_1 = [\mathbf I_{M-1} \quad \mathbf 0] </math> and <math> \mathbf J_2 = [\mathbf 0 \quad \mathbf I_{M-1}] </math>.
# Solve for <math display="inline">\mathbf{P}
</math> in <math display="inline">\mathbf S_2 = \mathbf S_1  \mathbf{P}
</math> (see the remark above).
# Compute the eigenvalues <math display="inline">\lambda_1, \lambda_2, \ldots, \lambda_K
</math> of <math display="inline">\mathbf{P}
</math>.
# The phases of the eigenvalues <math display="inline">\lambda_k = \alpha_k \mathrm{e}^{j \omega_k}
</math> provide the radial frequencies <math display="inline">\omega_k
</math>, i.e., <math display="inline">\omega_k = \arg \lambda_k
</math> .

=== Notes ===

==== Choice of selection matrices ====
In the derivation above, the selection matrices <math> \mathbf J_1 = [\mathbf I_{M-1} \quad \mathbf 0] </math> and <math> \mathbf J_2 = [\mathbf 0 \quad \mathbf I_{M-1}] </math> were used. However, any appropriate matrices <math display="inline"> \mathbf J_1 </math> and <math> \mathbf{J}_2 </math> may be used as long as the rotational invariance <math display="inline"> ( </math>i.e., <math> \mathbf J_2 \mathbf a(\omega_k) = e^{-j\omega_k} \mathbf J_1 \mathbf a(\omega_k) </math><math display="inline"> ) </math>, or some generalization of it (see below) holds; accordingly, the matrices <math display="inline"> \mathbf J_1 \mathbf A </math> and <math display="inline"> \mathbf J_2 \mathbf A </math> may contain any rows of <math display="inline"> \mathbf A </math>.

==== Generalized rotational invariance ====
The rotational invariance used in the derivation may be generalized. So far, the matrix <math> \mathbf H</math> has been defined to be a diagonal matrix that stores the sought-after complex exponentials on its main diagonal. However, <math> \mathbf H</math> may also exhibit some other structure.<ref>{{Cite journal |last1=Hu |first1=Anzhong |last2=Lv |first2=Tiejun |last3=Gao |first3=Hui |last4=Zhang |first4=Zhang |last5=Yang |first5=Shaoshi |date=2014 |title=An ESPRIT-Based Approach for 2-D Localization of Incoherently Distributed Sources in Massive MIMO Systems |url=https://ieeexplore.ieee.org/document/6777542 |journal=IEEE Journal of Selected Topics in Signal Processing |volume=8 |issue=5 |pages=996–1011 |doi=10.1109/JSTSP.2014.2313409 |arxiv=1403.5352 |bibcode=2014ISTSP...8..996H |s2cid=11664051 |issn=1932-4553}}</ref> For instance, it may be an upper triangular matrix. In this case, <math display="inline"> \mathbf{P} = \mathbf{F}^{-1} \mathbf H \mathbf{F}</math>constitutes a '''triangularization''' of <math> \mathbf P</math>.

==See also==
*[[Multiple signal classification]]
*[[Generalized pencil-of-function method]]
*[[Independent component analysis]]

==References==
{{reflist}}

==Further reading==
* {{Citation | last1=Paulraj | first1=A. | last2=Roy | first2=R. | last3=Kailath | first3=T. | title=Nineteenth Asilomar Conference on Circuits, Systems and Computers | isbn=978-0-8186-0729-5 | year=1985 | chapter=Estimation Of Signal Parameters Via Rotational Invariance Techniques - Esprit | doi=10.1109/ACSSC.1985.671426 | pages=83–89| s2cid=2293566 }}.
* {{Cite journal | last1=Roy | first1=R. | last2=Kailath | first2=T. | title=Esprit - Estimation Of Signal Parameters Via Rotational Invariance Techniques | year=1989 | journal=IEEE Transactions on Acoustics, Speech, and Signal Processing | url=http://www.vtvt.ece.vt.edu/research/references/uwb/ranging_mobile_location/esprit.pdf | volume=37 | issue=7 | pages=984–995 | doi=10.1109/29.32276 | s2cid=14254482 | access-date=2011-07-25 | archive-date=2020-09-26 | archive-url=https://web.archive.org/web/20200926212409/https://www.vtvt.ece.vt.edu/research/references/uwb/ranging_mobile_location/esprit.pdf | url-status=dead }}.
*{{cite journal|first1=A. M.|last1= Ibrahim|first2= M. I.|last2= Marei|first3= S. F.|last3= Mekhamer|first4= M. M.|last4= Mansour | journal=Electric Power Components and Systems 
|volume= 39| issue= 1|year= 2011 
 |title=An Artificial Neural Network Based Protection Approach Using Total Least Square Estimation of Signal Parameters via the Rotational Invariance Technique for Flexible AC Transmission System Compensated Transmission Lines|doi= 10.1080/15325008.2010.513363 |pages= 64–79|s2cid= 109581436}}
* Haardt, M., Zoltowski, M. D., Mathews, C. P., & Nossek, J. (1995, May). 2D unitary ESPRIT for efficient 2D parameter estimation. In icassp (pp.&nbsp;2096-2099). IEEE.

[[Category:Signal estimation]]
[[Category:Trigonometry]]
[[Category:Wave mechanics]]