Report on Two-Step Knowledge-Aided Iterative ESPRIT Algorithm

TL;DR

This paper introduces a two-step iterative ESPRIT algorithm that leverages prior knowledge and covariance matrix structure to enhance the accuracy of DOA estimation in signal processing.

Contribution

It presents a novel two-step knowledge-aided iterative ESPRIT method that improves DOA estimation accuracy by incorporating prior information and covariance matrix structure.

Findings

Improved DOA estimation accuracy over prior methods

Effective utilization of prior knowledge in covariance matrix estimation

Simulation results demonstrate superior performance

Abstract

In this work, we propose a subspace-based algorithm for direction-of-arrival (DOA) estimation, referred to as two-step knowledge-aided iterative estimation of signal parameters via rotational invariance techniques (ESPRIT) method (Two-Step KAI-ESPRIT), which achieves more accurate estimates than those of prior art. The proposed Two-Step KAI-ESPRIT improves the estimation of the covariance matrix of the input data by incorporating prior knowledge of signals and by exploiting knowledge of the structure of the covariance matrix and its perturbation terms. Simulation results illustrate the improvement achieved by the proposed method.

Tables

Table 1. TABLE I: Proposed Two-Step KAI-ESPRIT Algorithm

Inputs:

$M$,  $d$,  $\lambda$,  $N$,  $P$

Received vectors $𝒙\left(1\right)$,  $𝒙\left(2\right)$,$\cdots$, $𝒙\left(N\right)$

Prior knowledge $\to$ known DOAs:  ${\theta }_1$,  ${\theta }_2$, $\cdots$,  ${\theta }_q$

Outputs:

Estimates ${\hat{\theta }}_{q+1}^{\left(2\right)}$,  ${\hat{\theta }}_{q+2}^{\left(2\right)}$,$\cdots$,  ${\hat{\theta }}_P^{\left(2\right)}$

First step:

$\widehat{𝐑}=\frac{1}{N}\sum _{i=1}^N𝒙\left(i\right){𝒙}^H\left(i\right)$

$\{\hat{\theta }_1^{}{}_{}{}^{\left(1\right)},\hat{\theta }_2^{}{}_{}{}^{\left(1\right)},\cdots ,\hat{\theta }_P^{}{}_{}{}^{\left(1\right)}\}\underleftarrow{ESPRIT}$ $(\widehat{𝐑},P,d,\lambda )$

Second step:

${\widehat{𝐀}}^{\left(1\right)}=[𝐚\left({\hat{\theta }}_1^{\left(1\right)}\right),𝐚\left({\hat{\theta }}_2^{\left(1\right)}\right),\cdots ,𝐚\left({\hat{\theta }}_P^{\left(1\right)}\right)]$

compute for n = 1

${\widehat{𝐐}}_A^{\left(n\right)}={\widehat{𝐀}}^{\left(n\right)}{\left({\widehat{𝐀}}^{\left(n\right)H}{\widehat{𝐀}}^{\left(n\right)}\right)}^{-1}{\widehat{𝐀}}^{\left(n\right)H}$       (1)

${\widehat{𝐐}}_A^{\left(n\right)⟂}={𝐈}_M-{\widehat{𝐐}}_A^{\left(n\right)}$                           (2)

${𝐕}^{\left(n\right)}={\widehat{𝐐}}_A^{\left(n\right)}\widehat{𝐑}{\widehat{𝐐}}_A^{\left(n\right)⟂}$

for $\mu =\text{0}:\text{increment\hspace{0.33em}:\hspace{0.33em}1}$

compute   ${\widehat{𝐑}}^{\left(n+1\right)}=\widehat{𝐑}-\mu \left({𝐕}^{\left(n\right)}+{𝐕}^{\left(n\right)H}\right)$,   n = 1

$\{\hat{\theta }_{q+1}^{}{}_{}{}^{\left(2\right)},\hat{\theta }_{q+2}^{}{}_{}{}^{\left(2\right)},\cdots ,\hat{\theta }_P^{}{}_{}{}^{\left(2\right)}\}\underleftarrow{ESPRIT}$ $({\widehat{𝐑}}^{\left(2\right)},P,d,\lambda )$

