A Robust Algorithm for DOA Estimation of Coherent Sources with UCA

Ye Tian^1,2, Yonghui Huang¹, Xiaoxu Zhang^1,2, and Meiyan Lin^1,2

¹National Space Science Center, Chinese Academy of Sciences, Beijing, 100190, China
tianye171@mails.ucas.ac.cn, yonghui@nssc.ac.cn, zhangxiaoxu19@mails.ucas.ac.cn, linmeiyan18@mails.ucas.ac.cn

²University of Chinese Academy of Sciences, Beijing, 100049, China

Submitted On: February 14, 2022; Accepted On: August 10, 2022

Abstract

Direction of arrival (DOA) estimation of coherent sources with a uniform circular array (UCA) is an intractable problem. The method-of-direction-estimation (MODE) algorithm has strong superiority in handling coherent sources compared with the classical MUSIC, and ESPRIT algorithms. However, MODE is sensitive to source numbers and does not work well in the UCA scenario. In order to improve the performance of MODE, a robust DOA estimation method named UCA-PUMA (principal-eigenvector-utilization-for-modal-analysis) is proposed. The complicated non-Vandermonde structured steering vector of UCA is transformed into a virtual Vandermonde structured steering vector in mode space. The proposed method gives a closed-form solution compared with the original UCA-MODE algorithm. The performance of the UCA-PUMA method is evaluated by simulations. Simulation results demonstrate that the UCA-PUMA is more robust to source numbers than the UCA-MODE, and coherent sources can be handled without spatial smoothing. In addition, the UCA-PUMA fully takes advantage of the UCA, which is able to discriminate sources coming from a ${360}^{\circ }$ azimuthal field of view.

Index Terms: PUMA, coherent, DOA estimation, UCA, MODE.

I. INTRODUCTION

Direction-of-Arrival (DOA) estimation is attracting considerable critical attention from the array signal processing community. Several classical algorithms have been developed such as MUSIC [1], ESPRIT [2, 3], and Maximum Likelihood (ML) estimation [4, 5] in the past decades. Meanwhile, various array geometries have been considered, such as the uniform linear array (ULA), uniform circular array (UCA), and uniform rectangular array (URA), etc. Among many array configurations, the UCA has attracted much attention due to its advantages such as offering a ${360}^{\circ }$ azimuthal field of view, being easy to set up, etc. However, the widespread multipath effect of electromagnetic waves leads to coherent sources impinging on the array [6]. Some efforts have been tried in [7, 8] to deal with the problem of DOA estimation of coherent sources. Unfortunately, this problem becomes even trickier in the context of UCA [9], because of the non-Vandermonde structure of its steeringvector.

Thanks to the mode space transformation method [10], the UCA can be treated as a virtual ULA and the azimuthal isotropy of UCA are retained. The mode space transformation, also known as beam space transformation, is widely used to solve DOA estimation problems in UCA. In [11, 12], the mode space transformation method has been adopted to improve the accuracy of DOA estimation in the UCA. Also in [13], the authors have proposed a low complexity sparse beamspace DOA estimation method for UCA. The above two works focus on the one-dimensional DOA estimation of uncorrelated sources in the UCA, and the mode space transformation technique is used to improve the estimationperformance.

Furthermore, the mode space transformation technique can also play a great role in the context of DOA estimation of coherent sources in UCA. Since the covariance matrix is rank deficient in the scenarios of coherent sources, spatial smoothing techniques have been developed to address this problem [14]. Notably, the spatial smoothing technique requires the covariance matrix to be equipped with a Toeplitz structure, which is achieved by the ULA [15] or the URA [16]. Obviously, the mode space transformation technique provides the opportunity to convert the UCA into a virtual ULA, hence a similar method is also applied in the UCA scenario [17, 18]. However, in the process of implementing the spatial smoothing technique, the whole ULA or virtual ULA is divided into multiple sub-arrays, which reduces the effective array aperture and leads to severe performance degradation [19].

Particularly, the authors in [20] have proposed a ML-based approach to deal with the DOA estimation of coherent sources. Furthermore, with sophisticated mathematical proofs, this paper has shown that the proposed method does not require spatial smoothing in dealing with coherent sources. The MODE [21] is an ML type algorithm which is only required moderate computation. The MODE takes advantage of the autoregressive moving average model (ARMA) of the snapshot vector when the array steering vector has a Vandermonde structure realized by the ULA. In [22], the MODE algorithm has been extended to the UCA scenario by making full use of the mode space transformation. Nevertheless, the UCA-MODE algorithm is not robust to source number. The UCA-MODE performs well in the even source number scenario, but its performance is severely degraded in the odd case. Moreover, it performs even worse in the regime where coherent sources are hybrid with uncorrelated ones.

In this letter, we focus on the robust realization of the MODE algorithm for DOA estimation in UCA based on the recently proposed principal-eigenvector-utilization-for-modal-analysis (PUMA) method [23] combined with the mode space transformation of UCA. In this scheme, the problem of sensitivity to source number is overcome, and coherent sources can be directly handled without spatial averaging. The proposed method is named UCA-PUMA. It is worth noting that the UCA-PUMA can also deal with the scenario of complicated hybrid sources. Computer simulation is performed, and the results demonstrate that UCA-PUMA has a strong superiority over UCA-MODE [22], UCA-Smoothing [18], UCA-ESPRIT [10, 24] and UCA-Root-MUSIC [25]. The main contributions of this work are summarized as follows:

1. A robust DOA estimation algorithm, named UCA-PUMA, is proposed for coherent or hybrid sources impinging the UCA. Meanwhile, the UCA-PUMA does not require standard spatial smoothing steps that lead to loss of effective aperture.

2. The proposed UCA-PUMA algorithm takes advantage of the mode space transformation technique of UCA, and DOAs are obtained by finding the roots of a polynomial.

3. Computer simulations in various scenarios are carried out to demonstrate the superior performance of the proposed UCA-PUMA.

