Fast and Efficient Two-Dimensional DOA Estimation for Signals with Known Waveforms Using Uniform Circular Array

Abstract

This paper addresses the two-dimensional (2D) direction-of-arrival (DOA) estimation issue for signals with known waveforms but unknown amplitudes using uniform circular array (UCA). Unlike maximum likelihood (ML) methods such as decoupled maximum likelihood (DEML), parallel decomposition (PADEC) and so forth, which estimate DOA by spectrum peak search, we propose an efficient interferometer-based method with known waveforms. The proposed method first estimates spatial signature matrix based on ML method whose each column contains 2D DOA information corresponding to a source. Then, an interferometer procedure is performed to obtain 2D DOA estimate of each source from a closed-form solution separately and in parallel. Several simulation results show that the proposed method can achieve a significantly improvement in performance that coincides with the ML method as well as Cramér-Rao Bound (CRB), especially under some poor conditions, such as low SNR, fewer sensors or small snapshots. In addition, the performance will not degrade as the number of sources increases if the source signals are uncorrelated with each other. Meanwhile, it reduces a great amount of computational complexity without loss of too much accuracy. Thus, the proposed method is quite suitable for 2D DOA estimation with known waveforms in practice.

1. Introduction

Direction of arrival (DOA) estimation has been applied to a variety of civil and military fields, such as communications, sonar, air traffic control, electronic reconnaissance, etc. in the past decades [1,2]. A great many DOA estimation algorithms have been proposed, including beamforming-based methods [3], subspace-based methods [4,5,6] and sparsity-based methods [7,8,9]. These techniques generally consider the source signals to be non-cooperative signals and mainly utilize spatial properties such as the auto-correlation matrix, without taking into account any prior information. However, abundant prior information on signals of interest may be available to achieve better performance in mobile wireless communications, active radar or sonar, and many other applications; for example, known signal waveforms can be exploited to significantly improve DOA estimation accuracy, eliminate interfering signals and simplify computational complexity [10].

A series of algorithms have been developed in recent decades to deal with the problem of DOA estimation with known waveforms, including decoupled maximum likelihood (DEML) [11], coherent decoupled maximum likelihood (CDEML) [12], white coherent decoupled maximum likelihood (WCDEML) [13], parallel decomposition (PADEC) [14], subarray beamforming-based DOA (SBDOA) [15], linear operators (LP) [16] and more recent research [17,18,19,20,21,22,23]. The previous papers have proven that the Cramér-Rao Bound (CRB) and root-mean-square error (RMSE) of DOA estimation for signals with known waveforms are much lower than those with unknown waveforms [10,11].

Note that these algorithms mainly focus on using a uniform linear array (ULA) [15], sparse linear array (SLA) [18] or nonuniform linear array [19] to estimate signals’ one-dimensional (1D) DOA. However, a linear array only provides ${180}^{\circ }$ azimuthal coverage, while a planar array is a more attractive configuration known for its ability to provide ${360}^{\circ }$ azimuthal coverage and elevation angle information in addition. Moreover, when it comes to using a planar array such as uniform circular array (UCA) to estimate the above two-dimensional (2D) DOA, some algorithms specially designed for ULA, such as [6,14,15,16], may become invalid or need some complicated transformations [24].

On the other hand, the existing DOA estimators with known waveforms are often involved in eigenvalue decomposition [14], finding polynomial roots [13] or spectrum peak search [22], the computational complexity of which is too burdensome to realize in practice. In fact, interferometer-based methods are more suitable and widely used for DOA estimation with advantages of higher efficiency, broader flexibility and easier realization [25,26,27,28,29]. Up to now, there has been only one study concerning an interferometer-based DOA estimator with a known waveform [30]; however, it was only designed for one source using ULA.

In view of the above two aspects, we take advantage of the interferometer method to solve the 2D DOA estimation problem considering multiple sources with known waveforms using UCA. Extension to other planar arrays, such as uniform rectangular array, L-shape array, etc., is similar and straightforward. The main idea of our method is to exploit auto-correlation and cross-correlation matrix of known waveforms and array output vectors to estimate spatial signature matrix at first. Due to the fact that each column of spatial signature matrix is colinear with steering vector of one source, each 2D DOA is estimated through an interferometer procedure separately and in parallel. The proposed method is computationally efficient and easy to implement in practice, as all 2D DOA estimates are obtained in closed form from a series of least-squared (LS) solutions without searching the parameter space. Furthermore, we also derive the corresponding CRB of 2D DOA with known waveforms to prove the effectiveness of our method. Simulation results show that the performance of the proposed method achieves the ML method as well as the CRB under most circumstances, with the benefit of much less computational demand.

The rest of the paper is arranged as follows. Section 2 gives the model of the array output signals of UCA and known waveforms, along with some necessary assumptions. Section 3 introduces the ML and interferometer-based DOA estimator with known waveforms, respectively. Section 4 recalls CRB of 2D DOA without known waveforms and derives CRB with known waveforms. Section 5 carries out several simulations to verify the proposed method. Conclusions are summarized in Section 6.

