An Efficient Direct Position Determination Method for Multiple Strictly Noncircular Sources

Abstract

This paper focuses on the localization methods for multiple sources received by widely separated arrays. The conventional two-step methods extract measurement parameters and then estimate the positions from them. In the contrast to the conventional two-step methods, direct position determination (DPD) localizes transmitters directly from original sensor outputs without estimating intermediate parameters, resulting in higher location accuracy and avoiding the data association. Existing subspace data fusion (SDF)-based DPD developed in the frequency domain is computationally attractive in the presence of multiple transmitters, whereas it does not use special properties of signals. This paper proposes an improved SDF-based DPD algorithm for strictly noncircular sources. We first derive the property of strictly noncircular signals in the frequency domain. On this basis, the observed frequency-domain vectors at all arrays are concatenated and extended by exploiting the noncircular property, producing extended noise subspaces. Fusing the extended noise subspaces of all frequency components and then performing a unitary transformation, we obtain a cost function for each source location, which is formulated as the smallest eigenvalue of a real-valued matrix. To avoid the exhaustive grid search and solve this nonlinear function efficiently, we devise a Newton-type iterative method using matrix Eigen-perturbation theory. Simulation results demonstrate that the proposed DPD using Newton-type iteration substantially reduces the running time, and its performance is superior to other localization methods for both near-field and far-field noncircular sources.

1. Introduction

Noncircular complex signals (e.g., binary-phase-shift-keying (BPSK), multiple-amplitude-shift-keying (MASK) and offset-quadrature-phase-shift-keying (OQPSK) modulated signals) are extensively employed in modern communication systems. During the last few years, the algorithms for noncircular complex sources with application to signal processing [1,2,3,4,5,6,7] have received an upsurge of attention. Different from circular signals, the unconjugated statistical property of noncircular complex signals adds to the amount of available information, which can help to enhance the performance. Since many communication systems offer location-based services, the problem of accurate localization for noncircular sources is worth further investigation.

Conventional localization methods employ two-step processing [8,9,10,11], where the measurement parameters (e.g., direction of arrival (DOA), time of arrival (TOA) and frequency difference of arrival (FDOA)) are first extracted, and then, the source positions are estimated. In the first step, parameter estimation algorithms for noncircular (NC) signals such as the NC-MUSIC algorithm (multiple signal classification algorithm for noncircular signals) [5] can improve the accuracy, due to the exploitation of the noncircular property. However, the localization algorithm in the second step cannot utilize the properties of signals. Moreover, when multiple sources exist, the two-step localization system confronts the association problem of deciding which of the multiple measurements reported by the observers corresponds to which source. If the measurements are not correctly related to the transmitter, additional errors are produced. As a consequence, these conventional two-step methods are suboptimal and cannot guarantee high localization precision for noncircular sources.

Compared with the conventional two-step methods, direct position determination (DPD) [12,13,14,15,16,17,18,19,20,21,22] is a single-step localization method without computing the intermediate parameters and augments the position estimation with the constraint that all measurements correspond to the same geolocation of the transmitter. Therefore, DPD can avoid the association problem, and its location accuracy has shown to be significantly higher than that of the conventional two-step methods, especially under low signal-to-noise ratio (SNR) conditions [13]. However, the DPD technique often requires the transmission of the sampling data to a central processing location and yields a large amount of computations [16].

Nowadays, DPD algorithms applied to the scenario of widely-separated arrays have been intensively investigated [13,14,15,16,17,18]. A maximum likelihood (ML)-based DPD algorithm was presented in [16], where the locations of multiple sources are decoupled into several lower dimensional optimization problems with the information of uncorrelated waveforms. However, in passive localization systems, signal waveforms are always unknown at the receiver. When solving ML estimators for multiple sources without known waveforms, there is a variety of stray parameters, which requires a substantial computational effort. To this end, Amar and Weiss proposed an iterative algorithm for DPD of multiple unknown signals to reduce the complexity [17]. Although this ML estimator can approach the associated Cramér–Rao bound, its iterative procedure is complicated as it has to perform a grid search during each iteration.

To make the localization computationally attractive for multiple transmitters, an alternative method is the subspace data fusion (SDF). It has lower computational complexity than ML methods, due to the fact that all source positions are estimated from a MUSIC-like cost function that depends on the parameters of the same dimension as that for only one source. Two SDF-based DPDs were developed in [16] and [18], respectively. One is based on the time domain and implicitly uses the array responses [18], and the other processes the frequency-domain observations and relies on the assumption that the envelopes of the signals are the same at all observers, up to delay and amplitude caused by the propagation channel [16]. Therefore, the latter can exploit the location information embedded in both array responses and TOAs, leading to higher accuracy. However, these SDF-based DPD algorithms were designed for general sources (i.e., circular signals) and did not consider the property of noncircular signals. For this reason, we proposed an SDF-based DPD for strictly noncircular sources observed by a moving array in our early work [19], which exploits the time-domain property of noncircular signals. It can be easily adapted to the scenario of widely-separated arrays. Note that similar to the SDF-based DPD in [18], the DPD developed in our previous work utilizes only the array responses, neglecting the location information in the propagation time between the transmitter and the observer.

In light of the aforementioned related works, in this study, we consider the scenario of widely-separated arrays and assume a line of sight (LOS) propagation of multiple signals with unknown waveforms and unknown complex attenuation at each observer array. The purpose of this study is to develop an efficient SDF-based DPD method in the frequency domain for noncircular sources. The contributions of this paper can be summarized as follows:

• We derive the frequency-domain property of strictly noncircular signals and propose an improved SDF estimator. Based on this noncircular property, we establish an extended frequency-domain observation received by all arrays and compute the extended subspaces, which are implicitly related to array responses and TOAs. Fusing the extended subspaces of all frequency components, a cost function is formulated as the smallest eigenvalue of a symmetric real-valued matrix for each source location, due to a unitary transformation. Therefore, the real-valued eigen-decomposition is required instead of complex computations. Compared with the primitive SDF-based DPDs, this improved SDF estimator retains its superiority for requiring low-dimensional optimization and has higher robustness to noise that comes from exploiting noncircularity.
• We devise a Newton-type iterative algorithm to efficiently solve the prescribed cost function based on matrix Eigen-perturbation theory. It substantially reduces the computations of the straightforward implementation of the optimization for each position, which is always accomplished via a two- or three-dimensional grid search.

The remainder of this paper is organized as follows. Section 2 presents the signal model and formulates the problem. In Section 3, we propose the extended SDF DPD estimator for noncircular sources and devise a Newton-type iterative solution. In Section 4, we show the simulation results and make the discussion. Section 5 concludes the study.

Throughout the paper, upper case and lower case boldface letters will represent matrices and column vectors, respectively. For convenience, we list the notations used in this paper:

${\{\cdot \}}^{*}$	Conjugate.
${\{\cdot \}}^{\mathrm{T}}$	Transpose.
${\{\cdot \}}^{\mathrm{H}}$	Conjugate transpose.
$\mathrm{blkdiag}\{\cdot \}$	Composition of the block diagonal matrix.
$\mathrm{diag}\{\cdot \}$	Composition of the diagonal matrix.
$\mathrm{vec}\{\cdot \}$	The “vectorization” operator that turns a matrix into a vector by stacking the columns of the matrix, one below another.
$\otimes$	Kronecker matrix product.
$\mathbb{E}[\cdot ]$	Expectation.
$\mathrm{tr}\{\cdot \}$	Trace.
$\mathrm{Re}\{\cdot \}$	Real part.
$\mathrm{Im}\{\cdot \}$	Imaginary part.
${[\cdot ]}_n$	The n-th element of a vector.
${[\cdot ]}_{nm}$	The n,m-th entry of a matrix.
${\Vert \cdot \Vert }_2$	Euclidean norm.
${ℂ}^{N\times M}$	Set of the $N\times M$ complex matrices.
${ℝ}^{N\times M}$	Set of the $N\times M$ real matrices.