The remainder of this paper is organized as follows. Section 2 describes the UCA signal model and mode space transformation. In Section 3, the proposed UCA-PUMA algorithm is introduced. The simulation results and related discussions are included in Section 4. Finally, Section 5 concludes the paper.

Notations: In this letter, superscripts ${(\cdot )}^{-1}$, ${(\cdot )}^{\ast }$, ${(\cdot )}^T$, and ${(\cdot )}^H$ denote the inverse operation, complex conjugate, transpose, and conjugate transpose, respectively. $\mathrm{diag}\{\cdot \}$ and $\mathrm{tr}\{\cdot \}$ are diagonal matrix and trace operators, respectively. Boldface lowercase letters such as $\mathbf{a}$, $\mathbf{b}$ denote vectors, and boldface uppercase letters such as $\mathbf{A}$, $\mathbf{B}$ denote matrices. ${\mathbf{I}}_N$ is the $N\times N$ identity matrix. $\mathbf{\otimes }$ is the Kronecker product operator of matrices. $\mathrm{\angle }z$ means taking the argument of the complex number $z$.

II. SYSTEM MODEL

As shown in Fig. 1 and 2, consider $D$ far-field narrowband sources from directions $\mathit{\phi }={[{\phi }_1,{\phi }_2,\mathrm{\cdots },{\phi }_D]}^T$ impinging on a UCA, which consists of $N$ identical antenna elements uniformly distributed over a circle with radius $R$. The angle coordinate of the $n$-th antenna element is given by

$${\alpha }_n=\frac{2\pi (\mathit{\text{n}}-1)}{\mathit{\text{N}}}.$$

(1)

Figure 1: System model of UCA.

Figure 2: Block diagram of the DOA estimation system.

The manifold matrix $\mathbf{A}$ of the UCA is expressed as

$$\mathbf{A}=[\mathbf{a}({\phi }_1),\mathbf{a}({\phi }_2),\mathbf{\cdots },\mathbf{a}({\phi }_D)]\in {ℂ}^{N\times D},$$

(2)

where

$$\mathbf{a}({\phi }_d)=\left[\begin{array}{c}\hfill \exp \left[j2\pi \widetilde{\mathit{\text{R}}}\cos ({\phi }_d-{\alpha }_1)\right]\hfill \\ \hfill \exp \left[j2\pi \widetilde{\mathit{\text{R}}}\cos ({\phi }_d-{\alpha }_2)\right]\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill \exp \left[j2\pi \widetilde{\mathit{\text{R}}}\cos ({\phi }_d-{\alpha }_N)\right]\hfill \end{array}\right],$$

(3)

is the $d$-th steering vector in $\mathbf{A}$, and $\tilde{R}=R/\lambda$ is the radius normalized by wavelength. The $k$-th snapshot of the received signal is expressed as

$$\mathbf{x}(k)=\mathbf{As}(k)+\mathbf{n}(k),k=1,2,\mathrm{\cdots },K,$$

(4)

where $\mathbf{s}(k)$ is the source signal vector, and

$$\mathbf{n}(k)\sim 𝒞𝒩(0,{\sigma }_n^2{\mathbf{I}}_N),$$

(5)

is the additive white Gaussian noise vector, which is independent of the source signals.

When electromagnetic waves have multipath in space, coherent sources will appear at the receiving system. Without loss of generality, assume that there are $D_p$ uncorrelated sources, and there are $D_c$ sources that are coherent with previous $D_p$ sources. The total number of sources is $D=D_p+D_c$. Then the covariance matrix of the $D$ incoming sources is given by

$${\mathbf{R}}_{\mathbf{ss}}=𝔼\left[{\mathbf{ss}}^H\right]\in {ℂ}^{D\times D}.$$

(6)

The above assumption implies that ${\mathbf{R}}_{\mathbf{ss}}$ is rank deficient, i.e. its rank is $D_p<D$. Moreover, the covariance matrix of the array signal which is given by

$${\mathbf{R}}_{\mathbf{xx}}=𝔼\left[{\mathbf{xx}}^H\right]={\mathbf{AR}}_{\mathbf{ss}}{\mathbf{A}}^H+{\sigma }_n^2{\mathbf{I}}_N\in {ℂ}^{N\times N}.$$

(7)

By means of Cholesky factorization, the covariance matrix of the incoming sources is written as

$${\mathbf{R}}_{\mathbf{ss}}={\mathbf{TT}}^H,$$

(8)

with $\mathbf{P}\in {ℂ}^{D_p\times D}$. Then the following relation can be immediately obtained from 8, which is

$$\mathbf{T}=\left[\begin{array}{c}\hfill {\mathbf{I}}_{D_p}\hfill \\ \hfill \mathbf{C}\hfill \end{array}\right]\mathbf{L}=\mathbf{PL},$$

(9)

where ${\mathbf{LL}}^H=\mathbf{Q}\in {ℂ}^{D_p\times D_p}$ and $\mathbf{C}\in {ℂ}^{D_c\times D}$. Based on the above notations, ${\mathbf{R}}_{\mathbf{xx}}$ can be rewritten as

${\mathbf{R}}_{\mathbf{xx}}$	$={\mathbf{ATT}}^H{\mathbf{A}}^H+{\sigma }_n^2{\mathbf{I}}_N,$	(10)
$=(\mathbf{AP})\mathbf{Q}({\mathbf{P}}^H{\mathbf{A}}^H)+{\sigma }_n^2{\mathbf{I}}_N,$
$={\mathbf{A}}_P{\mathbf{QA}}_P^H+{\sigma }_n^2{\mathbf{I}}_N.$

By taking the eigenvalue decomposition of ${\mathbf{R}}_{\mathbf{xx}}$, we have

$${\mathbf{R}}_{\mathbf{xx}}={\mathbf{E}}_{D_p}{\mathbf{\Lambda }}_{D_p}{\mathbf{E}}_{D_p}^H+{\mathbf{E}}_n{\mathbf{\Lambda }}_n{\mathbf{E}}_n^H,$$

(11)