${\widehat{𝐀}}^{\left(2\right)}=[𝐚\left({\theta }_1\right),\dots ,𝐚\left({\theta }_q\right),𝐚\left({\hat{\theta }}_{q+1}^{\left(2\right)}\right),\dots ,𝐚\left({\hat{\theta }}_P^{\left(2\right)}\right)]$

repeat (1) and (2),  n = 2

$U\left(\mu \right)=\ln \det \left({\widehat{𝐐}}_A^{\left(2\right)}\widehat{𝐑}{\widehat{𝐐}}_A^{\left(2\right)}+\frac{\mathrm{Trace}\left\{{\widehat{𝐐}}_A^{⟂\left(2\right)}\widehat{𝐑}\right\}}{M-P}{\widehat{𝐐}}_A^{⟂\left(2\right)}\right),$

end for

${\mu }_{\mathrm{opt}}=\arg \min U\left(\mu \right)$

$\mathrm{DOAs}=\{{\theta }_1,\cdots ,{\theta }_q,{\hat{\theta }}_{q+1}\left({\mu }_{\mathrm{𝑜𝑝𝑡}}\right),\cdots ,{\hat{\theta }}_P\left({\mu }_{\mathrm{𝑜𝑝𝑡}}\right)\}$

Table 2. TABLE II: Computational complexity applying the SVD

Algorithm	Multiplications
MUSIC [61]	$\frac{180}{\Delta }\left[M^2+M\left(2-P\right)-P\right]+8M{N}^2$
root-MUSIC[62]	$2M^3-M^{2}P+8M{N}^2$
AVF [78]	$\frac{180}{\Delta }\left[M^2\left(3P+1\right)+M\left(4P-2\right)+P+2\right]$
	$+M^{2}N$
CG [79]	$\frac{180}{\Delta }\left[M^2\left(P+1\right)+M\left(6P+2\right)+P+1\right]+M^{2}N$
ESPRIT[63]	$2M^{2}P+M\left(P^2-2P+8N^2\right)+8P^3-P^2$
Two-Step	$\tau [3M^3+M^2(3P+2)+M(\frac{5}{2}{P}^2-\frac{3}{2}P+8N^2)$
KAI-ESPRIT	$+P^2(\frac{17}{2}P+\frac{1}{2})+1]$
(Proposed)	$+[2M^3+M^2\left(3P\right)+M(\frac{5}{2}{P}^2-\frac{3}{2}P+8N^2)$
	$+P^2(\frac{17}{2}P+\frac{1}{2})]$

Table 3. TABLE III: Computational complexity applying the SVD

Algorithm	Additions
MUSIC [61]	$\frac{180}{\Delta }\left[M^2+M\left(1-P\right)-2\right]+8M{N}^2$
root-MUSIC[62]	$2M^3-M^{2}P+M\left(8N^2-2\right)+1$
AVF [78]	$\frac{180}{\Delta }\left[4P{M}^2+M\left(2P-3\right)-3P+2\right]$
	$+M^2\left(N-1\right)$
CG [79]	$\frac{180}{\Delta }\left[M^2\left(P+1\right)+M\left(5P+1\right)-3P-2\right]$
	$+M^2\left(N-1\right)$
ESPRIT[63]	$2M^{2}P+M\left(P^2-4P+8N^2\right)+8P^3$
	$-P^2+2P$
Two-Step	$\tau [3M^3+M^2(3P-1)$
	$+M\left(\frac{5}{2}{P}^2-\frac{9}{2}P+8N^2+2\right)$
KAI-ESPRIT	$+P(8P^2-2P-\frac{5}{2})]$
(Proposed)	$+[2M^3+M^2(3P-3)$
	$+M\left(\frac{5}{2}{P}^2-\frac{9}{2}P+8N^2+1\right)$
	$+P(8P^2-2P-\frac{5}{2})]$