2. Problem Formulation

2.1. Property of the Noncircular Signal

A signal is normally circular when an arbitrary rotation cannot change its first-order and second-order statistical characteristics, namely rotation invariance. If a signal’s property of rotation invariance is not satisfied, then we take it as a noncircular signal. For a straightforward description, let $s(n)$ be a zero-mean signal, $s(n)$ is considered to be circular with the condition “$\mathbb{E}[|s(n){|}^2]\ne 0$ and $\mathbb{E}[s^2(n)]=0$”, whereas $s(n)$ is referred to as being noncircular with the condition “$\mathbb{E}[|s(n{)|}^2]\ne 0$ and $\mathbb{E}[s^2(n)]\ne 0$”. Specifically, a strictly noncircular signal can be generated through a phase shift of a real-valued signal [23]. According to this property, we assume $s(n)$ to be a strictly noncircular signal, and therefore, it can be generated as [23]:

$$s(n)=\stackrel{⌣}{s}(n)e^{\mathrm{i}\varphi },$$

(1)

where $\stackrel{⌣}{s}(n)$ is a real-valued signal and $\varphi$ signifies the initial phase of $s(n)$.

Using this time-domain property of a strictly noncircular signal, we can express the discrete Fourier transform (DFT) of $s(n)$ as follows:

$$\begin{array}{ll}\overline{s}(j)& ={\displaystyle \sum _{n=0}^{J-1}s(n)\exp (-\mathrm{i}2\mathsf{\pi }nj/J)}\\ & ={\displaystyle \sum _{n=0}^{J-1}\stackrel{⌣}{s}(n)e^{\mathrm{i}\varphi }\exp (-\mathrm{i}2\mathsf{\pi }nj/J)}j=0,1,\dots ,J-1,\end{array}$$

(2)

where $\overline{s}(j)$ is the q-th DFT coefficient of $s(n)$. The conjugate of $\overline{s}(j)$ is given by:

$${\overline{s}}^{*}(j)={\displaystyle \sum _{n=0}^{J-1}\stackrel{⌣}{s}(n)e^{-\mathrm{i}\varphi }\exp (\mathrm{i}2\mathsf{\pi }nj/J)}0,1,\dots ,J-1.$$

(3)

Using $\exp (\mathrm{i}2\mathsf{\pi }nj/J)=\exp (-\mathrm{i}2\mathsf{\pi }n(J-j)/J)$, we rewrite (3) by:

$$\begin{array}{ll}{\overline{s}}^{*}(j)& ={\displaystyle \sum _{n=0}^{J-1}\stackrel{⌣}{s}(n)e^{-\mathrm{i}\varphi }\exp (-\mathrm{i}2\mathsf{\pi }n(J-j)/J)}\\ & =e^{-\mathrm{i}2\varphi }{\displaystyle \sum _{n=0}^{J-1}\stackrel{⌣}{s}(n)e^{\mathrm{i}\varphi }\exp (-\mathrm{i}2\mathsf{\pi }n(J-j)/J)}j=0,1,\dots ,J-1.\end{array}$$

(4)

Based on (2), we have:

$$\overline{s}(J-j)={\displaystyle \sum _{n=0}^{J-1}\stackrel{⌣}{s}(n)e^{\mathrm{i}\varphi }\exp (-\mathrm{i}2\mathsf{\pi }n(J-j)/J)}.$$

(5)

Combing (4) and (5), it can be checked that:

$${\overline{s}}^{*}(j)=\overline{s}(J-j)e^{-\mathrm{i}2\varphi }j=0,1,\dots ,J-1.$$

(6)

The equation related to the frequency-domain signals as shown in (6) lays the foundation of subsequent theoretical development.

2.2. Frequency-Domain Signal Model and Problem Formulation

Let us consider $L$ stationary observers, each of whom is equipped with an antenna array comprising $M$ isotropic sensors. $Q$ transmitters are assumed to radiate uncorrelated narrowband and strictly noncircular signals in the field of interest. We denote ${\mathit{p}}_q\in {ℝ}^{D\times 1}$ as the position of the $q$-th transmitter for $q=1,2,\dots ,Q$. The complex envelope of the outputs of the $l$-th observer array at time $0\le t\le T$ is expressed as [16,17]:

$${\mathit{r}}_l(t)={\displaystyle \sum _{q=1}^Q{b}_{ql}{\mathit{a}}_l({\mathit{p}}_q)s_q(t-{\tau }_l({\mathit{p}}_q)-t_{q0})}+{\mathit{n}}_l(t)$$

(7)

for $l=1,2,\dots ,L.$ Here, $s_q(t-{\tau }_l({\mathit{p}}_q)-t_{q0})$ is the complex envelope of the $q$-th source signal, transmitted at time $t_{q0}$ and delayed by ${\tau }_l({\mathit{p}}_q)$, which indicates the propagation time between the $q$-th transmitter and the $l$-th observer. We assume that the source signals are wide-sense stationary. $b_{ql}$ is the unknown complex attenuation coefficient representing the channel effect between the $q$-th transmitter and the $l$-th observer. ${\mathit{a}}_l({\mathit{p}}_q)$ is the spatial array response of the $l$-th observer array to the signal transmitted from the position ${\mathit{p}}_q$, which is related to the corresponding DOA. ${\mathit{n}}_l(t)$ is the white and circularly Gaussian noise mixed through the $l$-th observer array. The sources and noises are assumed to be uncorrelated and to have a mean of zero.

Similarly to [16], we partition ${\mathit{r}}_l(t)(0\le t\le T)$ into $K$ sections, and each length equals $T/K$. Then, sample the signals of each section at time intervals $j{T}_s(j=0,1,\dots ,J-1)$ where $T_s$ represents the sampling period. Thus, the DFT of the $J$ samples in the $k$-th section can be expressed as [16]:

$${\overline{\mathit{r}}}_l(j,k)={\displaystyle \sum _{q=1}^Q{b}_{ql}{\mathit{a}}_l({\mathit{p}}_q){\overline{s}}_q(j,k)\exp (-\mathrm{i}2\mathsf{\pi }{\tau }_l({\mathit{p}}_q)j/J{T}_s)}+{\overline{\mathit{n}}}_l(j,k)j=0,1,\dots ,J-1,$$

(8)

where ${\overline{s}}_q(j,k)$ and ${\overline{\mathit{n}}}_l(j,k)$ are the $j$-th DFT coefficients of $s_q(t-t_{q0})$ and ${\mathit{n}}_l(t)$ in the $k$-th section, respectively. For convenience of presentation, we define:

$${\overline{a}}_l(j,{\mathit{p}}_q)={\mathit{a}}_l({\mathit{p}}_q)\exp (-\mathrm{i}2\mathsf{\pi }{\tau }_l({\mathit{p}}_q)j/J{T}_s),$$

(9)

and rewrite ${\overline{\mathit{r}}}_l(j,k)$ by:

$${\overline{\mathit{r}}}_l(j,k)={\displaystyle \sum _{q=1}^Q{b}_{ql}{\overline{\mathit{a}}}_l(j,{\mathit{p}}_q){\overline{s}}_q(j,k)}+{\overline{\mathit{n}}}_l(j,k)j=0,1,\dots ,J-1.$$

(10)

We now composite the j-th DFT coefficients of all the outputs of $L$ observer arrays and obtain:

$$\begin{array}{ll}\overline{r}(j,k)& ={[{\overline{r}}_1^{\mathrm{T}}(j,k),{\overline{r}}_2^{\mathrm{T}}(j,k),\dots ,{\overline{r}}_L^{\mathrm{T}}(j,k)]}^{\mathrm{T}}\\ & ={\displaystyle \sum _{q=1}^Q\Gamma (j,{\mathit{p}}_q)b_q{\overline{s}}_q(j,k)}+\overline{n}(j,k)\\ & ={\displaystyle \sum _{q=1}^Q{\alpha }_q(j){\overline{s}}_q(j,k)}+\overline{n}(j,k),\end{array}$$

(11)

where:

$$\overline{n}(j,k)={[{\overline{n}}_1^{\mathrm{T}}(j,k),{\overline{n}}_2^{\mathrm{T}}(j,k),\dots ,{\overline{n}}_L^{\mathrm{T}}(j,k)]}^{\mathrm{T}},$$

(12)

$$\Gamma (j,{\mathit{p}}_q)=\mathrm{blkdiag}\{{\overline{a}}_1(j,{\mathit{p}}_q),{\overline{a}}_2(j,{\mathit{p}}_q),\dots ,{\overline{a}}_L(j,{\mathit{p}}_q)\},$$

(13)

$$b_q={[b_{q1},b_{q2},\dots ,b_{qL}]}^{\mathrm{T}},$$

(14)

$${\alpha }_q(j)=\Gamma (j,{\mathit{p}}_q)b_q.$$

(15)

Let us denote $\Phi (j)=[{\alpha }_1(j),{\alpha }_2(j),\dots ,{\alpha }_Q(j)]\in {ℂ}^{LM\times Q}$ and $\overline{s}(j,k)={[{\overline{s}}_1(j,k),{\overline{s}}_2(j,k),\dots ,{\overline{s}}_Q(j,k)]}^{\mathrm{T}}$. Then, $\overline{r}(j,k)$ can be further expressed in the matrix form as:

$$\overline{r}(j,k)=\Phi (j)\overline{\mathit{s}}(j,k)+\overline{n}(j,k).$$

(16)

Therefore, the covariance matrix of $\overline{r}(j,k)$ can be written as:

$${\mathit{R}}_{\mathit{r}}(j)=\mathbb{E}[\overline{r}(j,k){\overline{r}}^{\mathrm{H}}(j,k)]=\Phi (j){\mathit{R}}_s(j){\Phi }^{\mathrm{H}}(j)+{\sigma }_{\mathrm{n}}^2{I}_{LM}.$$

(17)

where ${\mathit{R}}_s(j)=\mathbb{E}[\overline{s}(j,k){\overline{s}}^{\mathrm{H}}(j,k)]$, ${\sigma }_{\mathrm{n}}^2$ denotes the noise power and $I_{LM}$ signifies the $LM\times LM$ identity matrix.

Note that $\Phi (j)$ is associated with both DOAs and TOAs, and therefore, it contains the location information. The SDF DPD based on the frequency domain [16] relies on the general signal model as shown in (17), which cannot take full advantage of the signal properties. To this end, we come to exploit the noncircularity of signals derived in Section 2.1.

As the conjugate of $\overline{r}(j,k)$ can be expressed as:

$${\overline{\mathit{r}}}^{*}(j,k)={\Phi }^{*}(j){\overline{s}}^{*}(j,k)+{\overline{\mathit{n}}}^{*}(j,k),$$

(18)

applying (6) to (18) leads to:

$${\overline{\mathit{r}}}^{*}(j,k)={\Phi }^{*}(j){\Delta }_{\varphi }\overline{s}(J-j,k)+{\overline{\mathit{n}}}^{*}(j,k),$$

(19)

where ${\Delta }_{\varphi }=\mathrm{diag}\{e^{-\mathrm{i}2{\varphi }_1},e^{-\mathrm{i}2{\varphi }_2},\dots ,e^{-\mathrm{i}2{\varphi }_Q}\}$ with ${\varphi }_q$ being the initial phase of the q-th strictly noncircular source. Replacing $j$ with $J-j$ in (19), it follows that:

$${\overline{\mathit{r}}}^{*}(J-j,k)={\Phi }^{*}(J-j){\Delta }_{\varphi }\overline{s}(j,k)+{\overline{\mathit{n}}}^{*}(J-j,k),$$

(20)

for $j=0,1,\dots ,J-1$. Note that ${\overline{\mathit{r}}}^{*}(J,k)={\overline{\mathit{r}}}^{*}(0,k)$ as the DFT coefficient is periodic with the period $J$.

By exploiting the unconjugated property of noncircular signals, we combine (16) and (20), and therefore construct the extended observation:

$$\tilde{r}(j,k)={[{\overline{r}}^{\mathrm{T}}(j,k),{\overline{r}}^{\mathrm{H}}(J-j,k)]}^{\mathrm{T}}=\tilde{\Phi }(j)\overline{s}(j,k)+\tilde{n}(j,k),$$

(21)

where $\tilde{n}(j,k)={[{\overline{n}}^{\mathrm{T}}(j,k),{\overline{n}}^{\mathrm{H}}(J-j,k)]}^{\mathrm{T}}$ and $\tilde{\Phi }(j)={[{\Phi }^{\mathrm{T}}(j),{\Delta }_{\varphi }^{\mathrm{T}}{\Phi }^{\mathrm{H}}(J-j)]}^{\mathrm{T}}\in {ℂ}^{2LM\times Q}$, whose q-th column is given by:

$${\tilde{\alpha }}_q(j)=\tilde{\Gamma }(j,{\mathit{p}}_q){\tilde{b}}_q,$$

(22)

with:

$$\tilde{\Gamma }(j,{\mathit{p}}_q)=\mathrm{blkdiag}\left\{\Gamma (j,{\mathit{p}}_q),{\Gamma }^{*}(J-j,{\mathit{p}}_q)\right\},$$

(23)

$${\tilde{b}}_q={[b_q^{\mathrm{T}},b_q^{\mathrm{H}}{e}^{-\mathrm{i}2{\varphi }_q}]}^{\mathrm{T}}.$$

(24)

Then, the covariance matrix of the extended observation $\tilde{r}(j,k)$ has the following form:

$${\mathit{r}}_{\tilde{r}}(j)=\mathbb{E}[\tilde{r}(j,k){\tilde{r}}^{\mathrm{H}}(j,k)]=\tilde{\Phi }(j){\mathit{R}}_s(j){\tilde{\Phi }}^{\mathrm{H}}(j)+{\sigma }_{\mathrm{n}}^2{I}_{2LM}.$$

(25)

Until now, we have obtained the frequency-domain data model for our work as shown in (25). It is noteworthy that the dimension of the extended covariance matrix ${\mathit{R}}_{\tilde{r}}(j)$ is twice that of the covariance matrix ${\mathit{R}}_{\mathit{r}}(j)$ for circular signals, leading to more available information. To this end, we will apply ${\mathit{R}}_{\tilde{r}}(j)$ to our problem. The problem that we address now is, given ${\mathit{R}}_{\tilde{r}}(j)$ for $j=0,1,\dots ,J-1$, to directly estimate the locations of multiple transmitters without explicitly computing TOAs and DOAs.

3. Methods

3.1. Extended SDF

Based on (25), the eigen-decomposition of ${\mathit{R}}_{\tilde{r}}(j)$ yields:

$${\mathit{R}}_{\tilde{r}}(j)=U_k^{\mathrm{s}}(j){\Sigma }_k^{\mathrm{s}}(j)U_k^{\mathrm{sH}}(j)+{\sigma }_{\mathrm{n}}^2{U}_k^{\mathrm{n}}(j)U_k^{\mathrm{nH}}(j),$$

(26)