where ${\mathbf{E}}_{D_p}=[{\mathbf{e}}_1,\mathrm{\cdots },{\mathbf{e}}_{D_p}]$ denotes the signal subspace formed by $D_p$ eigenvectors corresponding to the $D_p$ principal largest eigenvalues.

The maximum likelihood method [4, 5] maximizes the following loss function

$$f(\phi )=\mathrm{tr}\{{\mathbf{A}}_P{[{\mathbf{A}}_P^H{\mathbf{A}}_P]}^{-1}{\mathbf{A}}_P^H\mathbf{W}\},$$

(12)

where

$\mathbf{W}$	$={\mathbf{E}}_{D_p}{\mathbf{\Gamma }}_{D_p}{\mathbf{E}}_{D_p}^H,$	(13a)
${\mathbf{\Gamma }}_{D_p}$	$=\mathrm{diag}\{{\gamma }_1,\mathrm{\cdots },{\gamma }_{D_p}\}$	(13b)
${\gamma }_d$	$={\displaystyle \frac{{({\lambda }_d-{\sigma }_n^2)}^2}{{\lambda }_d}}.$	(13c)

The above cost function (12) can be casted as

	$f(\phi )$	$=\mathrm{tr}\{\mathbf{AP}{[{\mathbf{P}}^H{\mathbf{A}}^H\mathbf{AP}]}^{-1}{\mathbf{P}}^H{\mathbf{A}}^H\mathbf{W}\},$		(14)
		$=\mathrm{tr}\{{[{\mathbf{P}}^H({\mathbf{A}}^H\mathbf{A})\mathbf{P}]}^{-1}{\mathbf{P}}^H({\mathbf{A}}^H\mathbf{WA})\mathbf{P}\},$

According to the extended Rayleigh quotient theorem, we obtain the following inequality

$$f(\phi )\le \sum _{p=1}^{D_p}{\xi }_p,$$

(15)

and when

$$\mathbf{PS}=\mathbf{K}=[{\mathbf{k}}_1,\mathrm{\cdots },{\mathbf{k}}_{D_p}],$$

(16)

is satisfied, the equal sign is obtained. Notably, $\mathbf{S}$ is an arbitrary nonsingular matrix, moreover, ${\{{\xi }_p\}}_{p=1}^{D_p}$ and ${\{{\mathbf{k}}_p\}}_{p=1}^{D_p}$ denote the $D_p$ principal eigenvalues and their associated eigenvectors of the matrix ${({\mathbf{A}}^H\mathbf{A})}^{-1}({\mathbf{A}}^H\mathbf{WA})$.

Since ${\mathbf{E}}_{D_p}$ and $\mathbf{A}$ span the same subspace, which implies ${\mathbf{E}}_{D_p}=\mathbf{AG}$ with the help of some matrix $\mathbf{G}\in {ℂ}^{N\times D_p}$. Note that

$${\mathbf{E}}_{D_p}^H{\mathbf{E}}_{D_p}=({\mathbf{E}}_{D_p}^H\mathbf{A})\mathbf{G}={\mathbf{I}}_{D_p},$$

(17)

implies that

$$\mathrm{rank}[{\mathbf{E}}_{D_p}^H\mathbf{A}]=D_p.$$

(18)

Notably, observe that

$$\mathrm{rank}[{({\mathbf{A}}^H\mathbf{A})}^{-1}({\mathbf{A}}^H\mathbf{WA})]=\mathrm{rank}({\mathbf{E}}_{D_p}^H\mathbf{A})=D_p.$$

(19)

In other words, the $D_p$ principal eigenvalues ${\{{\xi }_p\}}_{p=1}^{D_p}$ is also all eigenvalues of ${({\mathbf{A}}^H\mathbf{A})}^{-1}({\mathbf{A}}^H\mathbf{WA})$. Therefore we obtain that

$$\sum _{p=1}^{D_p}{\xi }_p=\mathrm{tr}\{{({\mathbf{A}}^H\mathbf{A})}^{-1}({\mathbf{A}}^H\mathbf{WA})\}.$$

(20)

To sum up, maximizing the cost function (12) is equivalent to minimizing the following function:

$$g(\phi )=\mathrm{tr}\{[{\mathbf{I}}_N-\mathbf{A}{({\mathbf{A}}^H\mathbf{A})}^{-1}({\mathbf{A}}^H)]\mathbf{W}\}.$$

(21)

A. UCA mode space transformation

DOA estimation in UCA based on the mode space transformation method is first proposed in [10]. Let $M$ denote the highest order mode that can be excited on a circle of radius $R$ at a reasonable strength which is given as

$$\mathit{\text{M}}=\lfloor 2\pi \widetilde{\mathit{\text{R}}}\rfloor ,$$

(22)

where $\lfloor \cdot \rfloor$ is the round-down operator. The $m$-th, ( $\left|m\right|\le \mathit{\text{M}}$ ) phase mode is excited by the normalized beamforming vector in terms of

$${\mathbf{b}}_m^H=\frac{1}{\mathit{\text{N}}}[e^{jm{\alpha }_1},e^{jm{\alpha }_2},\mathrm{\cdots },e^{jm{\alpha }_N}].$$

(23)

The resulting far field beam pattern of the UCA for mode $m$ is

	$f_m(\phi )$	$={\mathbf{b}}_m^H\mathbf{a}(\phi )$		(24)
		$={\displaystyle \frac{1}{\mathit{\text{N}}}}{\displaystyle \sum _{n=1}^N}{e}^{jm{\alpha }_n}{e}^{j2\pi \widetilde{\mathit{\text{R}}}\cos ({\phi }_d-{\alpha }_n)}$
		$=j^m{J}_m(2\pi \widetilde{\mathit{\text{R}}})e^{jm\phi }$
		$+{\displaystyle \sum _{c=1}^{\mathrm{\infty }}}\left[j^p{J}_p(2\pi \widetilde{\mathit{\text{R}}})e^{-jp\phi }+j^q{J}_q(2\pi \widetilde{\mathit{\text{R}}})e^{-jq\phi }\right],$