Equations

$$\mathbf{x}(i)=\mathbf{A}s(i)+\mathbf{n}(i),i=1,2,\dots ,N,$$

$$\mathbf{a}({\theta }_n)=[1,e^{j2\pi \frac{d}{{\lambda }_c}\sin {\theta }_n},\dots ,e^{j2\pi (M-1)\frac{d}{{\lambda }_c}\sin {\theta }_n}{]}^T,$$

$$\mathbf{R}=\mathbb{E}\left[\mathbf{x}(i){\mathbf{x}}^H(i)\right]=\mathbf{A}{\mathbf{R}}_{ss}{\mathbf{A}}^H+{\sigma }_n^2{\mathbf{I}}_M,$$

$$\hat{\mathbf{R}}=\frac{1}{N}\sum _{i=1}^N\mathbf{x}(i){\mathbf{x}}^H(i),$$

$$\hat{\mathbf{R}}=\frac{1}{N}\sum _{i=1}^N(\mathbf{A}s(i)+\mathbf{n}(i))(\mathbf{A}s(i)+\mathbf{n}(i){)}^H$$

$$=\mathbf{A}\left\{\frac{1}{N}\sum _{i=1}^N\mathbf{s}(i){\mathbf{s}}^H(i)\right\}{\mathbf{A}}^H+\frac{1}{N}\sum _{i=1}^N\mathbf{n}(i){\mathbf{n}}^H(i)+$$

$$\mathbf{A}\left\{\frac{1}{N}\sum _{i=1}^N\mathbf{s}(i){\mathbf{n}}^H(i)\right\}+\left\{\frac{1}{N}\sum _{i=1}^N\mathbf{n}(i){\mathbf{s}}^H(i)\right\}{\mathbf{A}}^H$$

$$\hat{\mathbf{s}}(i)=\arg \underset{\begin{array}{c}\mathbf{s}\end{array}}{\min }\parallel \mathbf{x}(i)-\hat{\mathbf{A}}\mathbf{s}{\parallel }_2^2.$$

$$\hat{\mathbf{s}}(i)=({\hat{\mathbf{A}}}^H\hat{\mathbf{A}}{)}^{-1}\hat{\mathbf{A}}\mathbf{x}(i)$$

$$\hat{\mathbf{n}}(i)=\mathbf{x}(i)-\hat{\mathbf{A}}\hat{\mathbf{s}}(i).$$

$$\mathbf{V}$$

$$=\hat{\mathbf{A}}\{\frac{1}{N}\sum _{i=1}^N({\hat{\mathbf{A}}}^H\hat{\mathbf{A}}{)}^{-1}{\hat{\mathbf{A}}}^H\mathbf{x}(i)$$

$$\times ({\mathbf{x}}^H(i)-{\mathbf{x}}^H(i)\hat{\mathbf{A}}({\hat{\mathbf{A}}}^H\hat{\mathbf{A}}{)}^{-1}{\hat{\mathbf{A}}}^H)\}$$

$$={\hat{\mathbf{Q}}}_A\left\{\frac{1}{N}\sum _{i=1}^N\mathbf{x}(i){\mathbf{x}}^H(i)\left({\mathbf{I}}_M-{\hat{\mathbf{Q}}}_A\right)\right\}$$

$$={\hat{\mathbf{Q}}}_A\hat{\mathbf{R}}{\hat{\mathbf{Q}}}_A^{\perp },$$

$${\hat{\mathbf{Q}}}_A\triangleq \hat{\mathbf{A}}({\hat{\mathbf{A}}}^H\hat{\mathbf{A}}{)}^{-1}{\hat{\mathbf{A}}}^H$$

$${\hat{\mathbf{Q}}}_A^{\perp }\triangleq {\mathbf{I}}_M-{\hat{\mathbf{Q}}}_A$$