where ${\Sigma }_k^{\mathrm{s}}(j)$ is a diagonal matrix with diagonal elements of the $Q$ largest eigenvalues of ${\mathit{R}}_{\tilde{r}}(j)$. $U_k^{\mathrm{s}}(j)\in {ℂ}^{2LM\times Q}$ is the extended signal subspace comprising the eigenvectors corresponding to the $Q$ largest eigenvalues, and it spans the same space as $\tilde{\Phi }(j)$ [5,16,19]. $U_k^{\mathrm{n}}(j)\in {ℂ}^{2LM\times (2LM-Q)}$ is the extended noise subspace consisting of the eigenvectors corresponding to the $2LM-Q$ smallest eigenvalue ${\sigma }_{\mathrm{n}}^2$.

According to the subspace orthogonality principle [5,16,19], namely the columns of $\tilde{\Phi }(j)$ are orthogonal to the noise subspace of ${\mathit{R}}_{\tilde{r}}(j)$, we have:

$${\tilde{\alpha }}_q^{\mathrm{H}}(j)U_k^{\mathrm{n}}(j)U_k^{\mathrm{nH}}(j){\tilde{\alpha }}_q(j)=0j=0,1,\dots ,J-1$$

(27)

for $q=1,2,\dots ,Q$. Fusing $J$ DFT components leads to:

$${\displaystyle \sum _{j=0}^{J-1}{\tilde{\alpha }}_q^{\mathrm{H}}(j)U_k^{\mathrm{n}}(j)U_k^{\mathrm{nH}}(j){\tilde{\alpha }}_q(j)}=0.$$

(28)

Inserting (22) into (28), it follows that:

$${\tilde{b}}_q^{\mathrm{H}}Q({\mathit{p}}_q){\tilde{b}}_q=0,$$

(29)

where $Q({\mathit{p}}_q)$ is expressed as:

$$Q({\mathit{p}}_q)={\displaystyle \sum _{j=0}^{J-1}{\tilde{\Gamma }}^{\mathrm{H}}(j,{\mathit{p}}_q)U_k^{\mathrm{n}}(j)U_k^{\mathrm{nH}}(j)\tilde{\Gamma }(j,{\mathit{p}}_q)}\in {ℂ}^{2L\times 2L}.$$

(30)

It is noteworthy that ${\tilde{b}}_q$ in (29) is formulated as ${\tilde{b}}_q={[b_q^{\mathrm{T}},b_q^{\mathrm{H}}{e}^{-\mathrm{i}2{\varphi }_q}]}^{\mathrm{T}}$. Next, we will utilize this structure of ${\tilde{b}}_q$ to transform the complex-valued equation in (29) into a real-valued one.

We now rewrite (29) by:

$${({\tilde{b}}_q{e}^{\mathrm{i}{\varphi }_q})}^{\mathrm{H}}Q({\mathit{p}}_q){\tilde{b}}_q{e}^{\mathrm{i}{\varphi }_q}={{\tilde{b}}^{\prime }}_q^{\mathrm{H}}Q({\mathit{p}}_q){{\tilde{b}}^{\prime }}_q=0,$$

(31)

where ${{\tilde{b}}^{\prime }}_q={\tilde{b}}_q{e}^{\mathrm{i}{\varphi }_q}$. As ${\tilde{b}}_q={[b_q^{\mathrm{T}},b_q^{\mathrm{H}}{e}^{-\mathrm{i}2{\varphi }_q}]}^{\mathrm{T}}$, ${{\tilde{b}}^{\prime }}_q$ can be given by:

$${{\tilde{b}}^{\prime }}_q={[{b^{\prime }}_q^{\mathrm{T}},{b^{\prime }}_q^{\mathrm{H}}]}^{\mathrm{T}}$$

(32)

with ${b^{\prime }}_q=b_q{e}^{\mathrm{i}{\varphi }_q}.$ Based on (32), we perform a unitary transformation of ${{\tilde{b}}^{\prime }}_q$ to establish a real-valued vector:

$${\stackrel{⌣}{b}}_q=\frac{1}{\sqrt{2}}T{{\tilde{b}}^{\prime }}_q,$$

(33)

where $T$ is the unitary matrix:

$$T=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_L& I_L\\ -\mathrm{i}{I}_L& \mathrm{i}{I}_L\end{array}\right],$$

(34)

and thus, ${\stackrel{⌣}{b}}_q$ can be expressed in the form of (35):

$${\stackrel{⌣}{b}}_q={\left[\mathrm{Re}\{{b^{\prime }}_q^{\mathrm{T}}\},\mathrm{Im}\{{b^{\prime }}_q^{\mathrm{T}}\}\right]}^{\mathrm{T}}.$$

(35)

Since $T^{\mathrm{H}}T=I_{2L}$, it follows that:

$${{\tilde{b}}^{\prime }}_q=\sqrt{2}{T}^{\mathrm{H}}{\stackrel{⌣}{b}}_q.$$

(36)

Substituting (36) back into (31), we get:

$${\stackrel{⌣}{b}}_q^{\mathrm{H}}{Q}^{\prime }({\mathit{p}}_q){\stackrel{⌣}{b}}_q=0,$$

(37)

where:

$$Q^{\prime }({\mathit{p}}_q)=TQ({\mathit{p}}_q)T^{\mathrm{H}}\in {ℂ}^{2L\times 2L}.$$

(38)

Given that ${\stackrel{⌣}{b}}_q$ is real-valued, we rewrite the left side of the equation in (37) by:

$$\begin{array}{l}{\stackrel{⌣}{b}}_q^{\mathrm{H}}{Q}^{\prime }({\mathit{p}}_q){\stackrel{⌣}{b}}_q\\ ={\stackrel{⌣}{b}}_q^{\mathrm{T}}(\mathrm{Re}\{Q^{\prime }({\mathit{p}}_q)\}+\mathrm{iIm}\{Q^{\prime }({\mathit{p}}_q)\}){\stackrel{⌣}{b}}_q\\ ={\stackrel{⌣}{b}}_q^{\mathrm{T}}\mathrm{Re}\{Q^{\prime }({\mathit{p}}_q)\}{\stackrel{⌣}{b}}_q+\mathrm{i}{\stackrel{⌣}{b}}_q^{\mathrm{T}}\mathrm{Im}\{Q^{\prime }({\mathit{p}}_q)\}{\stackrel{⌣}{b}}_q.\end{array}$$

(39)

Due to the fact that ${\stackrel{⌣}{b}}_q^{\mathrm{H}}{Q}^{\prime }({\mathit{p}}_q){\stackrel{⌣}{b}}_q$ is real-valued, we discard the imaginary term, $\mathrm{i}{\stackrel{⌣}{b}}_q^{\mathrm{T}}\mathrm{Im}\{Q^{\prime }({\mathit{p}}_q)\}{\stackrel{⌣}{b}}_q$, in (39). Then, ${\stackrel{⌣}{b}}_q^{\mathrm{H}}{Q}^{\prime }({\mathit{p}}_q){\stackrel{⌣}{b}}_q$ is equivalent to:

$${\stackrel{⌣}{b}}_q^{\mathrm{H}}{Q}^{\prime }({\mathit{p}}_q){\stackrel{⌣}{b}}_q={\stackrel{⌣}{b}}_q^{\mathrm{T}}\mathrm{Re}\{Q^{\prime }({\mathit{p}}_q)\}{\stackrel{⌣}{b}}_q.$$

(40)

According to the above derivation, the complex-valued equation in (29) is transformed into the following real-valued equation:

$${\stackrel{⌣}{b}}_q^{\mathrm{T}}\mathrm{Re}\{Q^{\prime }({\mathit{p}}_q)\}{\stackrel{⌣}{b}}_q=0.$$

(41)

