Two-Step Root-MUSIC for Direction of Arrival Estimation without EVD/SVD Computation

Abstract

Most popular techniques for super-resolution direction of arrival (DOA) estimation rely on an eigen-decomposition (EVD) or a singular value decomposition (SVD) computation to determine the signal/noise subspace, which is computationally expensive for real-time applications. A two-step root multiple signal classification (TS-root-MUSIC) algorithm is proposed to avoid the complex EVD/SVD computation using a uniform linear array (ULA) based on a mild assumption that the number of signals is less than half that of sensors. The ULA is divided into two subarrays, and three noise-free cross-correlation matrices are constructed using data collected by the two subarrays. A low-complexity linear operation is derived to obtain a rough noise subspace for a first-step DOA estimate. The performance is further enhanced in the second step by using the first-step result to renew the previous estimated noise subspace with a slightly increased complexity. The new technique can provide close root mean square error (RMSE) performance to root-MUSIC with reduced computational burden, which are verified by numerical simulations.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Subspace-based algorithms for estimating the direction of arrivals (DOAs) of multiple signals impinging on array of sensors have drawn considerable interests due to their super-resolution abilities and high estimation accuracies [1–6]. Over several decades, this topic has been extensively studied and varieties of techniques including multiple signal classification (MUSIC) [7], minimum norm (MN) [8], estimation of the signal parameters via rotational invariance techniques (ESPRIT) [9], subspace fitting (SF) [10], and maximum likelihood (ML) [11] have been developed. In subspace methods, the data matrix or a statistics matrix of the data is normally decomposed into the signal and noise subspaces. Traditionally, this decomposition is usually carried out using an eigen-decomposition (EVD) or a singular value decomposition (SVD) computation, which is computationally intensive and time-consuming. In applications like high-resolution passive sonar systems where the number of sensors is much larger than that of signals [12], the use of such methods is unattractive owing to their intensive computational implementation, and therefore, there is a need for techniques demanding less computations.

Several techniques that seek to determine signal subspace-based estimates without EVD/SVD have been developed. These include the propagator method (PM) [13], the bearing estimation without eigen-decomposition (BEWE) [14], subspace methods without eigen-decomposition (SWEDE) [15], and subspace-based method without eigen-decomposition (SUMWE) [16], in which the signal or the noise subspace is extracted by a linear operator from the array output data and DOA is finally estimated by a cost function similar to MUSIC [17] or root-MUSIC [18]. These methods are found to have significant computational saving over those that explicitly compute EVD or SVD. However, the root mean square error (RMSE) and resolution performances of linear operation-based techniques degrade severely when the signal-to-noise ratio (SNR) is low [19].

Another cross-correlation-based 2D DOA estimation (CODE) algorithm was recently proposed [20] based on an L-shaped array composed of two uniform linear arrays (ULAs). In the CODE method, the noise space is obtained through a linear operation of matrices formed from the forward cross-correlation matrix. It has been shown that CODE can provide good performance with low SNRs and small numbers of snapshots [16]. However, the accuracy of CODE is still inferior to potentially more computationally methods such as MUSIC and root-MUSIC.

The purpose of this paper is to propose an alternate low-complexity subspace-based algorithm without EVD/SVD computation to reduce the computational burden while providing acceptable accuracy for DOA estimation. We present a novel two-step root-MUSIC (TS-root-MUSIC) algorithm to achieve performances closer to the standard root-MUSIC while keeping computational complexity comparable with CODE and being lower than root-MUSIC. In particular, we consider a single ULA which can be possibly regarded as part of a more complex array such as an L-shaped one. We propose to use both the forward and backward cross-correlation matrices [21, 22] and performance enhancement by a further second-step process to improve the estimation accuracy. Numerical simulations are finally provided to illustrate the degree of performance enhancement.

This paper is organized as follows. Section 2 describes the data model of DOA estimation problem using a ULA as well as the background of a standard root-MUSIC algorithm. The proposed TS-root-MUSIC is presented in Section 3, in which the computation of forward and backward cross-correlation matrices, subspace extraction without EVD/SVD, and performance enhancement using the two-step process are discussed in detail. Section 4 summarizes the implementation of the new method and analyzes its computation complexity with comparisons to other methods including root-MUSIC and CODE. Finally, Section 5 conducts numerical simulation results to show that our algorithm costs a less computation burden than root-MUSIC while providing a comparable performance which is much better than CODE.

2. Signal Model and Standard Root-MUSIC

Assume that there are $p$ narrowband signals with unknown DOAs $\mathbf{\Theta }\triangleq {\left[{\theta }_1,{\theta }_2,\dots ,{\theta }_p\right]}^T$ simultaneously impinging from far-field on a ULA composed of $M$ sensors. Without loss of generality, assume that $M$ is an even number, that is, $M=2N$; the $M\times 1$ array output vector at time $t,t=1,2,\dots ,T$ can be written as $$\begin{array}{}\text{(1)}& \mathbf{x}\left(t\right)=\mathbf{A}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right),\end{array}$$where $\mathbf{s}\left(t\right)$ is the $p\times 1$ source waveform vector, $\mathbf{n}\left(t\right)$ is the $M\times 1$ additive white Gaussian noise (AWGN) vector, and $\mathbf{A}\triangleq \left[\mathbf{a}\left({\theta }_1\right),\mathbf{a}\left({\theta }_2\right),\dots ,\mathbf{a}\left({\theta }_p\right)\right]$ is the $M\times p$ array steering matrix. Define $z\triangleq e^{-j\cdot 2\pi \cdot d/\lambda \cdot \sin \theta }$, where $j\triangleq \sqrt{-1}$, $\lambda$ is the center wavelength and $d\le \lambda /2$ is the array interval; the $p\times 1$ steering vector can be written as $$\begin{array}{}\text{(2)}& \mathbf{a}\left(\theta \right)=\mathbf{a}\left(z\right)\triangleq {\left[1,z^1,\dots ,z^N,z^{N+1},z^{N+2},\dots ,z^{2N-1}\right]}^T.\end{array}$$