$${\hat{\mathbf{R}}}^{(n+1)}=\hat{\mathbf{R}}-\mu ({\mathbf{V}}^{(n)}+{\mathbf{V}}^{(n)H}),\text{n}=\text{1}$$

$$\begin{array}[]{ccc}\mathbf{\hat{R}}^{(n+1)}\end{array}=\left[\begin{array}[]{ccc}\mathbf{\hat{U}}_{s}&\mathbf{\hat{U}}_{n}\\ \end{array}\right]\left[\begin{array}[]{ccc}\mathbf{\hat{\Lambda}}_{s}&0\\ 0&\mathbf{\hat{\Lambda}}_{n}\end{array}\right]\left[\begin{array}[]{ccc}\mathbf{\hat{U}}_{s}^{H}\\ \mathbf{\hat{U}}_{n}^{H}\end{array}\right]$$

$${\mathbf{J}}_1{\hat{\mathbf{U}}}_s\mathbf{\Psi }\approx {\mathbf{J}}_2{\hat{\mathbf{U}}}_s$$

$$\hat{\mathbf{\Psi }}=\arg \underset{\begin{array}{c}\mathbf{\Psi }\end{array}}{\min }\parallel {\mathbf{J}}_2{\hat{\mathbf{U}}}_s-{\mathbf{J}}_1{\hat{\mathbf{U}}}_s\mathbf{\Psi }{\parallel }_F={\left({\mathbf{J}}_1{\hat{\mathbf{U}}}_s\right)}^{†}{\mathbf{J}}_2{\hat{\mathbf{U}}}_s,$$

$${\gamma }_n=\arg \left({\lambda }_n\right)$$

$${\hat{\theta }}_n=\arcsin \left(\frac{{\gamma }_n{\lambda }_c}{2\pi d}\right)$$

$$U(\mu )=\ln \det \left({\hat{\mathbf{Q}}}_A^{(2)}\hat{\mathbf{R}}{\hat{\mathbf{Q}}}_A^{(2)}+\frac{\mathrm{Trace}\{{\hat{\mathbf{Q}}}_A^{\perp (2)}\hat{\mathbf{R}}\}}{\mathrm{M}-\mathrm{P}}{\hat{\mathbf{Q}}}_A^{\perp (2)}\right)$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{MUSIC}\right)}\approx \frac{60M^3\Delta }{180M^2}\approx \frac{M}{30}$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{CG}\right)}\approx \frac{60M^3\Delta }{180M^{2}P}\approx \frac{M}{30P}$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{AVF}\right)}\approx \frac{60M^3\Delta }{180M^{2}3P}\approx \frac{M}{90P}$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{root}-\mathrm{MUSIC}\right)}\approx \frac{60M^3}{2M^3}\approx 30$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{original}\mathrm{ESPRIT}\right)}\approx \frac{60M^3}{2M^{2}P}\approx \frac{30M}{P}$$

$$\frac{\mathcal{O}\left(\mathrm{Two}-\mathrm{Step}\mathrm{KAI}\mathrm{ESPRIT}\right)}{\mathcal{O}\left(\mathrm{AVF}\right)}\approx \frac{60M^3\Delta }{180M^{2}4P}\approx \frac{M}{90P}$$

$$\text{centering}\mathrm{RMSE}=\sqrt{\frac{1}{LP}\sum _{l=1}^L\sum _{p=1}^P({\theta }_p-\hat{}{\theta }_p(l))},\text{@add@centering}$$

Full text

I Introduction