We now replace $U_k^{\mathrm{n}}(j)$ with ${\widehat{U}}_k^{\mathrm{n}}(j)$, which is computed by the eigen-decomposition of the estimated ${\mathit{R}}_{\tilde{r}}(j)$:

$${\widehat{R}}_{\tilde{r}}(j)=\frac{1}{K}{\displaystyle \sum _{k=1}^K\tilde{r}(j,k){\tilde{r}}^{\mathrm{H}}(j,k)}.$$

(42)

Then, the following optimization model presents the solution for ${\mathit{p}}_q$ and ${\stackrel{⌣}{b}}_q$:

$$\{{\widehat{p}}_q,{\widehat{\stackrel{⌣}{b}}}_q\}=\underset{p,\stackrel{⌣}{b}}{\mathrm{argmin}}{\stackrel{⌣}{b}}^{\mathrm{T}}\mathrm{Re}\{Q^{\prime }(p)\}\stackrel{⌣}{b}.$$

(43)

We observe that the estimation of ${\stackrel{⌣}{b}}_q$ is the eigenvector corresponding to the smallest eigenvalue of $\mathrm{Re}\{Q^{\prime }(p)\}$, and the estimation of ${\mathit{p}}_q$ is the point that minimizes the smallest eigenvalue of $\mathrm{Re}\{Q^{\prime }(p)\}$. Consequently, the position of the q-th transmitter can be obtained by:

$${\widehat{p}}_q=\underset{\mathit{p}}{\mathrm{argmin}}f(p)$$

(44)

for $q=1,2,\dots ,Q$, where:

$$f(p)={\lambda }_{\min }\left\{\mathrm{Re}\{Q^{\prime }(p)\}\right\}.$$

(45)

Here, ${\lambda }_{\min }\{\cdot \}$ signifies the smallest eigenvalue of the input matrix. As $f(p)$ is a nonlinear function in terms of the source position, the minimization of $f(p)$ is usually accomplished by a $D$-dimensional grid search.

3.2. Iterative Solution

When the area of interest is large and the grid step size is small, the grid search is computationally demanding. To locate the transmitters quickly, we will devise a Newton-type iterative method to efficiently solve the proposed cost function for each source location. As is known to all, the Newton-type iterative method provides an efficient way to solve nonlinear optimization problems [24]. It requires computing the gradient and Hessian matrix of the cost function, which is quite easy if the cost function is explicitly written. However, the cost function in (45) is the smallest eigenvalue of a symmetric matrix, whose gradient and Hessian matrix cannot be obtained straightforwardly. Therefore, we apply the matrix Eigen-perturbation theory to derive the gradient and Hessian matrix of (45) for each source position. For this purpose, a proposition needs to be introduced first.

Proposition: Assume a symmetric real-valued matrix $X\in {ℝ}^{N\times N}$ with ${\lambda }_{\min }={\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_N={\lambda }_{\max }$ being its eigenvalues and $e_1,e_2,\dots ,e_N$ being the corresponding eigenvectors. Denote $\delta X$ as the perturbation of $X$, and hence $\widehat{X}=X+\delta X$ where $\delta X$ and $\widehat{X}$ are also symmetric. Then, the eigenvalues of the perturbed matrix $\widehat{X}$, ${\widehat{\lambda }}_{\min }={\widehat{\lambda }}_1\le {\widehat{\lambda }}_2\le \cdots \le {\widehat{\lambda }}_N={\widehat{\lambda }}_{\max }$, can be expressed as [25]:

$${\widehat{\lambda }}_n={\lambda }_n+e_n^{\mathrm{T}}\delta X{e}_n+e_n^{\mathrm{T}}\delta X{E}_n\delta X{e}_n+o\left({\Vert \delta X\Vert }_2^2\right)(n=1,2,\dots ,N),$$

(46)

where:

$$E_n={\displaystyle \sum _{\begin{array}{c}m=1\\ m\ne n\end{array}}^N{({\lambda }_n-{\lambda }_m)}^{-1}{e}_m{e}_m^{\mathrm{T}}}(n=1,2,\dots ,N),$$

(47)

and $o\left({\Vert \delta X\Vert }_2^2\right)$ signifies the infinitesimal term of ${\Vert \delta X\Vert }_2^2$. This proposition reveals the relationship between the matrix perturbation and the eigenvalues of the disturbed matrix. The proof of (46) can be found in [25].

Now, let us denote $X(p)=\mathrm{Re}\{Q^{\prime }(p)\}$, whose eigenvalues and eigenvectors are ${\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_{2L}$ and $e_1,e_2,\dots ,e_{2L}$. Define $\delta p$ as the perturbation of $p$ and $\widehat{p}=p+\delta p$ as the perturbed position. Then, the eigenvalues and eigenvectors of $\widehat{X}(\widehat{p})$ are ${\widehat{\lambda }}_1\le {\widehat{\lambda }}_2\le \cdots \le {\widehat{\lambda }}_{2L}$ and ${\widehat{e}}_1,{\widehat{e}}_2,\dots ,{\widehat{e}}_{2L}$, respectively. Denoting $\delta X=\widehat{X}(\widehat{p})-X(p)$, we obtain the following equation according to the above proposition:

$${\widehat{\lambda }}_1={\lambda }_1+e_1^{\mathrm{T}}\delta X{e}_1+e_1^{\mathrm{T}}\delta X{E}_1\delta X{e}_1+o\left({\Vert \delta X\Vert }_2^2\right),$$

(48)

in which:

$$E_1={\displaystyle \sum _{m=2}^{2L}{({\lambda }_1-{\lambda }_m)}^{-1}{e}_m{e}_m^{\mathrm{T}}}.$$

(49)

Applying $\mathrm{vec}\{ABC\}=(C^{\mathrm{T}}\otimes A)\mathrm{vec}\{B\}$ and $\mathrm{tr}\{ABCD\}={\mathrm{vec}}^{\mathrm{T}}\{D^{\mathrm{T}}\}(C^{\mathrm{T}}\otimes A)\mathrm{vec}\{B\}$ (see [26]), (48) can be further expressed as a function of $\delta x=\mathrm{vec}\{\delta X\}$:

$$\begin{array}{ll}{\widehat{\lambda }}_1& ={\lambda }_1+\mathrm{vec}\{e_1^{\mathrm{T}}\delta X{e}_1\}+\mathrm{tr}\{e_1{e}_1^{\mathrm{T}}\delta X{E}_1\delta X\}+o\left({\Vert \delta X\Vert }_2^2\right)\\ & ={\lambda }_1+(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}})\delta x+\delta x^{\mathrm{T}}(E_1\otimes e_1{e}_1^{\mathrm{T}})\delta x+o\left({\Vert \delta x\Vert }_2^2\right).\end{array}$$

(50)

where the symmetry of $E_1$ and $\delta X$ is used.

Note that ${\widehat{\lambda }}_1$ is the cost function of the proposed estimator, and therefore, our purpose is to link $\delta p$ to ${\widehat{\lambda }}_1$. As (50) displays the relationship between $\delta x$ and ${\widehat{\lambda }}_1$, we proceed to derive the relationship between $\delta x$ and $\delta p$ in what follows.

Denote $x(p)=\mathrm{vec}\{X(p)\}$ and $\widehat{x}(\widehat{p})=\mathrm{vec}\{\widehat{X}(\widehat{p})\}$. Based on the result in [27], we can take the second-order Taylor expansion of $\widehat{x}(\widehat{p})$ around $p$ as:

$$\widehat{x}(\widehat{p})=x(p)+{\dot{X}}_{\mathit{p}}(p)\delta p+\frac{1}{2}(I_{4L^2}\otimes \delta {\mathit{p}}^{\mathrm{T}}){\ddot{X}}_{pp}(p)\delta p+o\left({\Vert \delta p\Vert }_2^2\right),$$