The $M\times M$ array covariance matrix can be written as $$\begin{array}{}\text{(3)}& \mathbf{R}=E\left[\mathbf{x}\left(t\right){\mathbf{x}}^H\left(t\right)\right]=\mathbf{A}E\left[\mathbf{s}\left(t\right){\mathbf{s}}^H\left(t\right)\right]{\mathbf{A}}^H+E\left[\mathbf{n}\left(t\right){\mathbf{n}}^H\left(t\right)\right]={\mathbf{A}\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}^H+{\sigma }_n^2\mathbf{I},\end{array}$$where ${\mathbf{R}}_{\mathrm{s}}\triangleq E\left\{\mathbf{s}\left(t\right){\mathbf{s}}^H\left(t\right)\right\}$ is the source covariance matrix, ${\sigma }_n^2$ is the power of AWGN, and $\mathbf{I}$ is the identity matrix. For full-rank AWGN, the EVD/SVD of matrix $\mathbf{R}$ can be written in a standard way as $$\begin{array}{}\text{(4)}& \mathbf{R}={\mathbf{E}}_{\mathrm{s}}{\mathbf{\Pi }}_{\mathrm{s}}{\mathbf{E}}_{\mathrm{s}}^H+{\mathbf{E}}_{\mathrm{n}}{\mathbf{\Pi }}_{\mathrm{n}}{\mathbf{E}}_{\mathrm{n}}^H,\end{array}$$where ${\mathbf{\Pi }}_{\mathrm{s}}$ and ${\mathbf{\Pi }}_{\mathrm{n}}$ are the two diagonal matrices of dimensions $p\times p$ and $\left(M-p\right)\times \left(M-p\right)$, respectively, and ${\mathbf{E}}_{\mathrm{s}}$ and ${\mathbf{E}}_{\mathrm{n}}$ are the so-called signal and noise matrices, which are of dimensions $M\times p$ and $M\times \left(M-p\right)$ and contain the eigenvectors relating to the $p$ significant and $M-p$ smallest eigenvalues of $\mathbf{R}$, respectively. In real-world applications, the theoretical $\mathbf{R}$ is unavailable, and it is usually estimated by $T$ snapshots of array observed data as $$\begin{array}{}\text{(5)}& \widehat{\mathbf{R}}=\frac{1}{T}\sum _{t=1}^N\mathbf{x}\left(t\right){\mathbf{x}}^H\left(t\right).\end{array}$$

Therefore, the EVD of the array covariance matrix is in fact given by $$\begin{array}{}\text{(6)}& \widehat{\mathbf{R}}={\widehat{\mathbf{E}}}_{\mathrm{s}}{\widehat{\mathbf{\Pi }}}_{\mathrm{s}}{\widehat{\mathbf{E}}}_{\mathrm{s}}^H+{\widehat{\mathbf{E}}}_{\mathrm{n}}{\widehat{\mathbf{\Pi }}}_{\mathrm{n}}{\widehat{\mathbf{E}}}_{\mathrm{n}}^H.\end{array}$$

Based on the special geometry of a ULA as well as the orthogonality between span (${\mathbf{E}}_{\mathrm{s}}$) and span (${\mathbf{E}}_{\mathrm{n}}$), the standard root-MUSIC algorithm transforms the tremendous spectral search step involved in MUSIC into a simplified polynomial rooting as [18] $$\begin{array}{}\text{(7)}& f_{\mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t}-\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}\left(z\right)\triangleq {\mathbf{a}}^T\left(z^{-1}\right){\widehat{\mathbf{P}}}_{\mathrm{n}}\mathbf{a}\left(z\right),\end{array}$$where ${\widehat{\mathbf{P}}}_{\mathrm{n}}\triangleq {\widehat{\mathbf{E}}}_{\mathrm{n}}{\widehat{\mathbf{E}}}_{\mathrm{n}}^H$ is the projection matrix of the noise subspace. For (7), there are generally $p$ pairs of roots that symmetrically distribute around the unit circle. By choosing the $p$ roots ${\widehat{z}}_i,i=1,2,\dots ,p$ which are situated the most closely to the unit circle, one can estimate signal DOAs by $$\begin{array}{}\text{(8)}& {\widehat{\theta }}_i=\arcsin \left(\frac{\lambda }{2\pi d}\angle {\widehat{z}}_i\right),\hspace{1em}i=1,2,\dots ,p.\end{array}$$

Note that the computational complexity of this polynomial rooting step is substantially lower than that of the spectral search [18]. However, as EVD computation requires about $\mathcal{O}\left(M^3\right)$ flops [20, 23], the complexity of root-MUSIC is still high, especially when large numbers of sensors are used [18].

3. The Proposed Algorithm

3.1. Noise Subspace Extraction without EVD/SVD Computation

As shown in Figure 1, we divide the ULA into two subarrays along its center position without overlapping such that there are $N$ elements in each subarray (Note that for the case $M=2N+1$, the center sensor is shared by the two subarrays and there will be $N+1$ sensors in each subarray.) Thus, both the output vectors of the two subarrays are of dimensions $N\times 1$, which can be written as $$\begin{array}{c}\text{(9)}& {\mathbf{x}}_1\left(t\right)={\mathbf{A}}_1\mathbf{s}\left(t\right)+{\mathbf{n}}_1\left(t\right),{\mathbf{x}}_2\left(t\right)={\mathbf{A}}_2\mathbf{s}\left(t\right)+{\mathbf{n}}_2\left(t\right),\end{array}$$where $$\begin{array}{c}\text{(10)}& {\mathbf{A}}_1=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ z_1^1& z_2^1& \cdots & z_p^1\\ ⋮& ⋮& \ddots & ⋮\\ z_1^{N-1}& z_2^{N-1}& \cdots & z_p^{N-1}\end{array}\right],{\mathbf{A}}_2=\left[\begin{array}{cccc}{z}_1^N& z_2^N& \cdots & z_p^N\\ z_1^{N+1}& z_2^{N+1}& \cdots & z_p^{N+1}\\ ⋮& ⋮& \ddots & ⋮\\ z_1^{2N-1}& z_2^{2N-1}& \cdots & z_p^{2N-1}\end{array}\right],\end{array}$$are the steering matrices of the two subarrays, respectively. Due to the Vandermonde structure, ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$ satisfy $$\begin{array}{}\text{(11)}& {\mathbf{J}\mathbf{A}}_1^{\ast }={\mathbf{A}}_1{\mathbf{B}}^{-\left(N-1\right)},\text{(12)}& {\mathbf{J}\mathbf{A}}_2^{\ast }={\mathbf{A}}_1{\mathbf{B}}^{-\left(2N-1\right)},\end{array}$$where $\mathbf{J}$ is a $N\times N$ exchange matrix with ones on its antidiagonal and zeros elsewhere and $\mathbf{B}$ is a $p\times p$ diagonal matrix, given by $$\begin{array}{}\text{(13)}& \mathbf{B}\triangleq \left[\begin{array}{cccc}{z}_1& 0& \cdots & 0\\ 0& z_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & z_p\end{array}\right].\end{array}$$