In array signal processing, direction-of-arrival (DOA) estimation is a key task in a broad range of important applications including radar and sonar systems, wireless communications and seismology [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 28, 17, 16, 19, 18, 20, 22, 21, 24, 25, 26, 31, 32, 29, 30, 23, 56, 35, 33, 34, 39, 37, 38, 39, 40, 41, 42, 43, 44, 46, 48, 49, 53, 54, 55, 56, 58, 59, 60]. Classical high-resolution methods for DOA estimation such as the multiple signal classification (MUSIC) method [61], the root-MUSIC algorithm [62], the estimation of signal parameters via rotational invariance techniques (ESPRIT) [63] and other recent subspace techniques [64, 65, 66] are based on estimating the signal and noise subspaces from the sample covariance matrix. The accuracy of the estimates of the covariance matrix is of fundamental importance in parameter estimation. In practical scenarios, only a limited number of samples is available and when the number of samples is comparable to the number of sensor array elements, the estimated and the true subspaces can significantly diverge. This problem has been dealt with using random matrix theory in [67, 68, 69], and the development of G-MUSIC, which considers the asymptotic situation when both the sample size and the number of array elements tend to infinity at the same rate. It is then deduced that the introduced method more accurately describes the case in which these two quantities are finite and similar in magnitude [69].

In this work, we take into account a different approach to improve the quality of the sample covariance matrix estimate in the finite sample size region. Inspired by the structural approach of [72, 73] and the use of prior knowledge about signals [70, 71], we develop an ESPRIT-based technique that exploits both prior knowledge about signals and the structure of the covariance matrix to improve the estimation accuracy. Our approach determines the value of a scaling factor that reduces the undesirable terms causing perturbations in the estimation of the signal and noise subspaces in an iterative manner, resulting in better estimates. This is done by choosing the set of DOA estimates that have higher likelihood of being the set of true DOAs. Furthermore, whereas in [72, 73], the undesirable terms are calculated based on steering vectors of estimates, in the proposed method this computation makes use not only of those estimates but also of the available prior knowledge of DOAs. Considering a practical scenario, this task can be achieved using signals coming from base stations or from static users in the system.

The remainder of this paper is organized as follows. Section II describes the system model and its parameters. Section III presents the proposed Two-Step KAI-ESPRIT algorithm. In section IV, we deal with the computational complexity analysis by means of counting the multiplications and additions involved in the processing. Section V illustrates and discusses the simulation results and finally, the concluding remarks are drawn in Section VI.

II System Model and Background

Let us assume that P uncorrelated narrowband signals from far-field sources impinge on a uniform linear array (ULA) of $M\ (M>\textit{P})$ sensor elements from directions ${\bm{\theta}}=[\theta_{1},\theta_{2},\ldots,\theta_{P}]^{T}$. We also consider that the sensors are spaced from each other by a distance $d\leq\frac{\lambda_{c}}{2}$, where $\lambda_{c}$ is the signal wavelength, and that without loss of generality, we have ${\frac{-\pi}{2}\leq\theta_{1}\leq\theta_{2}\ldots\leq\theta_{P}\leq\frac{-\pi}{2}}$.

The $i$th data snapshot of the $M$-dimensional array output vector can be modeled as

[TABLE]

where $\bm{s}(i)=[s_{1}(i),\ldots,s_{P}(i)]^{T}\in\mathbb{C}^{\mathit{P\times 1}}$ represents the zero-mean source data vector, $\bm{n}(i)\in\mathbb{C}^{\mathit{M\times 1}}$ is the vector of white circular complex Gaussian noise with zero mean and variance $\sigma_{n}^{2}$, and $N$ denotes the number of available snapshots. The Vandermonde matrix $\bm{A}(\bm{\Theta})=[\bm{a}(\theta_{1}),\ldots,\bm{a}(\theta_{P})]\in\mathbb{C}^{\mathit{M\times P}}$, known as the array manifold, contains the array steering vectors $\bm{a}(\theta_{j})$ corresponding to the $n$th source, which can be expressed as

[TABLE]

where $n=1,\ldots,P$. Using the fact that $\bm{s}(i)$ and $\bm{n}(i)$ are modeled as uncorrelated linearly independent variables, the $M\times M$ signal covariance matrix is calculated by

[TABLE]