(51)

where:

$${\dot{X}}_{\mathit{p}}(p)=\frac{\partial x(p)}{\partial {\mathit{p}}^{\mathrm{T}}}\in {ℝ}^{4L^2\times D}$$

(52)

and:

$${\ddot{X}}_{pp}(p)=\left[\begin{array}{c}\frac{{\partial }^2{[x(p)]}_1}{\partial p\partial p^{\mathrm{T}}}\\ \frac{{\partial }^2{[x(p)]}_2}{\partial p\partial p^{\mathrm{T}}}\\ ⋮\\ \frac{{\partial }^2{[x(p)]}_{4L^2}}{\partial p\partial p^{\mathrm{T}}}\end{array}\right]\in {ℝ}^{4L^{2}D\times D}.$$

(53)

For the further derivations of (52) and (53), see Appendix A and Appendix B, respectively.

As $\delta x=\widehat{x}(\widehat{p})-x(p)$, according to (51), we can relate $\delta x$ with $\delta p$:

$$\delta x={\dot{X}}_{\mathit{p}}(p)\delta p+\frac{1}{2}(I_{4L^2}\otimes \delta {\mathit{p}}^{\mathrm{T}}){\ddot{X}}_{pp}(p)\delta p+o\left({\Vert \delta p\Vert }_2^2\right).$$

(54)

Then, substituting (54) back into (50) yields:

$$\begin{array}{ll}{\widehat{\lambda }}_1& ={\lambda }_1+(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}})\delta x+\delta x^{\mathrm{T}}(E_1\otimes e_1{e}_1^{\mathrm{T}})\delta x+o\left({\Vert \delta x\Vert }_2^2\right)\\ & ={\lambda }_1+(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}}){\dot{X}}_{\mathit{p}}(p)\delta p+\\ & \frac{1}{2}(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}})(I_{4L^2}\otimes \delta {\mathit{p}}^{\mathrm{T}}){\ddot{X}}_{pp}(p)\delta p+\delta {\mathit{p}}^{\mathrm{T}}{\dot{X}}_{\mathit{p}}^{\mathrm{T}}(p)(E_1\otimes e_1{e}_1^{\mathrm{T}}){\dot{X}}_{\mathit{p}}(p)\delta p+o\left({\Vert \delta p\Vert }_2^2\right).\end{array}$$

(55)

Because:

$$(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}})(I_{4L^2}\otimes \delta {\mathit{p}}^{\mathrm{T}}){\ddot{X}}_{pp}(p)\delta p=\delta {\mathit{p}}^{\mathrm{T}}(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}}\otimes I_D){\ddot{X}}_{pp}(p)\delta p,$$

(56)

(55) is equivalent to:

$$\begin{array}{ll}{\widehat{\lambda }}_1& ={\lambda }_1+(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}}){\dot{X}}_{\mathit{p}}(p)\delta p+\frac{1}{2}\delta {\mathit{p}}^{\mathrm{T}}\left\{(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}}\otimes I_D){\ddot{X}}_{pp}(p)+2{\dot{X}}_{\mathit{p}}^{\mathrm{T}}(p)(E_1\otimes e_1{e}_1^{\mathrm{T}}){\dot{X}}_{\mathit{p}}(p)\right\}\delta p+o\left({\Vert \delta p\Vert }_2^2\right)\\ & ={\lambda }_1+g^{\mathrm{T}}(p)\delta p+\frac{1}{2}\delta {\mathit{p}}^{\mathrm{T}}H(p)\delta p+o\left({\Vert \delta p\Vert }_2^2\right),\end{array}$$

(57)

where $g(p)$ and $H(p)$ are given by:

$$g(p)={\dot{X}}_{\mathit{p}}^{\mathrm{T}}(p)(e_1\otimes e_1)$$

(58)

and:

$$H(p)=(e_1^{\mathrm{T}}\otimes e_1^{\mathrm{T}}\otimes I_D){\ddot{X}}_{pp}(p)+2{\dot{X}}_{\mathit{p}}^{\mathrm{T}}(p)(E_1\otimes e_1{e}_1^{\mathrm{T}}){\dot{X}}_{\mathit{p}}(p).$$

(59)

Finally, applying the result in (57) to the objective function in (45), we can express $f(\widehat{p})$ in terms of the second-order perturbation of $\delta p$ as follows:

$$f(\widehat{p})=f(p)+g^{\mathrm{T}}(p)\delta p+\frac{1}{2}\delta {\mathit{p}}^{\mathrm{T}}H(p)\delta p+o\left({\Vert \delta p\Vert }_2^2\right).$$

(60)

It can be seen from (60) that $g(p)$ and $H(p)$ are the gradient and Hessian matrix of $f(p)$, respectively. Therefore, we can get the estimation of ${\mathit{p}}_q$ following the Newton-type iterative procedure [28]:

$${\widehat{p}}_q^{(i+1)}={\widehat{p}}_q^{(i)}-{\mu }^i{H}^{-1}({\widehat{p}}_q^{(i)})g({\widehat{p}}_q^{(i)}),$$

(61)

where the superscript $i$ denotes the iteration number and $\mu$ is the step factor deciding the iteration step size. In this way, each transmitter can be located when an initial position is given. Note that the variable step size is used to solve the contradiction between convergence speed and misadjustment, which is usually encountered in the case of fixed step size.

To summarize this iterative process, we provide the procedure below:

Method: Eigen-Perturbation-Based Newton-Type Iterative Method

1: Set e as a positive number that is small enough. Then, choose a suitable $K\in {\mathbb{Z}}^{+}$ and $J\in {\mathbb{Z}}^{+}$ ${\mathbb{Z}}^{+}$ denotes the set of positive integers), and compute the DFT of the observed samples. 2: Use (42) to estimate the extended covariance matrix ${\widehat{R}}_{\tilde{r}}(j)$ for $j=0,1,\dots ,J-1.$. 3: Compute the eigenvectors corresponding to the $2LM-Q$ smallest eigenvalues of ${\widehat{R}}_{\tilde{r}}(j)$ for $j=0,1,\dots ,J-1$ and thus obtain the extended noise subspaces ${\widehat{U}}_k^{\mathrm{n}}(j)$ for $j=0,1,\dots ,J-1.$ 4: for q = 1,2,…,Q do 5: Initialize ${\widehat{p}}_q^{(0)}$, and make $i\leftarrow 0.$ 6: Insert ${\widehat{p}}_q^{(i)}$ into (58) and (59) to compute $g({\widehat{p}}_q^{(i)})$ and $H({\widehat{p}}_q^{(i)})$. 7: Compute ${\Delta }_i={\Vert g({\widehat{p}}_q^{(i)})\Vert }_2$. If ${\Delta }_i\le e,{\widehat{p}}_q^{(i)}$ is the estimated source position, and stop iterating; otherwise, continue. 8: Use (61) to update ${\widehat{p}}_q^{(i+1)}$. Make $i\leftarrow i+1$, and go to Step 6.

9: end for

Remark 1. Similarly to existing SDF-based DPDs, the method proposed above has an advantage over ML-based methods in the presence of multiple transmitters as it estimates the positions in a decoupled manner. Different from ML functions, which have only one global extremum, our cost function exhibits a few minima at the estimated source positions. Considering that the Newton-type iteration is a local search method and is bound to converge to a local minimum, we use it as a fast converging search procedure instead of the fine search. The initial guess will determine the position to which source the converging local minimum corresponds. Therefore, our iterative method for solving the proposed cost function relies on the initialization to resolve multiple sources, whereas the optimization methods for solving the ML functions require a good initial guess to avoid the local minima. For the initialization of our iterative method, the initial position of each transmitter can be effectively obtained through a coarse search as long as it can separate multiple sources. This can guarantee the accurate and fast localization.