The following notations are used in the paper: matrices and vectors are denoted by bold upper-case and lower-case letters, respectively. ${(·)}^T$, ${(·)}^{\ast }$, ${(·)}^H$ and ${(·)}^{-1}$ denote transpose, conjugate, conjugate transpose and inverse, respectively. $\mathrm{E}\{·\}$, $\mathrm{diag}\{·\}$, $\left|\right|·\left|\right|$ and ⊙ denote expectation, diagonalization, Euclidean norm and Hadamard product. $|·|$, $\mathrm{Re}\{·\}$ (or $\stackrel{⌢}{(·)}$), $\mathrm{Im}\{·\}$ (or $\stackrel{⌣}{(·)}$) and $\arg \{·\}$ denote modulus, real part, imaginary part and argument of a complex number, respectively. $I$ and $\mathbf{0}$ denote identity and zero matrix of proper size.

2. Signal Model

Consider a UCA with radius r and $M(M>K)$ isotropic sensors uniformly distributed over the circumference in the $xy$-plane, the geometry configuration of which is depicted in Figure 1. Assume that there are K co-frequency far-field narrowband signals impinging on the UCA from distinct angles $({\theta }_1,{\phi }_1),({\theta }_2,{\phi }_2),\cdots ,({\theta }_K,{\phi }_K)$, where ${\theta }_k\in [0,\pi /2)$ and ${\phi }_k\in [-\pi ,\pi )$ are elevation and azimuth angle of the k-th signal, respectively. The elevation angle is measured downward from the z-axis and the azimuth angle is measured counterclockwise from the x-axis.

The $M\times 1$ array output vector observed at the tth snapshot is modeled as:

$$x\left(t\right)=A(\theta ,\phi )s\left(t\right)+n\left(t\right),t=1,2,\cdots ,N$$

(1)

where $\theta ={[{\theta }_1,{\theta }_2,\cdots ,{\theta }_K]}^T$, $\phi ={[{\phi }_1,{\phi }_2,\cdots ,{\phi }_K]}^T$, $A(\theta ,\phi )=[a({\theta }_1,{\phi }_1),a({\theta }_2,{\phi }_2),\cdots$, $a({\theta }_K,{\phi }_K)]$ is the array manifold matrix, $a({\theta }_k,{\phi }_k)$ is the steering vector of the kth signal and ${\left\{a({\theta }_k,{\phi }_k)\right\}}_{k=1}^K$ is assumed to be linearly independent, i.e, $\mathrm{rank}\left(A\right(\theta ,\phi \left)\right)=K$, $s\left(t\right)={[s_1\left(t\right),s_2\left(t\right),\cdots ,s_K\left(t\right)]}^T$ is the source signal vector, $n\left(t\right)$ is the additive white Gaussian noise vector with zero-mean and assumed to be both spatially and temporally uncorrelated with the source signals, i.e., $\mathrm{E}\left\{n\left(t\right)n^H\left(t\right)\right\}={\sigma }_v^{2}I$, and N denotes the number of snapshots. According to geometry configuration of the UCA, $a({\theta }_k,{\phi }_k)$ is given by:

$$a({\theta }_k,{\phi }_k)={\left[a_1({\theta }_k,{\phi }_k),a_2({\theta }_k,{\phi }_k),\cdots ,a_M({\theta }_k,{\phi }_k)\right]}^T$$

(2)

where $a_m({\theta }_k,{\phi }_k)=e^{j(2\pi r/\lambda )sin{\theta }_{k}cos({\phi }_k-{\beta }_m)}$, ${\beta }_m=2\pi (m-1)/M$ and $\lambda$ denotes the wavelength.

We assume that the waveform of $s_k\left(t\right)$ is known, while its complex amplitude is unknown. Thus, $s_k\left(t\right)$ can be expressed as:

$$s_k\left(t\right)={\gamma }_k{y}_k\left(t\right),k=1,2,\cdots ,K$$

(3)

where $y_k\left(t\right)$ denotes the synchronous known waveform and ${\gamma }_k$ denotes the unknown complex amplitude due to propagation through different paths. Meanwhile, the number of signals K is known, since we even know their waveforms. Then, $x\left(t\right)$ becomes:

$$x\left(t\right)=As\left(t\right)+n\left(t\right)=A\mathsf{\Gamma }y\left(t\right)+n\left(t\right)=By\left(t\right)+n\left(t\right)$$

(4)

