# Estimation of signal parameters via rotational invariance techniques

Estimation of signal parameters via rotational invariant techniques (ESPRIT), is a technique to determine the parameters of a mixture of sinusoids in background noise. This technique was first proposed for frequency estimation. However, with the introduction of phased-array systems in everyday technology, it is also used for angle of arrival estimations.

## One-dimensional ESPRIT

At instance t{\textstyle t}, the M{\textstyle M} (complex-valued) output signals (measurements) y$_m$[t]{\textstyle y_{m}[t]}, m=1, …, M{\textstyle m=1, \ldots , M}, of the system are related to the K{\textstyle K} (complex-valued) input signals x$_k$[t]{\textstyle x_{k}[t]}, k=1, …, K{\textstyle k=1, \ldots , K}, as$y_{m}[t]=\sum _{k=1}^{K}a_{m,k}x_{k}[t]+n_{m}[t],$where n$_m$[t]{\textstyle n_{m}[t]} denotes the noise added by the system. The one-dimensional form of ESPRIT can be applied if the weights have the form a$_{m, k}$=e$^{−j(m−1)ωk}${\textstyle a_{m, k}=e^{-j(m-1)\omega _{k}}}, whose phases are integer multiples of some radial frequency $\omega _{k}$. This frequency only depends on the index of the system's input, i.e. k{\textstyle k}. The goal of ESPRIT is to estimate $\omega _{k}$'s, given the outputs y$_m$[t]{\textstyle y_{m}[t]} and the number of input signals, K{\textstyle K}. Since the radial frequencies are the actual objectives, a$_{m, k}${\textstyle a_{m, k}} is denoted as a$_m$(ω$_k$){\textstyle a_{m}(\omega _{k})}.
Collating the weights a$_m$(ω$_k$){\textstyle a_{m}(\omega _{k})} as $\mathbf {a} (\omega _{k})=[\,1\,\ e^{-j\omega _{k}}\,\ e^{-j2\omega _{k}}\,\ ...\,\ e^{-j(M-1)\omega _{k}}\,]^{\top }$ and the M{\textstyle M} output signals at instance t{\textstyle t} as $\mathbf {y} [t]=[\,y_{1}[t]\,\ y_{2}[t]\,\ ...\,\ y_{M}[t]\,]^{\top }$, $\mathbf {y} [t]=\sum _{k=1}^{K}\mathbf {a} (\omega _{k})x_{k}[t]+\mathbf {n} [t],$where $\mathbf {n} [t]=[\,n_{1}[t]\,\ n_{2}[t]\,\ ...\,\ n_{M}[t]\,]^{\top }$. Further, when the weight vectors $\mathbf {a} (\omega _{1}),\ \mathbf {a} (\omega _{2}),\ ...,\ \mathbf {a} (\omega _{K})$ are put into a Vandermonde matrix $\mathbf {A} =[\,\mathbf {a} (\omega _{1})\,\ \mathbf {a} (\omega _{2})\,\ ...\,\ \mathbf {a} (\omega _{K})\,]$, and the $K$ inputs at instance t{\textstyle t} into a vector $\mathbf {x} [t]=[\,x_{1}[t]\,\ ...\,\ x_{K}[t]\,]^{\top }$, we can write$\mathbf {y} [t]=\mathbf {A} \,\mathbf {x} [t]+\mathbf {n} [t].$With several measurements at instances t=1, 2, …, T{\textstyle t=1, 2, \dots , T} and the notations Y=[y[1] y[2] … y[T]]{\textstyle \mathbf {Y} =[\, \mathbf {y} [1]\, \ \mathbf {y} [2]\, \ \dots \, \ \mathbf {y} [T]\, ]}, X=[x[1] x[2] … x[T]]{\textstyle \mathbf {X} =[\, \mathbf {x} [1]\, \ \mathbf {x} [2]\, \ \dots \, \ \mathbf {x} [T]\, ]} and N=[n[1] n[2] … n[T]]{\textstyle \mathbf {N} =[\, \mathbf {n} [1]\, \ \mathbf {n} [2]\, \ \dots \, \ \mathbf {n} [T]\, ]}, the model equation becomes$\mathbf {Y} =\mathbf {A} \mathbf {X} +\mathbf {N} .$

### Dividing into virtual sub-arrays