where $p=cN-m$ and $q=cN+m$. In order to make the item $j^m{J}_m(2\pi \widetilde{\mathit{\text{R}}})e^{jm\phi }$ of (24) be the dominant one, the number of antenna markedasmathitalicN should meet the following condition

$$N>2M.$$

(25)

The property $J_{-m}(2\pi \widetilde{\mathit{\text{R}}})={(-1)}^m{J}_m(2\pi \widetilde{\mathit{\text{R}}})$ of Bessel functions is used, and the residual items are omitted [10], then the far field beam pattern UCA for mode $m$ can be expressed as

$$f_m(\phi )\approx j^{\left|m\right|}{J}_{\left|m\right|}(2\pi \widetilde{\mathit{\text{R}}})e^{jm\phi }\mathit{\hspace{1em}}\left|m\right|\le \mathit{\text{M}}.$$

(26)

The phase mode excitation matrix is defined as

$${\mathbf{W}}_B^H={\mathbf{C}}_J{\mathbf{B}}^H,$$

(27)

where

${\mathbf{C}}_J$	$=\mathrm{diag}\{j^{-M},\mathrm{\cdots },j^{-1},j^0,j^{-1},\mathrm{\cdots },j^{-M}\}$	(28a)
$\mathbf{B}$	$=\sqrt{N}[{\mathbf{b}}_{-M},\mathbf{\cdots },{\mathbf{b}}_0,\mathbf{\cdots },{\mathbf{b}}_M].$	(28b)

Then the resulting beamspace steering vector synthesized by ${\mathbf{W}}_B$ is given as

${\mathbf{a}}_B(\phi )$

$={\mathbf{W}}_B^H\mathbf{a}(\phi )=\left[\begin{array}{c}\hfill J_M(2\pi \widetilde{\mathit{\text{R}}})e^{-jM\phi }\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill J_1(2\pi \widetilde{\mathit{\text{R}}})e^{-j\phi }\hfill \\ \hfill J_0(2\pi \widetilde{\mathit{\text{R}}})\hfill \\ \hfill J_1(2\pi \widetilde{\mathit{\text{R}}})e^{j\phi }\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill J_M(2\pi \widetilde{\mathit{\text{R}}})e^{jM\phi }\hfill \end{array}\right].$

(29)

The above relation can be rewritten as

$${\mathbf{a}}_B(\phi )={\mathbf{Ja}}_v(\phi ),$$

(30)

where

	$\mathbf{J}=\mathrm{diag}\{J_M(2\pi \widetilde{\mathit{\text{R}}}),\mathrm{\cdots },J_1(2\pi \widetilde{\mathit{\text{R}}}),$		(31)
	$J_0(2\pi \widetilde{\mathit{\text{R}}}),J_1(2\pi \widetilde{\mathit{\text{R}}}),\mathrm{\cdots },J_M(2\pi \widetilde{\mathit{\text{R}}})\},$

is the diagonal matrix of Bessel functions, and

$${\mathbf{a}}_v(\phi )={[e^{-jM\phi },\mathrm{\cdots },e^{-j\phi },1,e^{j\phi },\mathrm{\cdots },e^{jM\phi }]}^T$$

(32)

is the Vandermonde structured steering vector of the virtual linear array.

To sum up, the total transforming matrix is defined as

$${\mathbf{T}}_B={\mathbf{J}}^{-1}{\mathbf{W}}_B^H,$$

(33)

and the transformed $k$-th snapshot vector is

	${\mathbf{y}}_v(k)$	$={\mathbf{T}}_B\mathbf{x}(k)$		(34)
		$={\mathbf{A}}_v\mathbf{s}(k)+{\mathbf{n}}_v(k),$

where ${\mathbf{A}}_v=[{\mathbf{a}}_v({\phi }_1),{\mathbf{a}}_v({\phi }_2),\mathbf{\cdots },{\mathbf{a}}_v({\phi }_D)]$ is the virtual manifold matrix with Vandermonde structure, and ${\mathbf{n}}_v(k)={\mathbf{T}}_B\mathbf{n}(k)$ is the transformed noise vector. The $K$ transformed snapshots can be packed into a matrix with the following form

$${\mathbf{Y}}_v=[{\mathbf{y}}_v(1),{\mathbf{y}}_v(2),\mathbf{\cdots },{\mathbf{y}}_v(K)]\in {ℂ}^{N_v\times K},$$

(35)

where $N_v=2M+1$ is the number of elements in the virtual linear array.

The sample covariance matrix ${\widehat{\mathbf{R}}}_v$ is calculated as

$${\widehat{\mathbf{R}}}_v=\frac{1}{K}{\mathbf{Y}}_v{\mathbf{Y}}_v^H.$$

(36)

However, the transformed noise vector ${\mathbf{n}}_v(k)$ is no longer a white noise vector; in this case, the signal and noise subspace can be solved based on the following Generalized Eigenvalue Decomposition (GEVD) problem:

$${\widehat{\mathbf{R}}}_v\mathbf{u}=\lambda ({\mathbf{T}}_B{\mathbf{T}}_B^H)\mathbf{u}.$$

(37)

Afterward, sort the $N_v$ generalized eigenvalues in descending order; the signal subspace is formed by the generalized eigenvectors corresponding to the first$D_p=D-D_c+1$ ($D_c$ is the number of coherent sources.) large eigenvalues. The signal subspace ${\mathbf{U}}_s$, noise subspace ${\mathbf{U}}_n$ and corresponding generalized eigenvalues ${\mathbf{\Lambda }}_s$, ${\mathbf{\Lambda }}_n$ have the following forms