3.3. Computational Complexity

The main computational costs of the proposed methods using exhaustive grid search and using Newton-type iteration and the SDF DPD method for circular sources based on the frequency domain [16] include four parts, i.e., computing the DFT, estimating the covariance matrix, performing the eigen-decomposition and solving the cost function.

Table 1. Computational complexity.

Method	Complexity				Comments
	Computing DFT	Estimating Covariance Matrix	Eigen-Decomposition	Solving Cost Function
SDF-based DPD in [16]	$O(2LMKJ{\log }_{2}J)$	$O(4L^2{M}^{2}KJ)$	$O(4L^3{M}^{3}J)$	$O((4L^2{M}^{2}QJ+6L^2{M}^{2}J+8L^3)Q{N}_p)$	- For $Q$ circular sources
Proposed DPD using exhaustive search	$O(2LMKJ{\log }_{2}J)$	$O(16L^2{M}^{2}KJ)$	$O(32L^3{M}^{3}J)$	$O\left(\left(\begin{array}{l}8L^3(M^2+M)J-\\ 8L^2(M+1)QJ+8L^3\end{array}\right)Q{N}_p\right)$	- For $Q$ strictly noncircular sources
Proposed DPD using Newton-type iteration	$O(2LMKJ{\log }_{2}J)$	$O(16L^2{M}^{2}KJ)$	$O(32L^3{M}^{3}J)$	$O\left(\left(\begin{array}{l}(8L^3(M^2+M)J-\\ 8L^2(M+1)QJ){(D+1)}^2+\\ 4L^2{D}^2+8L^3+4L^{2}D+\\ 4L^{}{D}^2+D^3+2LD\end{array}\right)Q{N}_{iter}\right)$	- For $Q$ strictly noncircular sources - Complexity for initialization is not included

lists the order of the number of real-valued multiply operations in each part consumed by the aforementioned methods, where $N_{\mathit{p}}$ denotes the number of position grid points and $N_{iter}$ represents the iteration number required for convergence.

As shown in Table 1, the SDF-based DPD in [16] is more computationally efficient than the proposed DPD using the exhaustive search because the dimension of the covariance matrix employed in our algorithm is twice as large as the dimension of the covariance matrix used in [16]. On the other hand, the complexity of the proposed DPD using Newton-type iteration is much lower than that of the exhaustive search implementation since $N_{\mathit{p}}$ is usually much larger than $N_{iter}$. For further comparison, we will examine their computer running times through the simulations described in Section 4.

4. Results

The purpose of this section is to present the simulation results and performance analysis relying on Monte Carlo simulations. We consider the scenario of three observers located at the positions [3,3]${}^{\mathrm{T}}(\mathrm{km})$, [3,−3]${}^{\mathrm{T}}(\mathrm{km})$ and ${[-3,3]}^{\mathrm{T}}(\mathrm{km})$. Each observer is equipped with a uniform linear array (ULA) composed of $M=3$ sensors. The adjacent sensors are spaced with $0.5\lambda$, where $\lambda$ denotes the wavelength. Without loss of generality, we assume two transmitters in the field of interest. They arrive at observers with attenuation coefficient vectors being $b_1=$ [10.94 + 0.34i,0.77 + 0.64i]^T and $b_2=$ [10.87 + 0.5i,0.77 + 0.64i]^T, respectively. The baseband signal waveforms are generated as narrowband strictly noncircular signals with identical power ${\sigma }_{\mathrm{s}}^2$, whose initial phases are ${\varphi }_1=\mathsf{\pi }/5(\mathrm{rad})$ and ${\varphi }_2=\mathsf{\pi }/3(\mathrm{rad})$. The noises are uncorrelated, circularly symmetric Gaussian random variables with constant power ${\sigma }_{\mathrm{n}}^2$. Unless stated otherwise, each observer collects signals in $K=50$ sections each of $J=16$ frequencies using a sample rate of 40k (samples/s). For comparison, we invoke the following five estimators:

• Proposed DPD estimator using the devised Newton-type iterative method.
• Proposed DPD estimator using the exhaustive search.
• SDF-based DPD estimator for general circular sources in the frequency domain [16] (denoted by FD-DPD).
• SDF-based DPD estimator for noncircular sources in time domain [19] (denoted by NC TD-DPD).
• Two-step processing estimator: DOA estimation using the NC-MUSIC algorithm [5] and TOA estimation using the ML criterion for noncircular signals [7] at each observer, and pseudo-linear weighted least square localization with the DOA and TOA estimates from all observers used as the data, where DOA and TOA estimates are assumed to be associated with the correct transmitter (denoted by two-step).

In the simulations, the exhaustive search of DPD estimators is implemented through a coarse search with a 0.01-km resolution and then a fine search with a 0.001-km resolution. For the Newton-type iterative method, the positions are initialized using a coarse search with a 0.1-km resolution, which is capable of resolving sources, as will be shown later on. The step factor is $\mu =0.9$. Under different experimental conditions, each simulation is conducted in 1000 Monte Carlo trials, and the location root mean square error (RMSE) is used to assess the localization performance:

$${\mathrm{RMSE}}_q=\sqrt{\frac{1}{1000}{\displaystyle \sum _{n=1}^{1000}{\Vert p_q-{\widehat{p}}_{n,q}\Vert }^2}}$$

(62)

for the q-th source, where ${\widehat{p}}_{n,q}^{}$ denotes the estimation of $p_q$ in the n-th Monte Carlo trial.

This section will show the performance of the proposed algorithms in seven examples. We begin with an analysis of the initial condition and the convergence performance of our Newton-type iterative method. Two near-field sources are assumed to be located at $p_1={[1,-1]}^{\mathrm{T}}(\mathrm{km})$ and $p_2={[0,2]}^{\mathrm{T}}(\mathrm{km})$, and the corresponding scenario is shown in Figure 1. For different SNRs, we assess the inverse cost function of the proposed estimator over the candidate positions with a 0.1-km resolution. Figure 2 displays the corresponding evaluation of the inverse cost function within a square area of $2\times 2(\mathrm{km}\times \mathrm{km})$ around the true positions of the two transmitters at SNR = 0 dB. It can be seen that the maximum of the inverse cost function appears near the true position of each source, and there is no pseudo peak in this area, which indicates that the coarse search with a 0.1-km resolution can provide the effective initial positions. Then, we evaluate the normalized Euclidean norm of the gradient for each source versus the number of iteration steps, when the SNR is −10 dB, 0 dB and 10 dB, respectively. As shown in Figure 3, fewer than six iterations are required for convergence in the whole range of SNRs.

In the second example, for the simulation conditions given above, we roughly evaluate the computer running time of each of the prescribed location algorithms based on an average of 1000 estimates, when the SNR is 0 dB. The computer machine used is a ThinkPad laptop equipped with a 2.50-GHz Intel Core CPU and 8 GB RAM. The results are listed in

Table 2. Average runtime for near-field sources. FD, frequency domain; TD, time domain; NC, noncircular.

Method	Runtime (s)
Proposed DPD (Newton-type Iterative Method)	0.0818
Proposed DPD (Exhaustive Grid Search)	0.7802
Proposed DPD (Nelder–Mead Simplex Search)	0.1231
FD-DPD	0.5274
NC TD-DPD	0.5070
Two-step	0.3120