The weight vector a(ω$_k$){\textstyle \mathbf {a} (\omega _{k})} has the property that adjacent entries are related.$[\mathbf {a} (\omega _{k})]_{m+1}=e^{-j\omega _{k}}[\mathbf {a} (\omega _{k})]_{m}$For the whole vector $\mathbf {a} (\omega _{k})$, the equation introduces two selection matrices $\mathbf {J} _{1}$and $\mathbf {J} _{2}$: $\mathbf {J} _{1}=[\mathbf {I} _{M-1}\quad \mathbf {0} ]$ and $\mathbf {J} _{2}=[\mathbf {0} \quad \mathbf {I} _{M-1}]$. Here, $\mathbf {I} _{M-1}$ is an identity matrix of size (M−1){\textstyle (M-1)} and $\mathbf {0}$ is a vector of zeros.
The vector $\mathbf {J} _{1}\mathbf {a} (\omega _{k})$ $[\mathbf {J} _{2}\mathbf {a} (\omega _{k})]$ contains all elements of $\mathbf {a} (\omega _{k})$ except the last [first] one. Thus, $\mathbf {J} _{2}\mathbf {a} (\omega _{k})=e^{-j\omega _{k}}\mathbf {J} _{1}\mathbf {a} (\omega _{k})$ and$\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}}.$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 Y{\textstyle \mathbf {Y} } is given as$\mathbf {Y} =\mathbf {U} \mathbf {\Sigma } \mathbf {V} ^{\dagger }$where U∈C$^{M×M}${\textstyle \mathbf {U} \in \mathbb {C} ^{M\times M}} and V∈C$^{T×T}${\textstyle \mathbf {V} \in \mathbb {C} ^{T\times T}} are unitary matrices and Σ{\textstyle \mathbf {\Sigma } } is a diagonal matrix of size M×T{\textstyle M\times T}, that holds the singular values from the largest (top left) in descending order. The operator †{\textstyle \dagger } denotes the complex-conjugate transpose (Hermitian transpose).
Let us assume that T≥M{\textstyle T\geq M}. Notice that we have K{\textstyle K} input signals. If there was no noise, there would only be K{\textstyle K} non-zero singular values. We assume that the K{\textstyle K} largest singular values stem from these input signals and other singular values are presumed to stem from noise. The matrices in the SVD of Y{\textstyle \mathbf {Y} } can be partitioned into submatrices, where some submatrices correspond to the signal subspace and some correspond to the noise subspace.${\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}},$where U$_S$∈C$^{M×K}${\textstyle \mathbf {U} _{\mathrm {S} }\in \mathbb {C} ^{M\times K}} and V$_S$∈C$^{N×K}${\textstyle \mathbf {V} _{\mathrm {S} }\in \mathbb {C} ^{N\times K}} contain the first K{\textstyle K} columns of U{\textstyle \mathbf {U} } and V{\textstyle \mathbf {V} }, respectively and Σ$_S$∈C$^{K×K}${\textstyle \mathbf {\Sigma } _{\mathrm {S} }\in \mathbb {C} ^{K\times K}}is a diagonal matrix comprising the K{\textstyle K} largest singular values.
Thus, The SVD can be written as$\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 },$where U$_S${\textstyle \mathbf {U} _{\mathrm {S} }}, ⁣ V$_S${\textstyle \mathbf {V} _{\mathrm {S} }}, and Σ$_S${\textstyle \mathbf {\Sigma } _{\mathrm {S} }} represent the contribution of the input signal x$_k$[t]{\textstyle x_{k}[t]} to Y{\textstyle \mathbf {Y} }. We term U$_S${\textstyle \mathbf {U} _{\mathrm {S} }} the signal subspace. In contrast, U$_N${\textstyle \mathbf {U} _{\mathrm {N} }}, V$_N${\textstyle \mathbf {V} _{\mathrm {N} }}, and Σ$_N${\textstyle \mathbf {\Sigma } _{\mathrm {N} }} represent the contribution of noise n$_m$[t]{\textstyle n_{m}[t]} to Y{\textstyle \mathbf {Y} }.
Hence, from the system model, we can write AX=U$_S$Σ$_S$V$_S$$^†${\textstyle \mathbf {A} \mathbf {X} =\mathbf {U} _{\mathrm {S} }\mathbf {\Sigma } _{\mathrm {S} }\mathbf {V} _{\mathrm {S} }^{\dagger }} and N=U$_N$Σ$_N$V$_N$$^†${\textstyle \mathbf {N} =\mathbf {U} _{\mathrm {N} }\mathbf {\Sigma } _{\mathrm {N} }\mathbf {V} _{\mathrm {N} }^{\dagger }}. Also, from the former, we can write$\mathbf {U} _{\mathrm {S} }=\mathbf {A} \mathbf {F} ,$where ${\mathbf {F} }=\mathbf {X} \mathbf {V} _{\mathrm {S} }\mathbf {\Sigma } _{\mathrm {S} }^{-1}$. In the sequel, it is only important that there exists such an invertible matrix F{\textstyle \mathbf {F} } and its actual content will not be important.
Note: The signal subspace can also be extracted from the spectral decomposition of the auto-correlation matrix of the measurements, which is estimated as$\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 }.$