${\mathbf{U}}_s$	$=[{\mathbf{u}}_1,\mathbf{\cdots },{\mathbf{u}}_{D_p}],$	(38a)
${\mathbf{U}}_n$	$=[{\mathbf{u}}_{D_p+1},\mathbf{\cdots },{\mathbf{u}}_{N_v}],$	(38b)
${\mathbf{\Lambda }}_s$	$=\mathrm{diag}\{{\lambda }_1,\mathrm{\cdots },{\lambda }_{D_p}\},$	(38c)
${\mathbf{\Lambda }}_n$	$=\mathrm{diag}\{{\lambda }_{D_p+1},\mathrm{\cdots },{\lambda }_{N_v}\}.$	(38d)

B. MODE algorithm

In this subsection, we briefly introduced the core idea of the MODE algorithm [5, 21, 22]. The MODE algorithm is an efficient implementation of ML type algorithm which explores the ARMA structure of the observed snapshot vector. The ML algorithm estimates DOA by minimizing the following cost function:

$$f({\mathbf{A}}_v)=\mathrm{tr}\{{\mathbf{\Pi }}_{{\mathbf{A}}_v}^{⟂}{\widehat{\mathbf{R}}}_v\},$$

(39)

where

	${\mathbf{\Pi }}_{{\mathbf{A}}_v}^{⟂}$	$={\mathbf{I}}_{N_v}-{\mathbf{\Pi }}_{{\mathbf{A}}_v}$		(40)
		$={\mathbf{I}}_{N_v}-{\mathbf{A}}_v{({\mathbf{A}}_v^H{\mathbf{A}}_v)}^{-1}{\mathbf{A}}_v^H,$

is the orthogonal projection of ${\mathbf{A}}_v$. Moreover, the cost function in (39) can be replaced with a parameter vector

$$\mathbf{g}={[g_1,g_2,\mathrm{\cdots },g_D]}^T\in {ℂ}^D,$$

(41)

and the connection between $\mathbf{g}$ and DOAs ${\{{\phi }_d\}}_{d=1}^D$ in terms of the following polynomial

$F(z)$	$=g_0{z}^D+g_1{z}^{D-1}+\mathrm{\cdots }+g_D$	(42a)
$=g_0{\displaystyle \prod _{d=1}^D}(z-e^{j{\phi }_d})=0.$	(42b)

Then a matrix $\mathbf{G}\in {ℂ}^{N_v\times (N_v-D)}$ is defined as

$\mathbf{G}=\left[\begin{array}{cccc}\hfill g_D^{\ast }\hfill & \hfill 0\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill \\ \hfill \mathrm{⋮}\hfill & \hfill g_D^{\ast }\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill g_0^{\ast }\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\ddots }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill g_0^{\ast }\hfill & \hfill \mathrm{⋮}\hfill & \hfill g_D^{\ast }\hfill \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill 0\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & \hfill g_0^{\ast }\hfill \end{array}\right],$

(43)

and a link between $\mathbf{G}$ and ${\mathbf{A}}_v$ is established as follows

$${\mathbf{\Pi }}_{{\mathbf{A}}_v}^{⟂}={\mathbf{\Pi }}_{\mathbf{G}}=\mathbf{G}{({\mathbf{G}}^H\mathbf{G})}^{-1}{\mathbf{G}}^H.$$

(44)

By substituting (44) into (39), a reparameterized cost function is cast as

$$f(\mathbf{g})=\mathrm{tr}\{{\mathbf{\Pi }}_{\mathbf{G}}{\mathbf{U}}_s\widehat{\mathbf{\Gamma }}{\mathbf{U}}_s^H\},$$

(45)

where

$\widehat{\mathbf{\Gamma }}$	$=\mathrm{diag}\{{\widehat{\gamma }}_1,\mathrm{\cdots },{\widehat{\gamma }}_{D_p}\},$	(46a)
${\widehat{\gamma }}_d$	$={\displaystyle \frac{{({\lambda }_d-{\widehat{\sigma }}_n^2)}^2}{{\lambda }_d}},$	(46b)
${\widehat{\sigma }}_n^2$	$={\displaystyle \frac{1}{N_v-D_p}}\mathrm{tr}({\mathbf{\Lambda }}_n).$	(46c)

III. PROPOSED ALGORITHM

Comparing the cost function (39) and (45), a complicated searching ${\mathbf{A}}_v$ problem is reduced to an efficient searching $\mathbf{g}$ problem based on the MODE method. However, additional assumptions on $\mathbf{g}$ such as conjugate symmetry is requested by the original MODE algorithm. This additional assumption causes the MODE algorithm less robust to the number of sources. Fortunately, in literatures [26, 27], it has been proved that the solution of cost function (45) is equivalent to the solution of the following weighted least square (WLS) problem:

$$\mathbf{g}=\arg \underset{\mathbf{g}}{\min }{(\mathbf{Fg}-\mathbf{h})}^H\widehat{\mathbf{W}}(\mathbf{Fg}-\mathbf{h}),$$

(47)

where $\mathbf{F}$, $\mathbf{h}$ and the weighting matrix $\widehat{\mathbf{W}}$ are given as

$$\mathbf{F}=\left[\begin{array}{c}\hfill {\mathbf{F}}_1\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathbf{F}}_d\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathbf{F}}_{D_p}\hfill \end{array}\right],\mathbf{h}=\left[\begin{array}{c}\hfill {\mathbf{h}}_1\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathbf{h}}_d\hfill \\ \hfill \mathrm{⋮}\hfill \\ \hfill {\mathbf{h}}_{D_p}\hfill \end{array}\right],$$

(48)

and ${\mathbf{F}}_d\in {ℂ}^{(N_v-D)\times D}$ is the $d$-th submatrix in $\mathbf{F}$:

${\mathbf{F}}_d=\left[\begin{array}{cccc}\hfill {({\mathbf{u}}_d)}_D\hfill & \hfill {({\mathbf{u}}_d)}_{D-1}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill {({\mathbf{u}}_d)}_1\hfill \\ \hfill {({\mathbf{u}}_d)}_{D+1}\hfill & \hfill {({\mathbf{u}}_d)}_D\hfill & \hfill \mathrm{\cdots }\hfill & \hfill {({\mathbf{u}}_d)}_2\hfill \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill {({\mathbf{u}}_d)}_{N_v-1}\hfill & \hfill {({\mathbf{u}}_d)}_{N_v-2}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill {({\mathbf{u}}_d)}_{N_v-D}\hfill \end{array}\right],$