. To make our results as comprehensive as possible, we also implement the fine search through the Nelder–Mead simplex search approach [29], which is a kind of local search method, and compare its average running time with that of the proposed Newton-type iterative method. We note that the complexity of our algorithm using the exhaustive grid search is larger than those of other approaches. As expected, our iterative method is the most computationally attractive among the six algorithms, and it accomplishes the location in about two thirds of the time required for the Nelder–Mead simplex search. Compared with the exhaustive search implementation, the running time of our iterative method is reduced by approximately ten times. This result agrees with the analysis in Section 3.3.

The third example will assess the accuracy of our algorithm for the near-field sources as shown in Figure 1. We plot the RMSE curves of the proposed DPD algorithms and other positioning algorithms versus SNR in Figure 4. It is obvious that the proposed DPD algorithms significantly outperform other comparison algorithms at low SNRs. As the SDF-based DPD in the time domain exploits only the DOA information, it exhibits the lowest accuracy even when the SNR is sufficiently high. These results indicate that both the noncircularity and the TOA information are helpful for enhancing accuracy. Furthermore, since the proposed DPD methods and two-step approach utilize the location information embedded in DOAs and TOAs for noncircular sources, the improved localization performance with respect to the traditional two-step method reveals the advantage of the DPD technique. As the RMSE curves of the proposed estimator using the Newton-type iterative method and that using Nelder–Mead simplex search overlap with each other, we provide their RMSE values in

Table 3. The estimated RMSEs of the proposed DPDs using the Newton-type iterative method and Nelder–Mead simplex search for near-field sources (km).

Source	Method	SNR (dB)
		−10	−6	−2	2	6	10
Source 1	Newton-type Iterative Method	0.167	0.079	0.049	0.028	0.017	0.011
	Nelder-Mead Simplex Search	0.166	0.078	0.050	0.029	0.017	0.012
Source 2	Newton-type Iterative Method	0.257	0.109	0.063	0.038	0.024	0.015
	Nelder–Mead Simplex Search	0.258	0.110	0.063	0.039	0.023	0.016

. It can be observed that the difference of the accuracy between two methods is not more than 0.001 km.

In the fourth example, we simulate the location performance with the number of sections $K$ varying from 20 to 100 when the SNR is 0 dB. The location geometry corresponds to that in the foregoing examples. As shown in Figure 5, the performance of our algorithm remains promising in the range of different numbers of sections. When the sample number becomes smaller, this performance improvement turns out to be more prominent.

In the fifth example, we study the location accuracy when the distance between the two transmitters changes. The two transmitters are assumed to be placed at $p_1={[0,-d/2]}^{\mathrm{T}}(\mathrm{km})$ and $p_2={[0,d/2]}^{\mathrm{T}}(\mathrm{km})$, respectively. The SNR is as high as 10 dB. As $d$ increases from 0.8 km to 2 km, we depict the RMSEs of the proposed DPDs and the existing SDF-based DPDs. Since the two-step approach can hardly separate the two transmitters with small distances, we do not show its RMSE curves. It can be seen from Figure 6 that our algorithms performs considerably better than other DPDs when the two sources are closely located. Specifically, the SDF-based DPD in the time domain fails to localize the second transmitter in the case of $d<1.6(\mathrm{km})$. This demonstrates that our algorithm exhibits superior resolution compared to the existing SDF-based DPDs.

In this example, we consider two far-field transmitters located at $p_1={[-12,-20]}^{\mathrm{T}}(\mathrm{km})$ and $p_2={[-20,-10]}^{\mathrm{T}}(\mathrm{km})$ as illustrated in Figure 7. To examine the initial positions of the far-field sources obtained by the coarse search with a 0.1-km resolution, we compute the inverse cost function of the proposed estimator over the candidate positions with a grid step size of 0.1 km for different SNRs. The results for SNR = 0 dB are shown in Figure 8. We note that the resolution of this coarse search is sufficiently high to resolve the far-field sources and thus to obtain the effective initialization.

Finally, we come to examine the running time and location accuracy of the proposed algorithms for the far-field sources as shown in Figure 7. When the SNR is 0 dB, we evaluate the average running time of each location algorithm. The corresponding results are shown in

Table 4. Average runtime for far-field sources.

Method	Runtime (s)
Proposed DPD (Newton-type Iterative Method)	0.0856
Proposed DPD (Exhaustive Grid Search)	0.7869
Proposed DPD (Nelder–Mead Simplex Search)	0.1330
FD-DPD	0.5144
NC TD-DPD	0.5072
Two-step	0.3132

, which confirms that the devised iterative method is less complex than other methods. Especially, in contrast to the implementation of the Nelder–Mead simplex search, our Newton-type iterative method serves as a faster local search procedure. Under the same geometry, Figure 9 displays the location accuracy for the two far-field sources versus SNR, indicating that the proposed algorithms still hold the lowest RMSEs. When the SNR is below 0 dB, the RMSEs of the proposed method are much lower than those of other location algorithms by at least 0.5 km. Moreover, we notice that the performance of the two-step approach (where the DOA and TOA estimations are designed for noncircular signals in the first step) is severely degraded at low SNRs, whereas it performs better than the SDF-based DPD in the frequency domain when the SNR is higher than 5 dB. Comparing this result with that in Figure 4, we find that the exploitation of the noncircular property is more crucial to the performance enhancement for the far-field sources than for the near-field sources. Besides, we show the RMSEs of the proposed estimators using the Newton-type iterative method and using the Nelder–Mead simplex search in

Table 5. The estimated RMSEs of the proposed DPDs using the Newton-type iterative method and Nelder–Mead simplex search for far-field sources (km).

Source	Method	SNR (dB)
		−5	0	5	10	15	20
Source 1	Newton-type Iterative Method	0.935	0.243	0.108	0.053	0.029	0.016
	Nelder-Mead Simplex Search	0.933	0.250	0.101	0.050	0.029	0.016
Source 2	Newton-type Iterative Method	1.168	0.325	0.135	0.068	0.035	0.020
	Nelder–Mead Simplex Search	1.160	0.324	0.141	0.068	0.036	0.021

. It can be seen that the accuracies of these two local search methods have no significant difference for the far-field sources.

Additionally, the results indicate that the location accuracy of the proposed DPD using the Newton-type iterative method is approximately the same as those of the Nelder–Mead simplex search and the exhaustive grid search. Considering that the running time of the Newton-type iterative method is much shorter than those of other methods (Table 2 and Table 4), this iterative method provides an efficient way to solve our DPD problem.

5. Conclusions

This study has investigated an efficient single-step localization method for noncircular sources received by widely-separated arrays. Through deriving and exploiting the potential of strictly noncircular signals in the frequency domain, the proposed algorithm fuses all the extended noise subspaces of all frequency components for the locations of multiple noncircular sources directly without estimating DOAs and TOAs. Relying on a unitary transformation, we have derived a cost function for each transmitter position, which is formulated as the smallest eigenvalue of a symmetric real-valued matrix. To further reduce the computation load of the exhaustive search for solving this cost function, we have devised a Newton-type iterative method based on the matrix Eigen-perturbation theory. Note that this iterative method can be easily extended to solve any estimator whose cost function is the eigenvalue of a specified matrix. Simulation results demonstrate that the proposed location method using Newton-type iteration is computationally efficient. Its running time is approximately one tenth of that of the exhaustive grid search and is nearly two thirds of the running time of the Nelder–Mead simplex search. Furthermore, our method shows greater performance robustness at low SNRs and under poor location geometry, compared with the conventional two-step approach and existing SDF DPD algorithms. Specifically, when the SNR is below 0 dB, the RMSEs of the proposed method are much lower than those of the existing SDF DPDs by at least 0.5 km for the simulated far-field sources.

To further improve the DPD performance and obtain the optimal location accurately and quickly, it is a possible direction to design a hybrid method combining an evolutionary algorithm and a fast local search approach for solving the ML-based DPD problem.