where $s\left(t\right)=\mathsf{\Gamma }y\left(t\right)$, $\mathsf{\Gamma }=\mathrm{diag}\left\{\gamma \right\}$ denotes the unknown amplitude matrix, $\gamma ={[{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_K]}^T$, $y\left(t\right)={[y_1\left(t\right),y_2\left(t\right),\cdots ,y_K\left(t\right)]}^T$ denotes the known waveform vector, $B=A\mathsf{\Gamma }$ denotes the spatial signature matrix and the argument $(\theta ,\phi )$ is dropped for convenience here and in the following.

Define cross-correlation matrix and auto-correlation matrix of the known waveforms and array output vectors as:

$$\left\{\begin{array}{c}{R}_{xy}\triangleq \mathrm{E}\left\{x\left(t\right)y^H\left(t\right)\right\}=\underset{N\to \infty }{lim}\frac{1}{N}\sum _{t=1}^{N}x\left(t\right)y^H\left(t\right)\hfill \\ R_{yy}\triangleq \mathrm{E}\left\{y\left(t\right)y^H\left(t\right)\right\}=\underset{N\to \infty }{lim}\frac{1}{N}\sum _{t=1}^{N}y\left(t\right)y^H\left(t\right)\hfill \end{array}\right.$$

(5)

We further assume that the K known waveforms ${\left\{y_k\left(t\right)\right\}}_{k=1}^K$ are not coherent, i.e., not completely correlated with each other, so that $R_{yy}$ is nonsingular, which is a common situation in communications. In addition, the source signals and noise are assumed to be uncorrelated, i.e., $R_{ny}=\mathbf{0}$, where $R_{ny}$ is defined similarly.

The problem of interest herein is to estimate the elevation angle $\theta$ and the azimuth angle $\phi$ from the array output vectors ${\left\{x\left(t\right)\right\}}_{t=1}^N$ and the known waveforms ${\left\{y\left(t\right)\right\}}_{t=1}^N$.

3. Direction-of-Arrival Estimators with Known Waveforms

In this section, we first derive the maximum likelihood (ML) DOA estimator for signals with known waveforms, then the interferometer-based DOA estimator with known waveforms is proposed intuitively.

3.1. Maximum Likelihood Estimator with Known Waveforms

We regard the source signals as deterministic process due to known waveforms, hence, likelihood function of ${\left\{x\left(t\right)\right\}}_{t=1}^N$ is:

$$\begin{array}{cc}\hfill P(\theta ,\phi ,\gamma ,{\sigma }_v^2)& ={\displaystyle \prod _{t=1}^N}{\displaystyle \frac{1}{{\left(\pi {\sigma }_v^2\right)}^M}}exp\left\{-{\displaystyle \frac{1}{{\sigma }_v^2}}\left|\right|x\left(t\right)-{By\left(t\right)\left|\right|}^2\right\}\hfill \\ & ={\displaystyle \frac{1}{{\left(\pi {\sigma }_v^2\right)}^{MN}}}exp\left\{-{\displaystyle \frac{1}{{\sigma }_v^2}}{\displaystyle \sum _{t=1}^N}\left|\right|x\left(t\right)-{By\left(t\right)\left|\right|}^2\right\}\hfill \end{array}$$

(6)

Since the noise in (1) is a zero-mean white Gaussian random process, it is obvious that maximizing the likelihood function with respect to $\theta$, $\phi$ and $\gamma$ is equivalent to minimizing the following cost function:

$$L(\theta ,\phi ,\gamma )=\sum _{t=1}^N\left|\right|x\left(t\right)-{By\left(t\right)\left|\right|}^2$$

(7)

In this case, the LS approach is the same as the ML method [31], so that the LS or ML estimate of $B$ is given as follows:

$$\begin{array}{cc}\hfill \widehat{B}& =\left({\displaystyle \sum _{t=1}^N}x\left(t\right)y^H\left(t\right)\right){\left({\displaystyle \sum _{t=1}^N}y\left(t\right)y^H\left(t\right)\right)}^{-1}\hfill \\ ={\widehat{R}}_{xy}{\widehat{R}}_{yy}^{-1}\hfill \end{array}$$

(8)

where ${\widehat{R}}_{xy}=1/N{\sum }_{t=1}^{N}x\left(t\right)y^H\left(t\right)$ and ${\widehat{R}}_{yy}=1/N{\sum }_{t=1}^{N}y\left(t\right)y^H\left(t\right)$ are estimates of auto-correlation and cross-correlation. As $R_{yy}$ is nonsingular, so does ${\widehat{R}}_{yy}$ when N is large enough. It is easy to find that $\widehat{B}$ is a asymptotic consistent estimate of $B$:

$$\begin{array}{cc}\hfill \underset{N\to \infty }{lim}\widehat{B}& =\underset{N\to \infty }{lim}{\widehat{R}}_{xy}{\widehat{R}}_{yy}^{-1}\hfill \\ =\underset{N\to \infty }{lim}\left[{\displaystyle \frac{1}{N}}{\displaystyle \sum _{t=1}^N}\left(By\left(t\right)+n\left(t\right)\right)y^H\left(t\right)\right]{\widehat{R}}_{yy}^{-1}\hfill \\ =\underset{N\to \infty }{lim}\left(B{\widehat{R}}_{yy}+{\widehat{R}}_{ny}\right){\widehat{R}}_{yy}^{-1}\hfill \\ =B\hfill \end{array}$$

(9)

where $\underset{N\to \infty }{lim}{\widehat{R}}_{ny}=R_{ny}=\mathbf{0}$ since the noise vectors and the known waveforms are uncorrelated.

Once we obtain $\widehat{B}$, based on above analysis and the structure of $B$, i.e., $b_k={\gamma }_{k}a({\theta }_k,{\phi }_k)$, where $b_k$ is the kth column of $B$, the problem of ML estimation of $\theta$ and $\phi$ is decoupled into K independent maximization problems as follows:

$$\left({\widehat{\theta }}_k,{\widehat{\phi }}_k\right)=\underset{\theta ,\phi }{argmax}{\left|a^H(\theta ,\phi ){\widehat{b}}_k\right|}^2,k=1,2,\cdots ,K$$

(10)

where ${\widehat{b}}_k$ is the kth column of $\widehat{B}$. Though the maximization in (10) can be performed by a 2D spectrum peak search, the computational cost is unattractive. To avoid the high complexity of 2D search, we propose a quite simple and efficient interferometer-based estimator, which reduces the search problem to a closed-form solution.

3.2. Interferometer-Based Estimator with Known Waveforms

The interferometer-based estimator with known waveforms herein comes after the cost function in (10). In other words, it is equivalent to the ML estimators in performance to some degree. Firstly, we can infer from (10) that vector ${\widehat{b}}_k$ and $a({\theta }_k,{\phi }_k)$ are colinear so that they only differ in complex amplitude ${\gamma }_k$ and share same phase-difference between different elements. Therefore, the phase-differences between any two of sensors for the kth source can be estimated by ${\widehat{b}}_k$. Finally, the estimated phase-differences are used to form a LS problem to obtain 2D DOA estimation in closed form. Note that our method is an extension of [26], as we bring in known waveforms and make use of baselines more comprehensively; therefore, our method will achieve much better performance in 2D DOA estimation.

Take the estimation of $({\theta }_k,{\phi }_k)$ as an example to demonstrate our method; the others are the same. Based on above analysis, we start with calculating the phase-difference (free of noise) between the pth and the qth sensor (baseline $(p,q)$) for the kth source as:

$$\begin{array}{cc}\hfill {\eta }_k^{p,q}& =\arg \left\{a_p({\theta }_k,{\phi }_k)a_q^{\ast }({\theta }_k,{\phi }_k)\right\}=\arg \left\{b_k\left(p\right)b_k^{\ast }\left(q\right)\right\}\hfill \\ =(2\pi r/\lambda )sin{\theta }_k\left(cos({\phi }_k-{\beta }_p)-cos({\phi }_k-{\beta }_q)\right)\hfill \\ ={\displaystyle \frac{2\pi r}{\lambda }}\left[\begin{array}{cc}cos\left({\beta }_p\right)-cos\left({\beta }_q\right)& sin\left({\beta }_p\right)-sin\left({\beta }_q\right)\end{array}\right]\left[\begin{array}{c}sin\left({\theta }_k\right)cos\left({\phi }_k\right)\\ sin\left({\theta }_k\right)sin\left({\phi }_k\right)\end{array}\right]\hfill \end{array}$$

(11)

where $b_k\left(p\right)$ denotes the pth element of $b_k$. To make sure that there is no phase ambiguity in ${\eta }_k^{p,q}$, i.e., ${\eta }_k^{p,q}\in [-\pi ,\pi )$, we assume that $r/\lambda \le 1/4$, since $|cos({\phi }_k-{\beta }_p)-cos({\phi }_k-{\beta }_q)|\le 2$. Otherwise, a great many techniques are available for solving the problem of phase ambiguity [32,33,34,35]. We refer the interested reader to these works for details and will not consider the ambiguity problem in this paper for simplicity.

Obviously, $({\theta }_k,{\phi }_k)$ can be solved by phase-difference equations constructed from (11) using the LS method. Assume that the sensor interval of every baseline $(q-p)$ is fixed at $=\mathsf{\Delta }$ or $\mathsf{\Delta }-M$. To achieve the best performance, we shall make use of baselines as much as possible. Thus, we can exploit baseline $(1,1+\mathsf{\Delta }),(2,2+\mathsf{\Delta }),\cdots ,(M-\mathsf{\Delta },M)$ and $(M-\mathsf{\Delta }+1,1),(M-\mathsf{\Delta }+2,2),\cdots ,(M,\mathsf{\Delta })$ to form the equations, define:

$${\eta }_k\left(\mathsf{\Delta }\right)\triangleq {[{\eta }_k^{1,1+\mathsf{\Delta }},{\eta }_k^{2,2+\mathsf{\Delta }},\cdots ,{\eta }_k^{M-\mathsf{\Delta },M},{\eta }_k^{M-\mathsf{\Delta }+1,1},{\eta }_k^{M-\mathsf{\Delta }+2,2},\cdots ,{\eta }_k^{M,\mathsf{\Delta }}]}^T\in {\mathbb{R}}^{M\times 1}$$

(12)

$$C\left(\mathsf{\Delta }\right)\triangleq \left[\begin{array}{cc}cos\left({\beta }_1\right)-cos\left({\beta }_{1+\mathsf{\Delta }}\right)& sin\left({\beta }_1\right)-sin\left({\beta }_{1+\mathsf{\Delta }}\right)\\ cos\left({\beta }_2\right)-cos\left({\beta }_{2+\mathsf{\Delta }}\right)& sin\left({\beta }_2\right)-sin\left({\beta }_{2+\mathsf{\Delta }}\right)\\ ⋮& ⋮\\ cos\left({\beta }_{M-\mathsf{\Delta }}\right)-cos\left({\beta }_M\right)& sin\left({\beta }_{M-\mathsf{\Delta }}\right)-sin\left({\beta }_M\right)\\ cos\left({\beta }_{M-\mathsf{\Delta }+1}\right)-cos\left({\beta }_1\right)& sin\left({\beta }_{M-\mathsf{\Delta }+1}\right)-sin\left({\beta }_1\right)\\ cos\left({\beta }_{M-\mathsf{\Delta }+2}\right)-cos\left({\beta }_2\right)& sin\left({\beta }_{M-\mathsf{\Delta }+2}\right)-sin\left({\beta }_2\right)\\ ⋮& ⋮\\ cos\left({\beta }_M\right)-cos\left({\beta }_{\mathsf{\Delta }}\right)& sin\left({\beta }_M\right)-sin\left({\beta }_{\mathsf{\Delta }}\right)\end{array}\right]\in {\mathbb{R}}^{M\times 2}$$

(13)

$$g_k\triangleq \frac{2\pi r}{\lambda }\left[\begin{array}{c}sin\left({\theta }_k\right)cos\left({\phi }_k\right)\\ sin\left({\theta }_k\right)sin\left({\phi }_k\right)\end{array}\right]$$

(14)

According to (11), we have a phase-difference equation set corresponding to the sensor interval $\mathsf{\Delta }$ as follows:

$${\eta }_k\left(\mathsf{\Delta }\right)=C\left(\mathsf{\Delta }\right)g_k$$

(15)

Consider the combination of all phase-difference equation sets corresponding to different sensor intervals $\mathsf{\Delta }=1,2,\cdots ,S$, where $S=M/2$ if the sensor number M is even, or $S=(M-1)/2$ if M is odd. The examples of baselines configuration at $M=$ 5, 6, 7 and 8 are shown in Figure 2. The rest of the cases are similar.

Define:

$${\overline{\eta }}_k\triangleq {[{\eta }_k^T\left(1\right),{\eta }_k^T\left(2\right),\cdots ,{\eta }_k^T\left(S\right)]}^T\in {\mathbb{R}}^{MS\times 1}$$

(16)

$$\overline{C}\triangleq {[C^T\left(1\right),C^T\left(2\right),\cdots ,C^T\left(S\right)]}^T\in {\mathbb{R}}^{MS\times 2}$$

(17)

Equation (15) will be further extended to:

$${\overline{\eta }}_k=\overline{C}{g}_k$$

(18)

Substitute the unambiguous phase-difference estimate ${\widehat{\eta }}_k^{p,q}=\arg \left\{{\widehat{b}}_k\left(p\right){\widehat{b}}_k^{\ast }\left(q\right)\right\}$ into (18), the LS solution of $g_k$ is given by:

$${\widehat{g}}_k={[{\widehat{g}}_k\left(1\right),{\widehat{g}}_k\left(2\right)]}^T={\left({\overline{C}}^T\overline{C}\right)}^{-1}{\overline{C}}^T{\widehat{\overline{\eta }}}_k$$

(19)

Finally, the estimates of $({\theta }_k,{\phi }_k)$ can be obtained as:

$$\left\{\begin{array}{c}{\widehat{\theta }}_k={sin}^{-1}\left(\frac{\lambda }{2\pi r}\sqrt{{\widehat{g}}_k^2\left(1\right)+{\widehat{g}}_k^2\left(2\right)}\right)\hfill \\ {\widehat{\phi }}_k=\arg \left\{{\widehat{g}}_k\left(1\right)+j{\widehat{g}}_k\left(2\right)\right\}\hfill \end{array}\right.$$

(20)

The procedure of the proposed 2D DOA estimator with known waveforms just requires a few multiplications and additions of complex numbers, where computational complexity is $MNS+O\left(4M\right)$ in the case of M sensors and N snapshots for each source. Moreover, 2D DOA estimates of K sources can be acquired separately and in parallel at the same time; hence, it is computationally efficient and easy to implement in practice.

4. Cramér-Rao Bound Derivation and Analysis

CRB is a theoretical lower bound for the covariance of any unbiased estimation. Therefore, it often serves as a benchmark to verify whether an algorithm achieves the optimal performance [36,37,38,39]. In this section, we shall give the CRB of 2D DOA estimation for signals without known waveforms from the previous works and derive the CRB with known waveforms, respectively.

According to [24,39], CRB of 2D DOA estimation for signals without known waveforms and unknown complex amplitudes, which is often referred as unconditional or stochastic CRB, is given by:

$${CRB}_u\left(\omega \right)=\frac{{\sigma }_v^2}{2N}{\left[\mathrm{Re}\left\{H\odot P_{+}^T\right\}\right]}^{-1}\in {\mathbb{C}}^{2K\times 2K}$$

(21)

where $\omega ={[{\theta }^T,{\phi }^T]}^T$, $H$ is defined as:

$$H\triangleq D^H\left[I-A{\left(A^{H}A\right)}^{-1}{A}^H\right]D\in {\mathbb{C}}^{2K\times 2K}$$

(22)

$$\left\{\begin{array}{c}D\triangleq \left[D_{\theta },D_{\phi }\right]\in {\mathbb{C}}^{M\times 2K}\hfill \\ D_{\theta }\triangleq \left[a_{{\theta }_1}^{\prime }({\theta }_1,{\phi }_1),a_{{\theta }_2}^{\prime }({\theta }_2,{\phi }_2),\cdots ,a_{{\theta }_k}^{\prime }({\theta }_K,{\phi }_K)\right]\in {\mathbb{C}}^{M\times K}\hfill \\ D_{\phi }\triangleq \left[a_{{\phi }_1}^{\prime }({\theta }_1,{\phi }_1),a_{{\phi }_2}^{\prime }({\theta }_2,{\phi }_2),\cdots ,a_{{\phi }_K}^{\prime }({\theta }_K,{\phi }_K)\right]\in {\mathbb{C}}^{M\times K}\hfill \end{array}\right.$$

(23)

$$\left\{\begin{array}{c}{a}_{\theta }^{\prime }({\theta }_k,{\phi }_k)\triangleq \frac{\partial a({\theta }_k,{\phi }_k)}{\partial {\theta }_k}={\left[\frac{\partial a_1({\theta }_k,{\phi }_k)}{\partial {\theta }_k},\frac{\partial a_2({\theta }_k,{\phi }_k)}{\partial {\theta }_k}\cdots ,\frac{\partial a_M({\theta }_k,{\phi }_k)}{\partial {\theta }_k}\right]}^T\in {\mathbb{C}}^{M\times 1}\hfill \\ a_{\phi }^{\prime }({\theta }_k,{\phi }_k)\triangleq \frac{\partial a({\theta }_k,{\phi }_k)}{\partial {\phi }_k}={\left[\frac{\partial a_1({\theta }_k,{\phi }_k)}{\partial {\phi }_k},\frac{\partial a_2({\theta }_k,{\phi }_k)}{\partial {\phi }_k}\cdots ,\frac{\partial a_M({\theta }_k,{\phi }_k)}{\partial {\phi }_k}\right]}^T\in {\mathbb{C}}^{M\times 1}\hfill \end{array}\right.$$

(24)

$P_{+}$ is defined as:

$$P_{+}\triangleq \left[\begin{array}{cc}P{A}^H{R}^{-1}AP& P{A}^H{R}^{-1}AP\\ P{A}^H{R}^{-1}AP& P{A}^H{R}^{-1}AP\end{array}\right]\in {\mathbb{C}}^{2K\times 2K}$$

(25)

where $R$ and $P$ are the array and signal covariance matrix, respectively:

$$R\triangleq \mathrm{E}\left\{x\left(t\right)x^H\left(t\right)\right\}=\underset{N\to \infty }{lim}\frac{1}{N}\sum _{t=1}^{N}x\left(t\right)x^H\left(t\right)=AP{A}^H+{\sigma }_v^{2}I\in {\mathbb{C}}^{M\times M}$$

(26)

$$P\triangleq \mathrm{E}\left\{s\left(t\right)s^H\left(t\right)\right\}=\underset{N\to \infty }{lim}\frac{1}{N}\sum _{t=1}^{N}s\left(t\right)s^H\left(t\right)\in {\mathbb{C}}^{K\times K}$$

(27)

CRB of 2D DOA estimation for signals with known waveforms and unknown complex amplitudes is a generalization of a similar result in [40] as:

$${CRB}_d\left(\omega \right)={\left[\mathsf{\Omega }-\mathrm{Re}\left\{{\left[\begin{array}{cc}\tilde{U}& \tilde{V}\end{array}\right]}^H{\tilde{W}}^{-1}\left[\begin{array}{cc}\tilde{U}& \tilde{V}\end{array}\right]\right\}\right]}^{-1}\in {\mathbb{C}}^{2K\times 2K}$$

(28)

where $\mathsf{\Omega }$ is defined as:

$$\mathsf{\Omega }\triangleq \frac{2N}{{\sigma }_v^2}\mathrm{Re}\left\{\left(D^{H}D\right)\odot P_{-}^T\right\}\in {\mathbb{C}}^{2K\times 2K}$$

(29)

$$P_{-}\triangleq \left[\begin{array}{cc}P& P\\ P& P\end{array}\right]\in {\mathbb{C}}^{2K\times 2K}$$

(30)

$\tilde{U}$, $\tilde{V}$ and $\tilde{W}$ are defined as:

$$\left\{\begin{array}{c}\tilde{U}\triangleq \frac{2}{{\sigma }_v^2}\sum _{t=1}^N{S}^H\left(t\right)A^H{D}_{\theta }S\left(t\right)=\frac{2N}{{\sigma }_v^2}\left(A^H{D}_{\theta }\right)\odot P^T\in {\mathbb{C}}^{K\times K}\hfill \\ \tilde{V}\triangleq \frac{2}{{\sigma }_v^2}\sum _{t=1}^N{S}^H\left(t\right)A^H{D}_{\phi }S\left(t\right)=\frac{2N}{{\sigma }_v^2}\left(A^H{D}_{\phi }\right)\odot P^T\in {\mathbb{C}}^{K\times K}\hfill \\ \tilde{W}\triangleq \frac{2}{{\sigma }_v^2}\sum _{t=1}^N{S}^H\left(t\right)A^{H}AS\left(t\right)=\frac{2N}{{\sigma }_v^2}\left(A^{H}A\right)\odot P^T\in {\mathbb{C}}^{K\times K}\hfill \end{array}\right.$$

(31)

$$\left\{\begin{array}{c}S\left(t\right)\triangleq \mathrm{diag}\left\{s\left(t\right)\right\}=\mathrm{diag}\{s_1\left(t\right),s_2\left(t\right),\cdots ,s_K\left(t\right)\}\in {\mathbb{C}}^{K\times K}\hfill \\ Y\left(t\right)\triangleq \mathrm{diag}\left\{y\left(t\right)\right\}=\mathrm{diag}\{y_1\left(t\right),y_2\left(t\right),\cdots ,y_K\left(t\right)\}\in {\mathbb{C}}^{K\times K}\hfill \end{array}\right.$$

(32)

The derivation details are shown in the Appendix A.

Here, we will provide some insights into ${CRB}_d\left(\omega \right)$. If the source signals are uncorrelated with each other, i.e., $P=\mathrm{diag}\{{\sigma }_{s1}^2,{\sigma }_{s2}^2,\cdots ,{\sigma }_{sK}^2\}$ is a diagonal matrix, where ${\sigma }_{sk}^2=\mathrm{E}\left\{|s_k{\left(t\right)|}^2\right\}$ is the power of $s_k\left(t\right)$, (28) for UCA will reduce to:

$$\begin{array}{cc}\hfill {CRB}_d({\theta }_k,{\phi }_k)& ={\displaystyle \frac{{\sigma }_v^2}{2N{\sigma }_{sk}^2}}{\mathrm{Re}}^{-1}\left\{\left[\begin{array}{cc}{a}_{{\theta }_k}^{\prime }({\theta }_k,{\phi }_k)a_{{\theta }_k}^{\prime }({\theta }_k,{\phi }_k)& 0\\ 0& a_{{\phi }_k}^{\prime }({\theta }_k,{\phi }_k)a_{{\phi }_k}^{\prime }({\theta }_k,{\phi }_k)\end{array}\right]\right\}\hfill \\ \hfill & ={\displaystyle \frac{{\lambda }^2{\sigma }_v^2}{4MN{\pi }^2{r}^2{\sigma }_{sk}^2}}{\mathrm{diag}}^{-1}\left\{{cos}^2\left({\theta }_k\right),{sin}^2\left({\theta }_k\right)\right\},k=1,2,\cdots ,K\hfill \end{array}$$

(33)

as ${a_{{\theta }_k}^{\prime }}^H({\theta }_k,{\phi }_k)a_{{\phi }_k}^{\prime }({\theta }_k,{\phi }_k)=0$, ${a_{{\theta }_k}^{\prime }}^H({\theta }_k,{\phi }_k)a({\theta }_k,{\phi }_k)=0$ and ${a_{{\phi }_k}^{\prime }}^H({\theta }_k,{\phi }_k)a({\theta }_k,{\phi }_k)=0$. That is to say, in such a case, the CRB of any one of multiple sources is the same as the CRB when it is the only source. Hence, the CRB of 2D DOA estimation of different sources are independent with each other and can be calculated separately. That also means that 2D DOA estimation performance will not degrade as the number of sources increases if the source signals are uncorrelated with each other.

5. Simulation Results and Analysis

In this section, several simulations are carried out to validate the performance and efficiency of the proposed method. All the simulations were performed using Matlab R2018a on the Mac OS Monterey platform with an Intel i7 9750H [email protected] GHz and 16 GB of RAM. We consider that there are three far-field narrowband signals impinging on a UCA with $r/\lambda =1/4$ from $(41.{91}^{\circ },-68.{92}^{\circ })$, $(53.{68}^{\circ },71.{97}^{\circ })$ and $(47.{56}^{\circ },135.{28}^{\circ })$. The signals and noise are both zero-mean complex white Gaussian sequences. Without loss of generality, powers of the source signals are identical and powers of known waveforms are normalized to unit. The unknown amplitude $\gamma$ is assumed to be a random complex: real parts are determined by signal-to-noise ratio (SNR), and the imaginary parts lie randomly in $[0,2\pi )$.

Four simulations were designed to investigate influence of different parameters on 2D DOA estimation performance, i.e., SNR, number of sensors M, number of snapshots N and correlation coefficient of signals. The source signals are assumed to be uncorrelated with each other in the first three simulations except the last one. In each simulation, the proposed interferometer-based estimator with known waveforms, the ML estimator with known waveforms and the multi-direction virtual array transformation algorithm (MVATA) without known waveforms [41] are presented for performance comparison along with the unconditional CRB and the CRB with known waveforms. For brevity, these terms are abbreviated as “Interferometer with KW”, “ML with KW”, “MVATA without KW”, “CRBu” and “CRBd” in the following figures and table. To evaluate the performance of the estimators, root mean squared error (RMSE) of azimuth and elevation angle estimation at each point in the following figures are measured from 2000 independent Monte Carlo runs, which is calculated as:

$$RMSE\left(\widehat{\theta }\right)=\sqrt{\mathrm{E}\left\{{({\widehat{\theta }}_1-{\theta }_1)}^2+{({\widehat{\theta }}_2-{\theta }_2)}^2+{({\widehat{\theta }}_3-{\theta }_3)}^2\right\}}$$

(34)

and $RMSE\left(\widehat{\phi }\right)$ is defined similarly. In addition, scatter plots of 2D DOA estimates of the proposed method are also depicted for each simulation.

5.1. Simulation 1: Performance of 2D DOA Estimation with Respect to SNR

In the first simulation, the performance with respect to SNR is studied. The SNR is varied form −21 dB to 9 dB with an interval of 3 dB, while the number of sensors M is fixed at 7 and the number of snapshots is fixed at 500. As shown in Figure 3 and Figure 4, the performance of the proposed method coincides with the ML estimator and the CRB over most areas of SNR, and only slightly deviates from the CRB when SNR is less than −15 dB.

In addition,

Table 1. Average time consumed per run of the three estimators.

Estimator	MVATA without KW	ML with KW	Interferometer with KW
Time consumed	18.4 ms	34.5 ms	0.59 ms

shows the average time consumed per run of three estimators, and it is obvious that the computational efficiency of the proposed method is much higher than the ML and MVATA estimator.

5.2. Simulation 2: Performance of 2D DOA Estimation with Respect to the Number of Sensors

The second simulation shows performance with respect to the number of sensors. It varies from 5 to 15, while the SNR is fixed at −6 dB and the number of snapshots is fixed at 500. As shown in Figure 5 and Figure 6, the two estimators with known waveforms have almost the same performance, which also meet with CRB. That means our method can acquire the best performance when the number of sensors ranges from 5 to 15, which covers most situations in reality. Moreover, the results in Figure 3 and Figure 5 show that the performance of the estimators with known waveforms outperforms the estimator without known waveforms, especially for lower SNR and fewer M.

5.3. Simulation 3: Performance of 2D DOA Estimation with Respect to the Number of Snapshots

The third simulation presents the performance with respect to the number of snapshots. It ranges from 100 to 1000 with an interval of 100, while SNR is fixed at −6 dB and the number of sensors is fixed at 7. As shown in Figure 7 and Figure 8, our method works very well even when the number of snapshots is 100.

5.4. Simulation 4: Performance of 2D DOA Estimation with Respect to Correlation Coefficient

The last simulation inspects the performance against correlation between signals. We assume the first two signals are correlated and the third signal is uncorrelated with the first two. Under this assumption, the signal covariance matrix $P$ is:

$$P=\left[\begin{array}{ccc}{\sigma }_s^2& \rho {\sigma }_s^2& 0\\ {\rho }^{\ast }{\sigma }_s^2& {\sigma }_s^2& 0\\ 0& 0& {\sigma }_s^2\end{array}\right]$$

(35)

where ${\sigma }_s^2$ is the power of the three signals and $\rho$ is the correlation coefficient of the first two signals, which is defined as:

$$\rho =\frac{\mathrm{E}\left\{s_1\left(t\right)s_2^{\ast }\left(t\right)\right\}}{\sqrt{\mathrm{E}\left\{|s_1{\left(t\right)|}^2\right\}\mathrm{E}\left\{|s_2{\left(t\right)|}^2\right\}}}=\underset{N\to \infty }{lim}\frac{1}{{\sigma }_s^2}\sum _{t=1}^N{s}_1\left(t\right)s_2^{\ast }\left(t\right)$$

(36)

In this simulation, the correlation coefficient of the first two signals is varied from 0 to 0.9, while SNR is fixed at −6 dB, the number of sensors is fixed at 7 and the number of snapshots is fixed at 500. As shown in Figure 9 and Figure 10, performance of the proposed method approaches to the CRB for weak and moderately correlated signals but degrades severely at medium and strong correlated signals. More specifically, the estimation accuracy of the first two correlated signals degrade as the correlation coefficient increases, while the estimation accuracy of the third uncorrelated signal keep unchanged. The reason for this phenomenon lies in the fact that the proposed method exploits the spatial signatures matrix $\widehat{B}$ to estimate the 2D DOA. The accuracy of $\widehat{B}$ degrades along with increasing correlation, which leads to noninvertiblity for the auto-correlation matrix ${\widehat{R}}_{yy}$. Therefore, our method is more suitable for source signals with weak or moderate correlation.

6. Conclusions

We have proposed a novel interferometer-based 2D DOA estimation method for signals with known waveforms and unknown amplitude in this paper. By introducing the known waveforms, the 2D DOA estimation performance of the proposed method outperforms conventional method without known waveforms. Under most circumstances, the proposed method achieves the same performance as the ML method and coincides with the corresponding CRB. More important, the proposed method is simple and computationally efficient, since the 2D DOA estimates of multiple sources are obtained through a set of closed-form solutions without loss of too much accuracy. Consequently, our method is a better choice for 2D DOA estimation with known waveforms for uncorrelated or moderately correlated signals in wireless communications, radar and sonar system. These conclusions have been verified by the simulation results.

Moreover, our method does not perform well for medium- or strong-correlated signals. Some possible further works may consider addressing the above problem as well as evaluating of theoretical lower bound on RMSE of the proposed method and extending the interferometer-based method to other planar arrays, such as uniform rectangular array or L-shape array.