(49)

${\mathbf{h}}_d=-{[{({\mathbf{u}}_d)}_{D+1},\mathrm{\cdots },{({\mathbf{u}}_d)}_{N_v}]}^T\in {ℂ}^{N_v-D},$

(50)

$$\widehat{\mathbf{W}}=\widehat{\mathbf{\Gamma }}\mathbf{\otimes }{({\mathbf{G}}^H\mathbf{G})}^{-1}.$$

(51)

Moreover, the solution to (47) is given by

$$\widehat{\mathbf{g}}={({\mathbf{F}}^H\widehat{\mathbf{W}}\mathbf{F})}^{-1}{\mathbf{F}}^H\widehat{\mathbf{W}}\mathbf{h}.$$

(52)

As we can see, there are no additional requirements in (52), and it is easy to update $\widehat{\mathbf{W}}$ and $\widehat{\mathbf{g}}$ iteratively. The DOAs can be calculated by

$${\widehat{\phi }}_d=\mathrm{\angle }{\widehat{z}}_d,$$

(53)

where ${\widehat{z}}_d$ is the $d$-th root of the polynomial (42). The detailed steps of the proposed algorithm are summarized in Algorithm 1.

---

Algorithm 1 UCA-PUMA Based DOA Estimation

---

Require:

Transformed Snapshots ${\mathbf{Y}}_v$,

Source Number $D$,

Coherent Source Number $D_c$,

Maximum Number of Iterations $N_{\mathrm{Iter}}$,

Tolerance $ϵ$;

Ensure:

Estimated DOAs $\{{\widehat{\phi }}_1,{\widehat{\phi }}_2,\mathrm{\cdots },{\widehat{\phi }}_D\}$

1: Calculate ${\widehat{\mathbf{R}}}_v$ and its GEVD via (36) and (37);

2: Calculate $\mathbf{F}$ and $\mathbf{h}$ via (48)-(50);

3: Initialize $\widehat{\mathbf{g}}$ as ${\widehat{\mathbf{g}}}_0={({\mathbf{FF}}^H)}^{-1}{\mathbf{F}}^H\mathbf{h}$;

4: for $i=1,2,\mathrm{\cdots },{\mathbf{N}}_{\mathrm{Iter}}$ do

5: ${\widehat{\mathbf{W}}}_{i+1}=\widehat{\mathbf{\Gamma }}\mathbf{\otimes }{({\mathbf{G}}_i^H{\mathbf{G}}_i)}^{-1}$;

6: ${\widehat{\mathbf{g}}}_{i+1}={({\mathbf{F}}^H{\widehat{\mathbf{W}}}_{i+1}\mathbf{F})}^{-1}{\mathbf{F}}^H{\widehat{\mathbf{W}}}_{i+1}\mathbf{h}$;

7: ${\mathbf{G}}_{i+1}$ is formed via (43);

8: if ${\parallel {\widehat{\mathbf{g}}}_{i+1}-{\widehat{\mathbf{g}}}_i\parallel }_2/{\parallel {\widehat{\mathbf{g}}}_i\parallel }_2<ϵ$ then

9: break

10: end if

11: end for

12: Calulate DOAs with ${\widehat{\mathbf{g}}}_{i+1}$ based on (53).

---

Figure 3: Source number D = 1.

IV. SIMULATION RESULTS

In this section, we evaluate the performance of the proposed UCA-PUMA algorithm in different scenarios by numerical simulations. The root-mean-square error (RMSE) is adopted to evaluate the estimated DOAs, which is defined as

$$\mathbf{RMSE}=\sqrt{\frac{1}{PD}\sum _{p=1}^P\sum _{d=1}^D{({\widehat{\phi }}_d(p)-{\phi }_d(p))}^2},$$

(54)

where $P=500$ is the number of Monte-Carlo trials.

Figure 4: Source number D = 2 (coherent sources).

Figure 5: Source number D = 3 (coherent sources).

Figure 6: Source number D = 3 (hybrid sources).

A. Single source

In this simulation, we evaluate the performance of the UCA-PUMA algorithm with respect to different SNRs and numbers of snapshots $K$ in the single source ($D=1$) scenario. The DOA of the single source is randomly selected from $\phi \in [0^{\circ },{360}^{\circ })$ in each Monte-Carlo trail. The array parameters are $N=7$, $R/\lambda =0.5$, and $M=3$. The UCA-MODE, UCA-Root-MUSIC, UCA-ESPRIT algorithms, and Cramér–Rao Lower Bound (CRLB) are compared with the UCA-PUMA algorithm.

The RMSE performance versus SNR (Snapshots K=128) and snapshots (SNR=10 dB) are plotted in Fig. 3 (a) and Fig. 3 (b), respectively. As shown in Fig.3, the proposed UCA-PUMA algorithm performs better than the UCA-Root-MUSIC, the UCA-ESPRIT, and the UCA-MODE. Moreover, the RMSE curve of the proposed UCA-PUMA is closest to CRLB for all SNRs and snapshots being simulated, which proves the superiority of our method in the case of a single source.

B. Coherent sources