### Estimation of radial frequencies

We have established two expressions so far: J$_2$A=J$_1$AH{\textstyle \mathbf {J} _{2}\mathbf {A} =\mathbf {J} _{1}\mathbf {A} \mathbf {H} } and U$_S$=AF{\textstyle \mathbf {U} _{\mathrm {S} }=\mathbf {A} \mathbf {F} }. Now, ${\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}}$where $\mathbf {S} _{1}=\mathbf {J} _{1}\mathbf {U} _{\mathrm {S} }$ and $\mathbf {S} _{2}=\mathbf {J} _{2}\mathbf {U} _{\mathrm {S} }$ denote the truncated signal sub spaces, and $\mathbf {P} =\mathbf {F} ^{-1}\mathbf {H} \mathbf {F} .$The above equation has the form of an eigenvalue decomposition, and the phases of the eigenvalues in the diagonal matrix $\mathbf {H}$ are used to estimate the radial frequencies.
Thus, after solving for $\mathbf {P}$ in the relation $\mathbf {S} _{2}=\mathbf {S} _{1}\mathbf {P}$, we would find the eigenvalues λ$_1$, λ$_2$, …, λ$_K${\textstyle \lambda _{1}, \ \lambda _{2}, \ \ldots , \ \lambda _{K}} of $\mathbf {P}$, where $\lambda _{k}=\alpha _{k}\mathrm {e} ^{j\omega _{k}}$, and the radial frequencies $\omega _{1},\ \omega _{2},\ \ldots ,\ \omega _{K}$ are estimated as the phases (argument) of the eigenvalues.
Remark: In general, $\mathbf {S} _{1}$ is not invertible. One can use the least squares estimate ${\mathbf {P} }=(\mathbf {S} _{1}^{\dagger }{\mathbf {S} _{1}})^{-1}\mathbf {S} _{1}^{\dagger }{\mathbf {S} _{2}}$. An alternative would be the total least squares estimate.

### Algorithm summary

Input: Measurements Y:=[y[1] y[2] … y[T]]{\textstyle \mathbf {Y} :=[\, \mathbf {y} [1]\, \ \mathbf {y} [2]\, \ \dots \, \ \mathbf {y} [T]\, ]}, the number of input signals K{\textstyle K} (estimate if not already known).
* Compute the singular value decomposition (SVD) of Y=UΣV$^†${\textstyle \mathbf {Y} =\mathbf {U} \mathbf {\Sigma } \mathbf {V} ^{\dagger }} and extract the signal subspace U$_S$∈C$^{M×K}${\textstyle \mathbf {U} _{\mathrm {S} }\in \mathbb {C} ^{M\times K}} as the first K{\textstyle K} columns of U{\textstyle \mathbf {U} }.
* Compute S$_1$=J$_1$U$_S${\textstyle \mathbf {S} _{1}=\mathbf {J} _{1}\mathbf {U} _{\mathrm {S} }} and S$_2$=J$_2$U$_S${\textstyle \mathbf {S} _{2}=\mathbf {J} _{2}\mathbf {U} _{\mathrm {S} }}, where $\mathbf {J} _{1}=[\mathbf {I} _{M-1}\quad \mathbf {0} ]$ and $\mathbf {J} _{2}=[\mathbf {0} \quad \mathbf {I} _{M-1}]$.
* Solve for P{\textstyle \mathbf {P} } in S$_2$=S$_1$P{\textstyle \mathbf {S} _{2}=\mathbf {S} _{1}\mathbf {P} } (see the remark above).
* Compute the eigenvalues λ$_1$, λ$_2$, …, λ$_K${\textstyle \lambda _{1}, \lambda _{2}, \ldots , \lambda _{K}} of P{\textstyle \mathbf {P} }.
* The phases of the eigenvalues λ$_k$=α$_k$e$^{jωk}${\textstyle \lambda _{k}=\alpha _{k}\mathrm {e} ^{j\omega _{k}}} provide the radial frequencies ω$_k${\textstyle \omega _{k}}, i.e., ω$_k$=arg⁡λ$_k${\textstyle \omega _{k}=\arg \lambda _{k}} .