where the superscript H and $\mathbb{E}[\cdot]$ in $\bm{R}_{ss}=\mathbb{E}[\bm{s}(i)\bm{s}^{H}(i)]$ and in $\mathbb{E}[\bm{n}(i)\bm{n}^{H}(i)]=\sigma_{n}^{2}\bm{I}_{M}$ denote the Hermitian transposition and the expectation operator and $\bm{I}_{M}$ stands for the $M\times M$ identity matrix. Since the true signal covariance matrix is unknown, it must be estimated and a widely-adopted approach is the sample average formula given by

[TABLE]

whose estimation accuracy is dependent on the number of snapshots.

III Proposed Two-Step KAI-ESPRIT Algorithm

In this section, we present the proposed two-step KAI-ESPRIT algorithm and detail its main features.

We can start by expanding (4) using (1) as follows:

[TABLE]

The first two terms of $\bm{\hat{R}}$ in (5) can be considered as estimates of the two summands of $\bm{R}$ given in (3), which represent the signal and the noise components, respectively. The last two terms in (5) are undesirable by-products, which can be seen as estimates for the correlation between the signal and the noise vectors. The system model under study is based on noise vectors which are zero-mean and also independent of the signal vectors. Thus, the signal and noise components are uncorrelated to each other. As a consequence, for a large enough number of samples N, the last two terms of (5) tend to zero. Nevertheless, in practice the number of available samples can be limited. In such situations, the last two terms in (5) may have significant values, which causes the deviation of the estimates of the signal and the noise subspaces from the true signal and noise ones. The key point of the proposed Two-Step KAI-ESPRIT algorithm is to modify the sample data covariance matrix in the second step based on the estimates obtained at the first step and the available known DOAs. The modified covariance matrix is computed by deriving a scaled version of the undesirable terms from $\bm{\hat{R}}$.

The steps of the proposed algorithm are listed in Table I. The algorithm starts by computing the sample data covariance matrix (4). Next, the DOAs are estimated using the ESPRIT algorithm. The superscript $(\cdot)^{(1)}$ refers to the estimation task performed in the first step. In the second step, the Vandermonde matrix is formed using the DOA estimates. Then, the amplitudes of the sources are estimated such that the square norm of the differences between the observation vector and the vector containing estimates and the available known DOAs is minimized. This problem can be formulated as

[TABLE]

The minimization of (6) is achieved using the least squares technique and the solution is described by

[TABLE]

The noise component is then estimated as the difference between the estimated signal and the observations made by the array, as given by

[TABLE]

After estimating the signal and noise vectors, the third term in (5) can be computed as

[TABLE]

where

[TABLE]

is an estimate of the projection matrix of the signal subspace, and

[TABLE]

is an estimate of the projection matrix of the noise subspace.

Lastly, the modified data covariance matrix is calculated by computing a scaled version of the estimated terms from the initial sample data covariance matrix as given

[TABLE]

The scaling or reliability factor $\mu$ increases from 0 to 1 incrementally, resulting in modified data covariance matrices. Each of them gives origin to new estimated DOAs denoted by the superscript $(\cdot)^{(2)}$ by using the ESPRIT algorithm, as briefly described ahead. Assuming the rank P is known, the eigenvalue decomposition (EVD) of the modified data covariance matrix (12) yields

[TABLE]

where the matrix $\mathbf{\hat{U}}_{s}$ $\in$ $\mathbb{C}^{\mathit{M\times P}}$ represents the signal subspace and the matrix $\mathbf{\hat{U}}_{n}$ $\in$ $\mathbb{C}^{\mathit{M\times(M-P)}}$ represents the the noise subspace respectively. The diagonal matrices $\mathbf{\hat{\Lambda}_{s}}$ and $\mathbf{\hat{\Lambda}_{n}}$ contain the P largest and the M-P smallest eigenvalues, respectively. We can form a twofold subarray configuration, as each row of the (Vandermonde) array steering matrix $\bm{A}(\bm{\Theta})$ corresponds to one particular sensor element of the antenna array. The subarrays are specified by two $\mathit{(s\times M)}$-dimensional selection matrices $\mathbf{J_{1}}$ and $\mathbf{J_{2}}$ which choose $\mathit{s}$ elements of the $\mathit{M}$ existing sensors respectively, where $\mathit{s}$ is in the range $\mathit{P\leq s<M}$. For maximum overlap, the matrix $\mathbf{J_{1}}$ selects the first $\mathit{s=M-1}$ elements and the matrix $\mathbf{J_{2}}$ selects the last $\mathit{s=M-1}$ rows of $\bm{A}(\bm{\Theta})$.