In this simulation, we focus on the complicated scenarios of multiple coherent sources. Moreover, the aforementioned classical method of processing coherent signals, namely spatial smoothing [18], is also used as a comparison. In the first example, two coherent sources whose DOAs $({\phi }_1$, ${\phi }_2)$ are randomly selected from $[0^{\circ },{360}^{\circ })$, and the separation between ${\phi }_1$ and ${\phi }_2$ is $\mathrm{\Delta }\phi =|{\phi }_1$-${\phi }_2|\ge 5^{\circ }$. As shown in Fig. 4 (a) and Fig. 4 (b), the RMSE of the UCA-PUMA, the Spatial-Smoothing, and the UCA-MODE are lower than ${10}^{\circ }$, and the RMSE of the UCA-MODE is slightly lower than the Spatial-Smoothing when the SNR is lower than 0 dB The proposed UCA-PUMA algorithm has the best performance among these methods, and converges to a tolerable value (${1.6}^{\circ }$) when SNR is greater than 5 dB Furthermore, we fix the SNR to a moderate value (10 dB), the UCA-PUMA still performs the best under various numbers of snapshots. The UCA-MODE and the Spatial-Smoothing behave similarly to the UCA-PUMA in this scenario, yet both the Spatial-Smoothing and the UCA-MODE algorithms show larger RMSEs over the UCA-PUMA. Meanwhile, the performance of the classical subspace-type algorithms (UCA-ESPRIT and UCA-Root-MUSIC) fail significantly.

In the second example, three fully coherent sources with DOAs being ${\phi }_1,{\phi }_2$ and ${\phi }_3$ are simulated; meanwhile the three DOAs are randomly selected from $[0^{\circ },{360}^{\circ })$ and the angular separation between two DOAs is $\mathrm{\Delta }\phi =|{\phi }_i$-${\phi }_j|\ge 5^{\circ }(i,j\in \{1,2,3\},i\ne j)$. The simulation results are shown in Fig. 5 (a) and Fig. 5 (b), respectively. In this scenario with the odd number of sources, the UCA-MODE gives RMSE more than ${10}^{\circ }$ which is severely degraded compared with the even source number scenario. There is also a decrease in the performance of the Spatial-Smoothing compared to the scenario of two sources. However, the UCA-PUMA algorithm still works well. Moreover, the RMSE of the Spatial-Smoothing converges to $4^{\circ }$ while the UCA-PUMA converges to ${2.2}^{\circ }$, and the reason is that the effective aperture of the UCA is reduced in the process of spatial smoothing.

Furthermore, the results depicted in Fig. 6 (a) and Fig. 6 (b) are in the regime where coherent sources are hybrid with uncorrelated ones (Hybrid Sources Scenario). Concretely, the first two sources are coherent, and the third one is uncorrelated from the previous. The DOAs of these three sources are selected in the same way as above. In this hybrid sources scenario, the RMSE of the UCA-MODE is roughly equivalent to ${20}^{\circ }$, and the RMSE of the Spatial-Smoothing converges to $5^{\circ }$ which is slightly higher than the result in the scenario of three coherent sources. Furthermore, the UCA-ESPRIT along with the UCA-Root-MUSIC is completely failed (RMSE $\approx {100}^{\circ }$). However, the UCA-PUMA gives considerable RMSE (${1.5}^{\circ }$) and performs best among the above algorithms whether in scenarios of fully coherent sources or hybrid sources. This also verifies that the UCA-PUMA algorithm is robust to source numbers and performs better than the UCA-MODE or the Spatial-Smoothing.

V. CONCLUSION

In this paper, a robust algorithm named UCA-PUMA for DOA estimation of coherent sources in UCA has been proposed. In order to take advantage of the ARMA model of the snapshot vectors, the non-Vandermonde structured steering vector of UCA is transformed into a virtual Vandermonde structured steering vector in mode space. After that, the recently developed PUMA algorithm performs well in UCA without spatial smoothing. Simulation results demonstrate that the proposed UCA-PUMA algorithm is robust to various numbers of sources. Moreover, the proposed UCA-PUMA exhibits significantly better performance than the algorithms including the UCA-MODE, the Spatial-Smoothing, the UCA-ESPRIT, and the UCA-Root-MUSIC.

REFERENCES

[1] R. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol. 34, no. 3, pp. 276-280, 1986.

[2] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984-995, 1989.

[3] F.-M. Han and X.-D. Zhang, “An ESPRIT-like algorithm for coherent DOA estimation,” IEEE Antennas and Wireless Propagation Letters, vol. 4, pp. 443-446, 2005.

[4] P. Stoica, R. L. Moses, B. Friedlander, and T. Soderstrom, “Maximum likelihood estimation of the parameters of multiple sinusoids from noisy measurements,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 3, pp. 378-392, 1989.

[5] P. Stoica and K. Sharman, “Maximum likelihood methods for direction-of-arrival estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. 38, no. 7, pp. 1132-1143, 1990.

[6] W. Xie, F. Wen, J. Liu, and Q. Wan, “Source association, DOA, and fading coefficients estimation for multipath signals,” IEEE Transactions on Signal Processing, vol. 65, no. 11, pp. 2773-2786, 2017.

[7] W. Du and R. L. Kirlin, “Improved spatial smoothing techniques for DOA estimation of coherent signals,” IEEE Transactions on signal processing, vol. 39, no. 5, pp. 1208-1210, 1991.

[8] C. Qi, Y. Wang, Y. Zhang, and Y. Han, ‘‘Spatial difference smoothing for DOA estimation of coherent signals,” IEEE Signal Processing Letters, vol. 12, no. 11, pp. 800-802, 2005.

[9] L. Ni, J. Zhang, and H. Chen, “Sparse construction decorrelation algorithm of uniform circular array,” in Proceedings of the 2020 4th International Conference on Digital Signal Processing, pp. 263-267, 2020.

[10] C. Mathews and M. Zoltowski, “Eigenstructure techniques for 2-D angle estimation with uniform circular arrays,” IEEE Trans. Signal Process., vol. 42, no. 9, pp. 2395-2407, 1994.

[11] Z. Xu, S. Wu, Z. Yu, and X. Guang, “A robust direction of arrival estimation method for uniform circular array,” Sensors, vol. 19, no. 20, p. 4427, 2019.

[12] Y. Ma, X. Wang, and X. Cao, “Coarray beamspace transformation based DOA estimation for uniform circular arrays,” in 2018 IEEE Radar Conference (RadarConf18), pp. 0792-0797, IEEE, 2018.