[figure omitted; refer to PDF]

The $N\times N$ forward cross-covariance matrix of the two subarrays can be computed by $$\begin{array}{}\text{(14)}& {\mathbf{R}}_{x_1{x}_2}^f=E\left[{\mathbf{x}}_1\left(t\right){\mathbf{x}}_2^H\left(t\right)\right]={\mathbf{A}}_{1}E\left[\mathbf{s}\left(t\right){\mathbf{s}}^H\left(t\right)\right]{\mathbf{A}}_2^H+E\left[{\mathbf{n}}_1\left(t\right){\mathbf{n}}_2^H\left(t\right)\right]={\mathbf{A}}_1{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H.\end{array}$$

Define two new $N\times 1$ virtual signal vectors as follows: $$\begin{array}{c}\text{(15)}& {\mathbf{y}}_1\left(t\right)\triangleq {\mathbf{J}\mathbf{x}}_1^{\ast }\left(t\right),{\mathbf{y}}_2\left(t\right)\triangleq {\mathbf{J}\mathbf{x}}_2^{\ast }\left(t\right).\end{array}$$

Using (11) and (12), the $N\times N$ forward cross-covariance matrices between ${\mathbf{y}}_1\left(t\right)$ and ${\mathbf{x}}_2\left(t\right)$ and the $N\times N$ backward cross-covariance matrix between ${\mathbf{y}}_2\left(t\right)$ and ${\mathbf{x}}_1\left(t\right)$ can be written as $$\begin{array}{c}\text{(16)}& {\mathbf{R}}_{y_1{x}_2}^f=E\left[{\mathbf{y}}_1\left(t\right){\mathbf{x}}_2^T\left(t\right)\right]=E\left\{\left[{\mathbf{A}}_1{\mathbf{B}}^{-\left(N-1\right)}{\mathbf{s}}^{\ast }\left(t\right)+{\mathbf{J}\mathbf{n}}_1^{\ast }\left(t\right)\right]\times {\left[{\mathbf{A}}_2\mathbf{s}\left(t\right)+{\mathbf{n}}_2\left(t\right)\right]}^T\right\}={\mathbf{J}\mathbf{A}}_1^{\ast }E\left[{\mathbf{s}}^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right]{\mathbf{A}}_2^T+\mathbf{J}E\left[{\mathbf{n}}_1^{\ast }\left(t\right){\mathbf{n}}_2^T\left(t\right)\right]={\mathbf{A}}_1{\mathbf{B}}^{-\left(N-1\right)}{\mathbf{R}}_s^{\ast }{\mathbf{A}}_2^T,{\mathbf{R}}_{y_2{x}_1}^b=E\left[{\mathbf{y}}_2\left(t\right){\mathbf{x}}_1^T\left(t\right)\right]=E\left\{\left[{\mathbf{A}}_1{\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{s}}^{\ast }\left(t\right)+{\mathbf{J}\mathbf{n}}_2^{\ast }\left(t\right)\right]\times {\left[{\mathbf{A}}_1\mathbf{s}\left(t\right)+{\mathbf{n}}_1\left(t\right)\right]}^T\right\}={\mathbf{J}\mathbf{A}}_2^{\ast }E\left[{\mathbf{s}}^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right]{\mathbf{A}}_1^T+\mathbf{J}E\left[{\mathbf{n}}_1^{\ast }\left(t\right){\mathbf{n}}_1^T\left(t\right)\right]={\mathbf{A}}_1{\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_1^T.\end{array}$$

Note that all the three theoretical cross-covariance matrices ${\mathbf{R}}_{x_1{x}_2}^f$, ${\mathbf{R}}_{y_1{x}_2}^f$, and ${\mathbf{R}}_{y_2{x}_1}^b$ are noise-free in nature, which makes it possible to extract a signal or noise subspace from those three matrices with only linear operations.

Under the assumption $p<N$, we can divide ${\mathbf{A}}_1$ into two submatrices along its row direction as $$\begin{array}{}\text{(17)}& {\mathbf{A}}_1=\left[\begin{array}{c}{\mathbf{A}}_{11}\\ {\mathbf{A}}_{12}\end{array}\right],\end{array}$$where ${\mathbf{A}}_{11}$ and ${\mathbf{A}}_{12}$ are the $p\times p$ and $\left(N-p\right)\times p$ submatrices consisting of the first $p$ and last $N-p$ rows of ${\mathbf{A}}_1$, respectively. Since ${\mathbf{A}}_1$ is a Vandermonde matrix with full rank, we have $\mathrm{rank}\left[{\mathbf{A}}_{11}\right]=p$. Therefore, the rows of ${\mathbf{A}}_{12}$ can be expressed as linear combinations of those of ${\mathbf{A}}_{11}$, and we further have $$\begin{array}{}\text{(18)}& {\mathbf{A}}_{12}-{\mathbf{C}}^H{\mathbf{A}}_{11}\triangleq {\mathbf{U}}^H{\mathbf{A}}_1=\mathbf{0},\end{array}$$where $\mathbf{C}$ is an unknown matrix of dimensions $p\times \left(N-p\right)$ and $\mathbf{U}={\left[-{\mathbf{C}}^H\mathbf{I}\right]}^H$ is of dimensions $N\times \left(N-p\right)$. Combining ${\mathbf{R}}_{x_1{x}_2}^f$, ${\mathbf{R}}_{y_1{x}_2}^f$, and ${\mathbf{R}}_{y_2{x}_1}^b$, we can form two new $N\times M$ noise-free cross-correlation matrices as follows: $$\begin{array}{}\text{(19)}& \mathbf{P}=\left[\begin{array}{cc}{\mathbf{R}}_{x_1{x}_2}^f& {\mathbf{R}}_{y_1{x}_2}^f\end{array}\right]={\mathbf{A}}_1\left[\begin{array}{cc}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H& {\mathbf{B}}^{-\left(N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_2^T\end{array}\right],\text{(20)}& \mathbf{Q}=\left[\begin{array}{cc}{\mathbf{R}}_{x_1{x}_2}^f& {\mathbf{R}}_{y_2{x}_1}^b\end{array}\right]={\mathbf{A}}_1\left[\begin{array}{cc}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H& {\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_1^T\end{array}\right].\end{array}$$

Inserting (18) into (19) and (20) leads to $$\begin{array}{}\text{(21)}& \mathbf{P}=\left[\begin{array}{c}{\mathbf{A}}_{11}\\ {\mathbf{C}}^H{\mathbf{A}}_{11}\end{array}\right]\left[\begin{array}{cc}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H& {\mathbf{B}}^{-\left(N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_2^T\end{array}\right]\triangleq \left[\begin{array}{c}{\mathbf{P}}_1\\ {\mathbf{P}}_2\end{array}\right],\text{(22)}& \mathbf{Q}=\left[\begin{array}{c}{\mathbf{A}}_{11}\\ {\mathbf{C}}^H{\mathbf{A}}_{11}\end{array}\right]\left[{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H\hspace{1em}{\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_1^T\right]\triangleq \left[\begin{array}{c}{\mathbf{Q}}_1\\ {\mathbf{Q}}_2\end{array}\right],\end{array}$$where ${\mathbf{P}}_1$ and ${\mathbf{Q}}_1$ are of dimensions $p\times M$ while ${\mathbf{P}}_2$ and ${\mathbf{Q}}_2$ are of dimensions $\left(N-p\right)\times M$, respectively, given by $$\begin{array}{}\text{(23)}& {\mathbf{P}}_1\triangleq \left[\begin{array}{c}{\mathbf{A}}_{11}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H\\ {\mathbf{A}}_{11}{\mathbf{B}}^{-\left(N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_2^T\end{array}\right],\text{(24)}& {\mathbf{P}}_2\triangleq \left[\begin{array}{c}{\mathbf{C}}^H{\mathbf{A}}_{11}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H\\ {\mathbf{C}}^H{\mathbf{A}}_{11}{\mathbf{B}}^{-\left(N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_2^T\end{array}\right],\text{(25)}& {\mathbf{Q}}_1\triangleq \left[\begin{array}{c}{\mathbf{A}}_{11}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H\\ {\mathbf{A}}_{11}{\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_2^T\end{array}\right],\text{(26)}& {\mathbf{Q}}_2\triangleq \left[\begin{array}{c}{\mathbf{C}}^H{\mathbf{A}}_{11}{\mathbf{R}}_{\mathrm{s}}{\mathbf{A}}_2^H\\ {\mathbf{C}}^H{\mathbf{A}}_{11}{\mathbf{B}}^{-\left(2N-1\right)}{\mathbf{R}}_{\mathrm{s}}^{\ast }{\mathbf{A}}_1^T\end{array}\right].\end{array}$$

Define a new $p\times 2M$ matrix ${\mathbf{D}}_1\triangleq \left[\begin{array}{cc}{\mathbf{P}}_1& {\mathbf{Q}}_1\end{array}\right]$ and a new $\left(N-p\right)\times 2M$ matrix ${\mathbf{D}}_2\triangleq \left[\begin{array}{cc}{\mathbf{P}}_2& {\mathbf{Q}}_2\end{array}\right]$. Because ${\mathbf{P}}_2={\mathbf{C}}^H{\mathbf{P}}_1$ and ${\mathbf{Q}}_2={\mathbf{C}}^H{\mathbf{Q}}_1$, we have ${\mathbf{D}}_2={\mathbf{C}}^H{\mathbf{D}}_1$. Since $\mathrm{rank}\left({\mathbf{D}}_1\right)=p$, we have $$\begin{array}{}\text{(27)}& \mathbf{C}={\left({\mathbf{D}}_1{\mathbf{D}}_1^H\right)}^{-1}{\mathbf{D}}_1{\mathbf{D}}_2^H.\end{array}$$

Now, we have computed matrix $\mathbf{C}$. According to (18), the projection matrix of the noise subspace can be extracted without EVD/SVD computation by $$\begin{array}{}\text{(28)}& {\mathbf{P}}_{\mathrm{n}}\triangleq \mathbf{U}{\left({\mathbf{U}}^H\mathbf{U}\right)}^{-1}{\mathbf{U}}^H,\end{array}$$with which we can construct a root-MUSIC-like polynomial $$\begin{array}{}\text{(29)}& f_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}1}\left(z\right)={\mathbf{a}}_1^T\left(\frac{1}{z}\right){\widehat{\mathbf{P}}}_{\mathrm{n}}{\mathbf{a}}_1\left(z\right),\end{array}$$where ${\mathbf{a}}_1\left(z\right)\triangleq {\left[1,z^1,z^2,\dots ,z^{N-1}\right]}^T$. The $p$ signal DOAs can be found by finding the phases of the $p$ roots of $f_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}1}\left(z\right)$ that lie closest to and inside the unit circle.

3.2. Performance Enhancement Using the Two-Step Process

As mentioned above, the ideal cross-covariance matrix ${\mathbf{R}}_{x_1{x}_2}^f$ is noise-free in theory. However, with $T\ll \infty$ snapshots in real-world applications, it is usually estimated by $$\begin{array}{}\text{(30)}& {\widehat{\mathbf{R}}}_{x_1{x}_2}^f=\frac{1}{T}\sum _{t=1}^T\left[{\widehat{\mathbf{x}}}_1\left(t\right){\widehat{\mathbf{x}}}_2^H\left(t\right)\right]=\frac{1}{T}\sum _{t=1}^T\left[{\mathbf{A}}_1\mathbf{s}\left(t\right)+{\mathbf{n}}_1\left(t\right)\right]{\left[{\mathbf{A}}_2\mathbf{s}\left(t\right)+{\mathbf{n}}_2\left(t\right)\right]}^H={\mathbf{A}}_1\left(\frac{1}{T}\sum _{t=1}^T\mathbf{s}\left(t\right){\mathbf{s}}^H\left(t\right)\right){\mathbf{A}}_2^H+\frac{1}{T}\left(\sum _{t=1}^T{\mathbf{n}}_1\left(t\right){\mathbf{n}}_2^H\left(t\right)\right)+{\mathbf{A}}_1\left(\frac{1}{T}\sum _{t=1}^N\mathbf{s}\left(t\right){\mathbf{n}}_2^H\left(t\right)\right)+\frac{1}{T}\left(\sum _{t=1}^T{\mathbf{n}}_1\left(t\right){\mathbf{s}}^H\left(t\right)\right){\mathbf{A}}_2^H.\end{array}$$

It can be observed by comparing (30) with (14) that the last three terms of (30) are undesirable by-products that can be taken as estimates for the correlation between the signal and noise vectors [24]. This is mainly caused by a finite number of snapshots, and only for a large enough $T$, the last three terms of (30) tend to zero. According to [24], the accuracy of ${\widehat{\mathrm{R}}}_{x_1{x}_2}^f$ can be enhanced by a two-step process as follows: $$\begin{array}{}\text{(31)}& {\widehat{\mathbf{R}}}_{x_1{x}_2}^{f,\left(2\right)}={\mathbf{R}}_{x_1{x}_2}^f-\gamma \left({\mathbf{\Delta }}_1+{\mathbf{\Delta }}_2+{\mathbf{\Delta }}_3\right),\end{array}$$where $\gamma$ is a real scaling factor between zero and one, which can be predetermined by minimizing the stochastic maximum likelihood (SML) objective function [24]. ${\mathbf{\Delta }}_1$, ${\mathbf{\Delta }}_2$, and ${\mathbf{\Delta }}_3$ are least square (LS) estimates of the last three terms in (30), given by $$\begin{array}{c}\text{(32)}& {\mathbf{\Delta }}_1\triangleq \frac{1}{T}\left(\sum _{t=1}^T{\mathbf{n}}_1\left(t\right){\mathbf{n}}_2^H\left(t\right)\right)={\widehat{\mathbf{P}}}_{{\mathbf{A}}_1}^{\perp }{\widehat{\mathbf{R}}}_{x_1{x}_2}^f{\widehat{\mathbf{P}}}_{{\mathbf{A}}_2}^{\perp },{\mathbf{\Delta }}_2\triangleq {\mathbf{A}}_1\left(\frac{1}{T}\sum _{t=1}^T\mathbf{s}\left(t\right){\mathbf{n}}_2^H\left(t\right)\right)={\widehat{\mathbf{P}}}_{{\mathbf{A}}_1}{\widehat{\mathbf{R}}}_{x_1{x}_2}^f{\widehat{\mathbf{P}}}_{{\mathbf{A}}_2}^{\perp },\\ {\mathbf{\Delta }}_3\triangleq \frac{1}{T}\left(\sum _{t=1}^T{\mathbf{n}}_1\left(t\right){\mathbf{s}}^H\left(t\right)\right){\mathbf{A}}_2^H={\widehat{\mathbf{P}}}_{{\mathbf{A}}_1}^{\perp }{\widehat{\mathbf{R}}}_{x_1{x}_2}^f{\widehat{\mathbf{P}}}_{{\mathbf{A}}_2},\end{array}$$where $$\begin{array}{}\text{(33)}& {\widehat{\mathbf{P}}}_{{\mathbf{A}}_i}\triangleq {\widehat{\mathbf{A}}}_i\left({\widehat{\mathbf{\Theta }}}^{\left(1\right)}\right){\left[{\widehat{\mathbf{A}}}_i^H\left({\widehat{\mathbf{\Theta }}}^{\left(1\right)}\right){\widehat{\mathbf{A}}}_i\left({\widehat{\mathbf{\Theta }}}^{\left(1\right)}\right)\right]}^{-1}{\widehat{\mathbf{A}}}_i^H\left({\widehat{\mathbf{\Theta }}}^{\left(1\right)}\right),\text{(34)}& {\widehat{\mathbf{P}}}_{{\mathbf{A}}_i}^{\perp }\triangleq \mathbf{I}-{\widehat{\mathbf{P}}}_{{\mathbf{A}}_i},\end{array}$$are the projection matrices of the signal and noise subspaces extracted from subarray $i,i=1,2$, respectively, and ${\widehat{\mathbf{\Theta }}}^{\left(1\right)}$ is a rough DOA estimation result reduced by using root-MUSIC with the extracted noise projection matrix ${\mathbf{P}}_{\mathrm{n}}$ given by (28).

Similarly, with $T\ll \infty$ snapshots, the other two cross-covariance matrices ${\mathbf{R}}_{y_1{x}_2}^f$ and ${\mathbf{R}}_{y_2{x}_1}^b$ are estimated by $$\begin{array}{}\text{(35)}& {\widehat{\mathbf{R}}}_{y_1{x}_2}^f=\frac{1}{T}\sum _{t=1}^T\left[\mathbf{J}{\widehat{\mathbf{x}}}_1^{\ast }\left(t\right){\widehat{\mathbf{x}}}_2^T\left(t\right)\right]={\mathbf{J}\mathbf{A}}_1^{\ast }\left\{\frac{1}{T}\sum _{t=1}^T{\mathbf{s}}^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right\}{\mathbf{A}}_2^T+\frac{1}{T}\left\{\sum _{t=1}^T{\mathbf{n}}_1^{\ast }\left(t\right){\mathbf{n}}_2^T\left(t\right)\right\}+{\mathbf{J}\mathbf{A}}_1^{\ast }\left\{\frac{1}{T}\sum _{t=1}^T{\mathbf{s}}^{\ast }\left(t\right){\mathbf{n}}_2^T\left(t\right)\right\}+\frac{1}{T}\left\{\sum _{t=1}^T{\mathbf{n}}_1^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right\}{\mathbf{A}}_2^T,\text{(36)}& {\widehat{\mathbf{R}}}_{y_2{x}_1}^b=\frac{1}{T}\sum _{t=1}^T\left[\mathbf{J}{\widehat{\mathbf{x}}}_2^{\ast }\left(t\right){\widehat{\mathbf{x}}}_1^T\left(t\right)\right]={\mathbf{J}\mathbf{A}}_2^{\ast }\left\{\frac{1}{T}\sum _{t=1}^T{\mathbf{s}}^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right\}{\mathbf{A}}_1^T+\frac{1}{T}\left\{\sum _{t=1}^T{\mathbf{n}}_2^{\ast }\left(t\right){\mathbf{n}}_1^T\left(t\right)\right\}+{\mathbf{J}\mathbf{A}}_2^{\ast }\left\{\frac{1}{T}\sum _{t=1}^T{\mathbf{s}}^{\ast }\left(t\right){\mathbf{n}}_1^T\left(t\right)\right\}+\frac{1}{T}\left\{\sum _{t=1}^T{\mathbf{n}}_2^{\ast }\left(t\right){\mathbf{s}}^T\left(t\right)\right\}{\mathbf{A}}_1^T.\end{array}$$

By comparing the last three items in ${\widehat{\mathbf{R}}}_{y_1{x}_2}^f$ and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^b$ with that in ${\widehat{\mathbf{R}}}_{x_1{x}_2}^f$, it is easy to obtain the two-step enhancements for ${\widehat{\mathbf{R}}}_{y_1{x}_2}^f$ and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^b$ as follows: $$\begin{array}{}\text{(37)}& {\widehat{\mathbf{R}}}_{y_1{x}_2}^{f,\left(2\right)}={\widehat{\mathbf{R}}}_{y_1{x}_2}^f-\gamma \left({\mathbf{\Delta }}_1^{\ast }+{\mathbf{J}\mathbf{\Delta }}_2^{\ast }+{\mathbf{\Delta }}_3^{\ast }\right),\text{(38)}& {\widehat{\mathbf{R}}}_{y_2{x}_1}^{b,\left(2\right)}={\widehat{\mathbf{R}}}_{y_2{x}_1}^b-\gamma \left({\mathbf{\Delta }}_1^T+{\mathbf{J}\mathbf{\Delta }}_3^T+{\mathbf{\Delta }}_2^T\right).\end{array}$$

With the three accuracy enhanced cross-covariance matrices, we can use ${\widehat{\mathbf{R}}}_{x_1{x}_2}^{f,\left(2\right)}$, ${\widehat{\mathbf{R}}}_{y_1{x}_2}^{f,\left(2\right)}$, and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^{b,\left(2\right)}$ instead of ${\widehat{\mathbf{R}}}_{x_1{x}_2}^f$, ${\widehat{\mathbf{R}}}_{y_1{x}_2}^f$, and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^b$ to compute a new noise projection matrix ${\widehat{\mathbf{P}}}_n^{\left(2\right)}$ by using (19), (20), (21), (22), (23), (24), (25), (26), (27), and (28), and obtain another polynomial as follows: $$\begin{array}{}\text{(39)}& f_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}2}\left(z\right)={\mathbf{a}}_1^T\left(\frac{1}{z}\right){\widehat{\mathbf{P}}}_{\mathrm{n}}^{\left(2\right)}{\mathbf{a}}_1\left(z\right),\end{array}$$

A more accurate DOA ${\widehat{\mathbf{\Theta }}}^{\left(2\right)}$ can be estimated by finding the phases of the $p$ roots of $f_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}2}\left(z\right)$ that lie closest to and inside the unit circle.

4. Summary and Complexity Analysis

With $T$ snapshots of the observed data, detailed steps for implementing the proposed two-step root-MUSIC algorithm for DOA estimation without EVD/SVD computation are summarized in Table 1.

Table 1

The proposed two-step root-MUSIC algorithm.

Step 1

(i) Estimate ${\widehat{\mathbf{R}}}_{x_1{x}_2}^f$, ${\widehat{\mathbf{R}}}_{y_1{x}_2}^f$, and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^b$ by (30), (35), and (36).

(ii) Compute $\widehat{\mathbf{P}}$ and $\widehat{\mathbf{Q}}$ by (21), (22), (23), (24), (25), (26), and (27).

(iii) Obtain the rough noise projection matrix ${\widehat{\mathbf{P}}}_{\mathrm{n}}$ by (28).

(iv) Output the rough DOA estimate ${\widehat{\mathbf{\Theta }}}^{\left(1\right)}$ by rooting (29).

Step 2

(i) Insert ${\widehat{\mathbf{\Theta }}}^{\left(1\right)}$ into (33) and (34) to compute ${\mathbf{P}}_{{\mathbf{A}}_i}$ and ${\widehat{\mathbf{P}}}_{{\mathbf{A}}_i}^{\perp },\hspace{1em}i=1,2$.

(ii) Obtain ${\widehat{\mathbf{R}}}_{x_1{x}_2}^{f,\left(2\right)}$, ${\widehat{\mathbf{R}}}_{y_1{x}_2}^{f,\left(2\right)}$, and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^{b,\left(2\right)}$ by (31), (37), and (38), respectively.

(iii) Renew the rough noise projection matrix ${\widehat{\mathbf{P}}}_{\mathrm{n}}$ by ${\widehat{\mathbf{P}}}_{\mathrm{n}}^{\left(2\right)}$.

(iv) Output final fine DOA estimate ${\widehat{\mathbf{\Theta }}}^{\left(2\right)}$ by rooting (39).

Table 2 compares the primary computational complexity of diffident algorithms including root-MUSIC, CODE [20], and the proposed method, where $M=2N$ is the sensor used for root-MUSIC for a fair comparison.

Table 2

Comparison of computational complexity.

Item	Covariance matrix	Subspace extracting	Polynomial rooting
Algorithm
Root-MUSIC	$$\mathcal{O}\left({TM}^2\right)$$	$$\mathcal{O}\left(M^3\right)$$	$$\mathcal{O}{\left(2M-1\right)}^3$$
CODE	$$\mathcal{O}\left({TN}^2\right)$$	$$\mathcal{O}\left({\left(N-p\right)}^3\right)$$	$$\mathcal{O}{\left(2N-1\right)}^3$$
Proposed, step 1	$$\mathcal{O}\left({TN}^2\right)$$	$$\mathcal{O}\left({\left(N-p\right)}^3\right)$$	$$\mathcal{O}{\left(2N-1\right)}^3$$
Proposed, step 2	$$2\times \mathcal{O}\left({TN}^2\right)$$	$$2\times \mathcal{O}\left({\left(N-p\right)}^3\right)$$	$$2\times \mathcal{O}{\left(2N-1\right)}^3$$

For root-MUSIC, estimation of the $M\times M$ covariance matrix $\widehat{\mathbf{R}}$ requires $\mathcal{O}\left({TM}^2\right)$ flops while its EVD costs $\mathcal{O}\left(M^3\right)$ flops [23]. As the polynomial in (7) is of order $2M-1$, the rooting step costs $\mathcal{O}{\left(2M-1\right)}^3$ flops [16, 20]. For both CODE and the proposed method, as ${\widehat{\mathbf{R}}}_{x_1{x}_2}^f$, ${\widehat{\mathbf{R}}}_{y_1{x}_2}^f$, and ${\widehat{\mathbf{R}}}_{y_2{x}_1}^b$ are all of reduced dimensions $N\times N$, estimations of the three cross-covariance matrices cost $\mathcal{O}\left({TN}^2\right)$ flops in total. According to (27) and (28), the computation of ${\widehat{\mathbf{P}}}_{\mathrm{n}}$ needs about $32p^3+40p^2\left(N-2p\right)+\mathcal{O}\left(p^3\right)+16p^2\left(N-2p\right)+8N{\left(N-2p\right)}^2+8{\left(N-p\right)}^2\left(N-2p\right)+\mathcal{O}\left(p^3\right)\approx \mathcal{O}\left({\left(N-p\right)}^3\right)$ flops since it occurs often in application of DOA estimation that $N\gg p$ [16, 20, 25]. Finally, as the polynomial order is reduced to half in the proposed method, the rooting step for the new method costs $\mathcal{O}{\left(2N-1\right)}^3$ flops.

Based on the above analysis and notice that $M=2N$, it can be concluded from Table 2 that the proposed method has a comparable computational complexity to CODE and is much lower than root-MUSIC.

5. Simulation Results

Numerical simulations are conducted to assess the performance of the proposed method with comparisons among CODE and root-MUSIC. Throughout the simulations, 500 independent Monte-Carlo trials have been used on a ULA with $M=2N$ sensors. For the RMSE comparison, the unconditional Cramer-Rao lower bound (CRLB) given in [26] is applied for a common reference, where the RMSE for the estimates of source incident angle $\theta$ is defined as $$\begin{array}{}\text{(40)}& \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}\triangleq \sqrt{\frac{1}{500}\sum _{i=1}^{500}{\left({\widehat{\theta }}_i-\theta \right)}^2},\end{array}$$where $\theta$ is the true DOA and ${\widehat{\theta }}_i$ represents the estimated value of the DOA of the $i\mathrm{t}\mathrm{h}$ trial. For $p$ sources, the SNR is defined as $$\begin{array}{}\text{(41)}& \mathrm{S}\mathrm{N}\mathrm{R}\triangleq 10{\log }_{10}\left[\frac{P_{\mathrm{a}\mathrm{v}\mathrm{g}}}{{\sigma }_n^2}\right],\end{array}$$where $P_{\mathrm{a}\mathrm{v}\mathrm{g}}=1/p{\sum }_{l=1}^p{P}_l$ denotes the average power of all sources, and $P_l=E\left[{\mathbf{s}}_l^2\left(t\right)\right]$ is the power of the $l\mathrm{t}\mathrm{h},l\in \left[1,p\right]$ source. According to the discussions in [24], the parameter $\gamma$ in step 2 is fixed as 0.5 throughout the simulations.

In the first simulation, we examine the maximum number of signals that can be detected by the proposed method. To this end, we plot the MUSIC spectrums computed by using the noise subspace given by CODE with that by our method in Figure 2, where $p=3$ sources at ${\theta }_1=-{20}^{°}$, ${\theta }_2=-{30}^{°}$, and ${\theta }_3={40}^{°}$ are detected by a half-wavelength ULA with $2N=8$ sensors. It is seen from Figure 2 that both the two methods can resolve the three sources well. It is also verified by Figure 2 that the maximum number of signals that can be detected by the proposed method is less than half that of the number of sensors, which is the same as CODE.

[figure omitted; refer to PDF]

In the second simulation, we evaluate the RMSE performance of DOA estimation by the proposed method and compare it with those by CODE [20] and standard root-MUSIC. For fair comparisons, we exploited the two-step process for both CODE and the proposed method. In the simulation, $p=2$ sources are estimated by using a nonuniform linear array (UNLA) with 10 sensors spaced by a half-wavelength.

First, we consider $p=2$ fixed sources at ${\theta }_1={30}^{°}$ and ${\theta }_2={36}^{°}$. In Figure 3, we plot the RMSEs as functions of the SNR, where the number of snapshot is fixed as $T=100$. In Figure 4, we plot the RMSEs as functions of the number of snapshot, where we set $\mathrm{S}\mathrm{N}\mathrm{R}=5\mathrm{d}\mathrm{B}$.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Next, to see more clearly the performance of the new method, we consider $p=2$ varying sources at ${\theta }_1={30}^{°}$ and ${\theta }_2={30}^{°}+\mathrm{\Delta }\theta$, where $\mathrm{\Delta }\theta$ varies from $1^{°}$ to ${10}^{°}$. We plot the RMSEs as functions of angle separation $\mathrm{\Delta }\theta$ in Figure 5, where we set $T=100$ and $\mathrm{S}\mathrm{N}\mathrm{R}=5\mathrm{d}\mathrm{B}$.

[figure omitted; refer to PDF]

It can be seen clearly from the three figures that the RMSE curves of the three methods including CODE, root-MUSIC, and the proposed TS-root-MUSIC decrease dramatically as the SNRs, the number of snapshots, and the values of angle separation increase. However, the RMSE curves of our method locate below that of CODE all along in all the three figures, which indicates that our method is superior to CODE. In addition, the proposed method in step 2 also outperforms CODE with the two-step process. This is because TS-root-MUSIC uses both the forward and backward cross-covariance matrices to exploit more information for DOA estimate. It is also seen from the three figures that both CODE and the proposed TS-root-MUSIC algorithm in step 2 shows a much better accuracy than those in step 1, which verifies the performance enhancement by the two-step process.

In the last simulation, we verify the efficiency of the proposed approach by comparing the simulation times of DOA estimates by CODE, root-MUSIC, and TS-root-MUSIC as functions of the number of sensors in Table 3. The simulated results are given by a PC with Intel® Core™ Duo T5870 2.0 GHz CPU and 1 GB RAM by running the Matlab codes in the same environment. It can be seen from Table 3 that the proposed method in both step 1 and step 2 costs lower simulation time than the standard root-MUSIC, especially for large numbers of sensors. Since TS-root-MUSIC in step 1 has a similar complexity as CODE, the new method costs a similar CPU time to CODE in step 1 while it costs approximately twice that in step 2. Considering that TS-root-MUSIC has a much better accuracy than CODE, our method trades off MSEs by complexity efficiently as compared to CODE.

Table 3

Comparison of simulation time in seconds.

Sensor numbers	$M=10$	$M=14$	$M=18$	$M=22$	$M=26$
Algorithm
Root-MUSIC	0.1143	0.1470	0.2165	0.3095	0.4419
CODE	0.0461	0.0470	0.0773	0.0804	0.0849
Proposed, step 1	0.0462	0.0473	0.0766	0.0810	0.0854
Proposed, step 2	0.0907	0.0933	0.1550	0.1646	0.1744

6. Conclusions

We have proposed a new computationally efficient algorithm for DOA estimation using a ULA, in which the ULA is divided into two subarrays and both the forward and backward cross-covariance matrices are exploited to obtain a signal subspace without EVD/SVD computation. A two-step process is also applied for further performance enhancement with a slightly increased complexity. Numerical simulations demonstrate that the new technique has a much lower computational complexity than root-MUSIC and CODE have while the RMSE performance of the new method is comparable to root-MUSIC and is much better than CODE.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (61501142) and the Science and Technology Program of WeiHai and project supported by the Discipline Construction Guiding Foundation in Harbin Institute of Technology (Weihai) (WH20160107).

References

[1] H. Krim, M. Viberg, "Two decades of array signal processing research: the parametric approach," IEEE Signal Processing Magazine, vol. 13 no. 4, pp. 67-94, DOI: 10.1109/79.526899, 1996.

[2] F. G. Yan, M. Jin, S. Liu, X. L. Qiao, "Real-valued MUSIC for efficient direction estimation with arbitrary array geometries," IEEE Transactions on Signal Processing, vol. 62 no. 6, pp. 1548-1560, DOI: 10.1109/TSP.2014.2298384, 2014.

[3] F. G. Yan, B. Cao, J. J. Rong, Y. Shen, M. Jin, "Spatial aliasing for efficient direction-of-arrival estimation based on steering vector reconstruction," EURASIP Journal on Advances in Signal Processing, vol. 2016,DOI: 10.1186/s13634-016-0419-y, 2016.

[4] F. Yan, Y. Shen, M. Jin, X. Qiao, "Computationally efficient direction finding using polynomial rooting with reduced-order and real-valued computations," Journal of Systems Engineering and Electronics, vol. 27 no. 4, pp. 739-745, DOI: 10.21629/jsee.2016.04.01, 2016.

[5] F.-G. Yan, Z.-K. Chen, M.-J. Sun, Y. Shen, M. Jin, "Two-dimensional direction-of-arrivals estimation based on one-dimensional search using rank deficiency principle," International Journal of Antennas and Propagation, vol. 2015,DOI: 10.1155/2015/127621, 2015.

[6] F. Yan, M. Jin, X. Qiao, "Source localization based on symmetrical MUSIC and its statistical performance analysis," Science China Information Sciences, vol. 56 no. 6,DOI: 10.1007/s11432-013-4841-6, 2013.

[7] R. Schmidt, "Multiple emitter location and signal parameter estimation," IEEE Transactions on Antennas and Propagation, vol. 34 no. 3, pp. 276-280, DOI: 10.1109/TAP.1986.1143830, 1986.

[8] Y. Ephraim, N. Merhav, H. L. Van Trees, "Min-norm interpretations and consistency of MUSIC, MODE and ML," IEEE Transactions on Signal Processing, vol. 43 no. 12, pp. 2937-2942, DOI: 10.1109/78.476437, 1995.

[9] R. Roy, T. Kailath, "ESPRIT-estimation of signal parameters via rotational invariance techniques," IEEE Transactions on Signal Processing, vol. 37 no. 7, pp. 984-995, DOI: 10.1109/29.32276, 1989.

[10] P. Pelin, "A fast minimization technique for subspace fitting with arbitrary array manifolds," IEEE Transactions on Signal Processing, vol. 49 no. 12, pp. 2935-2939, DOI: 10.1109/78.969502, 2001.

[11] J. W. Shin, Y. J. Lee, H. N. Kim, "Reduced-complexity maximum likelihood direction-of-arrival estimation based on spatial aliasing," IEEE Transactions on Signal Processing, vol. 62 no. 24, pp. 6568-6581, DOI: 10.1109/TSP.2014.2367454, 2014.

[12] M. Khosravi, M. Fakharzadeh, M. H. Bastani, "Large array null steering using compressed sensing," IEEE Signal Processing Letters, vol. 23 no. 8, pp. 1032-1036, DOI: 10.1109/LSP.2016.2580587, 2016.

[13] S. Marcos, A. Marsal, M. Benidir, "The propagator method for source bearing estimation," Signal Processing, vol. 42 no. 2, pp. 121-138, DOI: 10.1016/0165-1684(94)00122-G, 1995.

[14] P. Stoica, T. Söderström, "Statistical analysis of a subspace method for bearing estimation without eigendecomposition," IEE Proceedings F Radar and Signal Processing, vol. 139 no. 4, pp. 301-305, DOI: 10.1049/ip-f-2.1992.0042, 1992.

[15] A. Eriksson, P. Stoica, T. Söderström, "On-line subspace algorithms for tracking moving sources," IEEE Transactions on Signal Processing, vol. 42 no. 9, pp. 2319-2330, DOI: 10.1109/78.317854, 1994.

[16] J. Xin, A. Sano, "Computationally efficient subspace-based method for direction-of-arrival estimation without eigendecomposition," IEEE Transactions on Signal Processing, vol. 52 no. 4, pp. 876-893, DOI: 10.1109/TSP.2004.823469, 2004.

[17] F.-G. Yan, X.-W. Yan, J. Shi, J. Wang, S. Liu, M. Jin, Y. Shen, "MUSIC-like direction of arrival estimation based on virtual array transformation," Signal Processing, vol. 139, pp. 156-164, DOI: 10.1016/j.sigpro.2017.04.017, 2017.

[18] F.-G. Yan, Y. Shen, M. Jin, "Fast DOA estimation based on a split subspace decomposition on the array covariance matrix," Signal Processing, vol. 115,DOI: 10.1016/j.sigpro.2015.03.008, 2015.

[19] M. A. Hasan, M. R. Azimi-Sadjadi, A. A. Hasan, "Rational invariant subspace approximations with applications," IEEE Transactions on Signal Processing, vol. 48 no. 11, pp. 3032-3041, DOI: 10.1109/78.875461, 2000.

[20] G. Wang, J. Xin, N. Zheng, A. Sano, "Computationally efficient subspace-based method for two-dimensional direction estimation with L-shaped Array," IEEE Transactions on Signal Processing, vol. 59 no. 7, pp. 3197-3212, DOI: 10.1109/TSP.2011.2144591, 2011.

[21] Y. H. Choi, "ESPRIT-based coherent source localization with forward and backward vectors," IEEE Transactions on Signal Processing, vol. 58 no. 12, pp. 6416-6420, DOI: 10.1109/TSP.2010.2077634, 2010.

[22] Y.-H. Chen, C.-T. Chiang, "Kalman-based spatial domain forward-backward linear predictor for DOA estimation," IEEE Transactions on Aerospace and Electronic Systems, vol. 31 no. 1, pp. 474-479, DOI: 10.1109/7.366330, 1995.

[23] G. H. Golub, C. H. Van Loan, Matrix Computations, 1996.

[24] M. Shaghaghi, S. A. Vorobyov, "Subspace leakage analysis and improved DOA estimation with small sample size," IEEE Transactions on Signal Processing, vol. 63 no. 12, pp. 3251-3265, DOI: 10.1109/TSP.2015.2422675, 2015.

[25] N. Xi, L. Liping, "A computationally efficient subspace algorithm for 2-D DOA estimation with L-shaped array," IEEE Signal Processing Letters, vol. 21 no. 8, pp. 971-974, DOI: 10.1109/lsp.2014.2321791, 2014.

[26] P. Stoica, A. Nehorai, "Performance study of conditional and unconditional direction-of-arrival estimation," IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 38 no. 10, pp. 1783-1795, DOI: 10.1109/29.60109, 1990.