Since the matrices $\mathbf{J_{1}}$ and $\mathbf{J_{2}}$ have now been computed, we can estimate the operator $\mathbf{\Psi}$ by solving the approximation of the shift invariance equation (14) given by

[TABLE]

Using the least square (LS) method, which yields

[TABLE]

where $\parallel\cdot\parallel_{F}$ denotes the Frobenius norm and $\left(\cdot\right)^{\dagger}$ stands for the pseudo-inverse.

Lastly, the eigenvalues $\lambda_{n}$ of $\hat{\mathbf{\Psi}}$ contain the estimates of the spatial frequencies $\gamma_{n}$ computed as

[TABLE]

so that the DOAs can be calculated as

[TABLE]

where for (16) and (17) $\mathrm{n=1,\cdots,P}$.

Then, a new Vandermonde matrix is formed by the steering vectors of those newly estimated DOAs and the available known DOAs. By using the new matrix, it is possible to compute the new estimates of the projection matrices of the signal $\mathbf{\hat{Q}}_{A}^{(n)}$ and the noise $\mathbf{\hat{Q}}_{A}^{(n)\perp}$ subspaces, both for n=2.

Next, employing the new estimates of the projection matrices, the initial sample data matrix, $\bm{\hat{R}}$, and the number of sensors and sources, the stochastic maximum likelihood objective function [74] is computed for each value of $\mu$, as follows:

[TABLE]

The previous computation selects the set of unavailable DOA estimates that have a higher likelihood. Then, the set of estimated DOAs corresponding to the optimum value of $\mu$ that minimizes (18) is determined. Finally, the output of the proposed Two-Step KAI-ESPRIT algorithm is formed by the set of available known DOAs and the estimates of the unavailable DOAs, as described in Table I.

IV Computational Complexity Analysis

In this section, we evaluate the computational cost of the proposed Two-Step KAI-ESPRIT algorithm which is compared to the following subspace methods: MUSIC, Root-MUSIC, Conjugate Gradient (CG) and Auxiliary Vector Filtering (AVF). All of them made use of Singular Value Decomposition (SVD) of the sample covariance matrix (4). The computational cost in terms of number of multiplications and additions is depicted in Tables II and III, where $\mathrm{\Delta}$ is the search step and $\mathrm{\tau}=\frac{1}{increment}+1$, and this increment is defined in Table I.

As can be seen, the proposed algorithm shows a relatively high computational burden with $\mathit{O(\tau(3M^{3}))}$, where $\tau$ is typically an integer inside [1 20]. This is motivated by two reasons. The first is that the modified data covariance matrix (12) needs to be computed $\tau$ times. The second is the need for three matrix multiplications of order $\left[M\times M\right]$ that define the undesirable by-products(5), (9) to be subtracted from the sample covariance matrix (12). A brief comparison between the computational cost of the proposed algorithm and the others listed in Table II can be done considering the dominant terms in the required multiplications. For this purpose we suppose typical values of the increment $\tau=20$ and of the search step $\Delta=0.1$ degrees. For MUSIC, CG and AVF, which require peak searches, the comparisons yield

[TABLE]

[TABLE]

[TABLE]

For root-MUSIC and the original ESPRIT, which do not involve peak searches, we have

[TABLE]

[TABLE]