[13] D. Zhao, W. Tan, Z. Deng, and G. Li, ‘‘Low complexity sparse beamspace DOA estimation via single measurement vectors for uniform circular array,” EURASIP Journal on Advances in Signal Processing, vol. 2021, no. 1, pp. 1-20, 2021.

[14] T.-J. Shan, M. Wax, and T. Kailath, “On spatial smoothing for direction-of-arrival estimation of coherent signals,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 33, no. 4, pp. 806-811, 1985.

[15] J. Pan, M. Sun, Y. Wang, and X. Zhang, “An enhanced spatial smoothing technique with ESPRIT algorithm for direction of arrival estimation in coherent scenarios,” IEEE Transactions on Signal Processing, vol. 68, pp. 3635-3643, 2020.

[16] Y.-M. Chen, “On spatial smoothing for two-dimensional direction-of-arrival estimation of coherent signals,” IEEE Transactions on Signal Processing, vol. 45, no. 7, pp. 1689-1696, 1997.

[17] M. Wax and J. Sheinvald, “Direction finding of coherent signals via spatial smoothing for uniform circular arrays,” IEEE Transactions on Antennas and Propagation, vol. 42, no. 5, pp. 613-620,1994.

[18] Z. Ye, L. Xiang, and X. Xu, “DOA estimation with circular array via spatial averaging algorithm,” IEEE Antennas and Wireless Propagation Letters, vol. 6, pp. 74-76, 2007.

[19] K. M. Reddy and V. Reddy, “Analysis of spatial smoothing with uniform circular arrays,” IEEE Transactions on Signal Processing, vol. 47, no. 6, pp. 1726-1730, 1999.

[20] P. Stoica, B. Ottersten, M. Viberg, and R. L. Moses, “Maximum likelihood array processing for stochastic coherent sources,” IEEE Transactions on Signal Processing, vol. 44, no. 1, pp. 96-105, 1996.

[21] A. B. Gershman and P. Stoica, “New MODE-based techniques for direction finding with an improved threshold performance,” Signal Processing, vol. 76, no. 3, pp. 221-235, 1999.

[22] C. Ma, Y. Peng, and L. Tian, “An improved MODE algorithm for DOA estimation based on UCA,” in TENCON ’97 Brisbane - Australia. Proceedings of IEEE TENCON ’97. IEEE Region 10 Annual Conference. Speech and Image Technologies for Computing and Telecommunications (Cat. No. 97CH36162), vol. 2, pp. 811-814 vol. 2,1997.

[23] H. C. So, F. K. W. Chan, W. H. Lau, and C.-F. Chan, “An efficient approach for two-dimensional parameter estimation of a single-tone,” IEEE Trans. Signal Process., vol. 58, no. 4, pp. 1999-2009, 2010.

[24] H. Minghao, Y. Yixin, and Z. Xianda, “UCA-ESPRIT algorithm for 2-D angle estimation,” in WCC 2000 - ICSP 2000. 2000 5th International Conference on Signal Processing Proceedings. 16th World Computer Congress 2000, vol. 1, pp. 437-440, 2000.

[25] R. Goossens, H. Rogier, and S. Werbrouck, “UCA Root-MUSIC with sparse uniform circular arrays,” IEEE Trans. Signal Process., vol. 56, no. 8, pp. 4095-4099, 2008.

[26] C. Qian, L. Huang, M. Cao, J. Xie, and H. C. So, “PUMA: An improved realization of MODE for DOA estimation,” IEEE Trans. Aerosp. Electron. Syst., vol. 53, no. 5, pp. 2128-2139, 2017.

[27] C. Qian, L. Huang, N. D. Sidiropoulos, and H. C. So, “Enhanced PUMA for Direction-of-Arrival estimation and its performance analysis,” IEEE Trans. Signal Process., vol. 64, no. 16, pp. 4127-4137, 2016.

BIOGRAPHIES

Ye Tian was born in Liaoning Province, China, in 1996. He received the B.S. degree in Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu, China, in 2017. He is currently pursuing the Ph.D. degree in electromagnetic field and microwave technology with National Space Science Center (Beijing, China) of Chinese Academy of Sciences. His current research interests include array signal processing, DOA estimation and RF signal geolocation.

Yonghui Huang was born in Anshan, China. He received the B.Sc. degree in electronics engineering from Tsinghua University of Beijing, China in 1998. In 2001, he obtained the M.S. degree in aero-spacecraft design from University Chinese Academy of Sciences, China. He achieved the Ph.D. degree in wireless communication from Aalborg University of Aalborg, Denmark in 2008. He is currently a professor with National Space Science Center of Chinese Academy of Science in Beijing, China. From 2002 to 2011, He worked as a Postdoc. and research assistant in Aalborg University of Aalborg, Denmark. He is currently a researcher at the National Space Science Center of the Chinese Academy of Sciences. His current research interests include radio frequency machine learning, phased array antenna and transmitter linearization. He is the TPC member of IEEE CCET and IEEE WiSEE.

Xiaoxu Zhang was born in Fujian Province, China, in 1997. She received the B.S. degree in Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu, China, in 2019. She is currently pursuing the Ph.D. degree in Computer Science and Technology with the National Space Science Center at the Chinese Academy of Sciences, Beijing, China. Her current research interests include signal sorting, RF fingerprint recognition, and deep learning.

Meiyan Lin was born in Fujian Province, China, in 1995. She received the B.S. degree in School of Electronic and Information Engineering from Beijing Jiaotong University, Beijing, China, in 2018. She is currently pursuing the Ph.D. degree in National Space Science Center at Chinese Academy of Sciences, Beijing, China. Her main research interests include statistical signal processing and emitter localization.

ACES JOURNAL, Vol. 37, No. 6, 692–701.
doi: 10.13052/2022.ACES.J.370605
© 2021 River Publishers

Abstract

I. INTRODUCTION

II. SYSTEM MODEL

A. UCA mode space transformation

B. MODE algorithm

III. PROPOSED ALGORITHM

IV. SIMULATION RESULTS

A. Single source

B. Coherent sources

V. CONCLUSION

REFERENCES

BIOGRAPHIES