As can be noticed in (19),(20),(21),(22) and (23), the order of multiplications of the proposed algorithm applied to P=4 supposed signals impinging a ULA formed with M=40 sensors is approximately $1.3\times$ the order of multiplications of the MUSIC, $0.3\times$ the order of multiplications of the CG and $0.1\times$ the order of multiplications of the AVF algorithms, respectively. Therefore, in this particular scenario, the number of multiplications of Two-Step KAI ESPRIT can be considered to be approximately equal or less the number of multiplications of the algorithms that require peak search which are considered in this work. It can also be seen that the number of multiplications of the proposed algorithm is roughly $30\frac{M}{P}\times$ the number of the multiplications of the original ESPRIT and $(30\frac{M}{P}\times$ the number of multiplications of root-MUSIC, which do not require extensive peak searches. To sum up, the comparison with the algorithms that require peak search allow us to consider that the relatively high computational burden, which is associated with extra matrix multiplications and the increment $\tau$ applied to cancel the undesirable by-products, is not a too high cost to be paid for the improved performance achieved. Similar results can be obtained for the order of additions, except for the comparison with the AVF algorithm, which yields

[TABLE]

what means that in the same scenario described before the order of multiplications of our proposed algorithm is approximately $(0.08\times)$ the order of multiplications of the AVF algorithm.

V Simulations

In this section, we examine the performance of the proposed Two-Step KAI-ESPRIT algorithm in terms of probability of resolution and RMSE and compare them to the conventional ESPRIT [63], the Iterative ESPRIT (IESPRIT), which is also developed here by combining the approach in [73] that exploits knowledge of the structure of the covariance matrix and its perturbation terms and the standard ESPRIT, and the Knowledge-Aided ESPRIT (KA-ESPRIT) [70] using general linear combination [75]. We employ a ULA with * M=40* sensors, inter-element spacing $\Delta=\frac{\lambda_{c}}{2}$ and assume there are four uncorrelated complex Gaussian signals with equal power impinging on the array. The closely-spaced sources are separated by $2^{o}$, at $\mathrm{(}13^{o},15^{o},17^{o},19^{o})$, and the number of available snapshots is N=10. Furthermore, we presume a priori knowledge of the two last true DOAs at $\mathrm{(}17^{o},19^{o})$ only in the proposed Two-Step KAI-ESPRIT and in the KA-ESPRIT. In Fig. 1, we show the probability of resolution versus SNR. We take into account the criterion [76], in which two sources with DOA $\theta_{1}$ and $\theta_{2}$ are said to be resolved if their respective estimates $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are such that both $\left|\hat{\theta}_{1}-\theta_{1}\right|$ and $\left|\hat{\theta}_{2}-\theta_{2}\right|$ are less than $\left|\theta_{1}-\theta_{2}\right|/2$. The proposed Two-Step KAI-ESPRIT algorithm outperforms KA-ESPRIT [63, 70], IESPRIT and the standard ESPRIT.

In Fig. 2, it is shown the RMSE of the two supposedly unknown DOAs versus SNR. For the computations we adopted the expression

[TABLE]

where L=number of trials. Alternatively, Fig. 2 can be expressed in terms of dB as shown in Fig. 3, where the term CRB refers to the square root of the deterministic Cramér-Rao bound [77]. The results show the superior performance of the proposed Two-Step KAI-ESPRIT algorithm in terms of RMSE and also of probability of resolution. In particular, the proposed technique can obtain a higher probability of resolution and a lower RMSE than existing techniques for lower values of SNR. According to Figs. 1 and 2, Two-Step KAI-ESPRIT can obtain for the same performance in RMSE or probability of resolution gains in SNR that range from $0.5$ to $3.0$ dB.

VI Conclusions

We have proposed in this work the Two-Step KAI-ESPRIT algorithm which exploits prior knowledge of source signals and the structure of the covariance matrix and its perturbations. The proposed Two-Step KAI-ESPRIT algorithm can obtain significant gains in RMSE or probability of resolution performance over previously reported techniques, and has excellent potential for applications with short data records in large-scale antenna systems for wireless communications, radar and other large sensor arrays. The relatively high computational burden required, which is associated with extra matrix multiplications and the increment applied to reduce the undesirable by-products can be justified for the superior performance achieved. Future work will consider approaches to reducing computational cost.