Angular Parameters Estimation of Multiple Incoherently Distributed Sources Generating Noncircular Signals

TL;DR

This paper presents a new, less complex method for accurately estimating the directions and angular spreads of multiple incoherently-distributed non-circular sources, outperforming existing techniques in accuracy and robustness.

Contribution

The paper introduces a novel decoupled estimation method for DOAs and angular spreads of non-circular sources, reducing complexity and improving performance.

Findings

Outperforms state-of-the-art techniques in accuracy

Provides robust estimation in challenging scenarios

Reduces computational complexity

Abstract

We introduce a new method for the estimation of the angular parameters [i.e., central directions of arrival (DOAs) and angular spreads] of multiple non-circular and incoherently-distributed (ID) sources and thoroughly analyze its performance. By decoupling the estimation of the central DOAs from that of the angular spreads, we reduce significantly the complexity of the proposed technique. The latter outperforms most well-known state-of-the-art techniques in terms of estimation accuracy and robustness.

Equations

$$\delta (.)$$

$$(.{)}^{\ast }$$

$$\angle (.)$$

$$|.|$$

$$(.{)}^T$$

$$(.{)}^H$$

$$\simeq$$

$${\text{argmin}}_K\{.\}$$

$$\text{tr}\{\mathbf{A}\}$$

$$\text{diag}\{\mathbf{v}\}$$

$${\Vert .\Vert }_{\text{Fro}}$$

$$\Re \{.\}$$

$$\mathbb{E}\{.\}$$

$$\frac{{\partial }^n(.)}{\partial (.{)}^n}$$

$$\text{eig}\{\mathbf{A}\}$$

$$\text{varodot}$$

$${\mathbf{I}}_p$$

$$0_{p\times q}$$

$$\text{Toeplitz}\{\mathbf{v}\}$$

$$\text{Hankel}\{{\mathbf{v}}_1,{\mathbf{v}}_2\}$$

$$a_l(\theta )$$

$${f_l(\theta )}_{l=0,2,\dots ,(L-1)}$$

$${\bar{\Theta }}_k$$

$${\bar{\sigma }}_k$$

$${\rho }_k(\theta ,{\bar{\mathbf{\psi }}}_k)$$

$$p_{kk}(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k)$$

$$p_{kk}^{\prime }(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k)$$

$$p_{k{k}^{\prime }}(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k,{\bar{\mathbf{\psi }}}_{k^{\prime }})$$

$$p_{k{k}^{\prime }}^{\prime }(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k,{\bar{\mathbf{\psi }}}_{k^{\prime }})$$

$$x_l(n)$$

$$a_l(\theta )$$

$$\bar{\mathbf{\Theta }}$$

$$\bar{\mathbf{\sigma }}$$

$$\mathbf{x}(n)$$

$$\mathbf{w}(n)$$

$$\mathbf{a}(\theta )$$

$$p_{kk}(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k)$$

$$=$$

$$p_{kk}^{\prime }(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k)$$

$$=$$

$$p_{kk}^{\prime }(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k)$$

$${\mathbf{R}}_{\mathbf{xx}}$$

$${\mathbf{R}}_{\mathbf{xx}}^{\prime }$$

$$p_{k{k}^{\prime }}(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k,{\bar{\mathbf{\psi }}}_{k^{\prime }})$$

$$p_{k{k}^{\prime }}^{\prime }(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k,{\bar{\mathbf{\psi }}}_{k^{\prime }})$$

$$p_{k{k}^{\prime }}(\theta ,{\theta }^{\prime };{\bar{\mathbf{\psi }}}_k,{\bar{\mathbf{\psi }}}_{k^{\prime }})$$

Full text

Keywords: Angular spread estimation, central DOAs estimation, multiple incoherently distributed sources, noncircular signals, stochastic Cramér-Rao lower bound (CRLB).

I Introduction

Direction of arrivals estimation for multiple plane waves impinging on an arbitrary array of sensors has received a significant amount of attention over the last several decades [References]. It has typically found many applications in different areas such as modern wireless communication systems [References], audio/speech processing systems [References], radar and sonar [References], just to name a few. In most applications, however, DOA estimation methods are based on the point-source model which postulates that the signals are generated from far-field point sources and travel along a single path to the receiving antenna array. Using this simplified model, many DOA estimators have been developed for both temporally uncorrelated [References-References] and correlated [References,References] signals. However, in real-world surroundings, especially in typical urban environments, multipath propagation made by a cluster of reflections close to each mobile causes angular spreading [References]. In other words, the signal radiated by each source hits the antenna array via different paths with different angles. In this more realistic model, the source is viewed by the array as spatially distributed, i.e., with a central DOA and an angular spread. The latter influences the quality of the communication link and represents an important characteristic for spatial diversity schemes [References, References]. DOA estimation becomes more challenging in presence of local scattering [References, References] because the latter affects the signal spatial distribution. In this context, some studies have shown that classical point-source estimation methods suffer from severe performance degradation when applied to the distributed-source scenario [References, References]. This observation has prompted an increasing interest, over the few recent years, in developing DOA estimation algorithms that can handle both point and scattered sources in order to improve direction finding capabilities in real-world propagation environments.

Depending on the nature of scattering, signal components arriving from different directions exhibit varying degrees of correlation. Hence, we distinguish two different types of the propagation channel. The first one is when the received signal components originated from a source and scattered at different angles are delayed and scaled replicas of the same signal. This feature is known in the literature as “coherent source distribution” or “coherently-distributed (CD) source” [References]. The second type of the propagation channel corresponds to the fact that the signal components of a source impinging from different scatterers at different angles are uncorrelated. This is termed in the literature as “incoherent source distribution” or “incoherently-distributed (ID) source” [References, References]. Therefore, for uncorrelated CD sources, each source contributes rank-one component to the spatial covariance matrix and, as such, the rank of the noise-free covariance matrix is equal to the number of sources [References]. Consequently, many classical DOA estimation methods based on the simplistic point-source model can be easily extended to CD sources. Particulary, authors proposed in [References] an efficient DSPE algorithm for estimating the angular parameters of CD sources. This method enables a decoupled estimation of the DOAs from that of the angular spreads of sources with small angular spread. However, for ID sources, the whole observation space is occupied by signal components, and the noise subspace is generally degenerate [References]. Therefore, the rank of the noise-free covariance matrix is different from the number of sources; it even increases with the angular spread. This makes the trivial generalization of traditional point-source subspace-based methods to the ID case not feasible. To sidestep this problem, tremendous efforts have been directed to developing new angular parameters estimators that are specifically tailored to ID sources. In particular, techniques that are able to handle a single ID source were developed in [References-References].

Many estimators were also developed to estimate the angular parameters of multiple ID sources. In fact, a class of subspace methods were proposed in [References, References, References] wherein the effective dimension of the signal subspace is defined as the number of the first eigenvalues (of the noise-free covariance matrix) that reflect most of the signal energy. More computationally attractive approaches that are based on the beamforming techniques were later introduced in [References, References]. Despite their good performance, all these estimators assume the angular distributions to be perfectly known and identical to all the sources. Methods which are able to handle the multi-source case with known but different angular distributions were also proposed in [References, References]. Recently, a robust version of the generalized Capon principle [References] (RGC) has been developed in [References] which, in contrast to all existing approaches, does not need the a priori knowledge of the angular distributions. Moreover, the latter does not need to be the same for all the sources. This robust approach is, however, statistically less efficient than the aforementioned subspace-based (high-resolution) methods [References, References, References], especially in the presence of closely-spaced ID sources.

More recently, A. Zoubir et al. proposed an efficient subspace-based (ESB) algorithm [References] to estimate the angular parameters of multiple ID circular sources. ESB enjoys a good trade-off between estimation performance and computational complexity. In order to alleviate the computational burden stemming from eigendecomposing the covariance matrix, ESB exploits the properties of its inverse and estimates the angular parameters using a 2-D search. Both the statistical efficiency and high-resolution capabilities of the subspace-based techniques are maintained and, most interestingly, ESB is not limited to a particular antenna array geometry or to a specific type of scatterers’ angular distribution. Yet, it still requires the angular distribution to be perfectly known and identical for all the sources on the top of being derived specifically for circular sources. Recently, a new method for tracking the central DOAs assuming multiple ID mobile sources has been also proposed in [References]. It is based on a simple covariance fitting optimization technique [References] to estimate the central DOAs and the Kalman filter to model the dynamic property of directional changes for the moving sources. Despite its efficiency, this method requires the sources’ angular distributions to be perfectly known and is derived for circular sources only.

Noncircular signals, however, such as binary-phase-shift-keying (BPSK) and offset quadrature-phase shift-keying (OQPSK)-modulated signals, are also frequently encountered in digital communications. Therefore, there has been a recent surge of interest in deriving new algorithms that are able to properly handle noncircular signals as well [References-References]. These estimators extract additional information about the angular parameters from the unconjugated spatial covariance matrix that is non-zero for noncircular sources, in contrast to circular ones. From this perspective, we have been also able to propose a robust technique which is able to handle both temporally and spatially correlated sources in presence of noncircular signals [References]. By accounting for both signals’ noncircularity and temporal correlation, the proposed estimator was indeed shown to offer huge performance enhancements with respect to the main state-of-the-art techniques. Yet, all the aforementioned estimators [References-References] are applicable for the point-source model only. And, to the best of the authors’ knowledge, no contribution has dealt so far with the problem of angular parameters estimation (i.e., central DOAs and angular spreads) of multiple noncircular ID sources.

Motivated by these fact, we tackle in this paper for the very first time the problem of estimating the angular parameters of ID noncircular sources. We propose a new method that allows decoupling the estimation of each central DOA from its associated angular spread in the presence of noncircular sources. This method will be derived by going through three different stages resulting in two versions of the proposed estimator. The first one is a new 2-D search algorithm that extends ESB from circular to noncircular sources. And the second is a robust version that estimates the angular parameters by means of two successive one-dimensional (1-D) parameter searches. Towards this goal, we will use unstructured models for the conjugated and unconjugated noise-free covariance matrices that depend on the unknown angular spreads only. Most interestingly, such unstructured models are totally oblivious to the angular distributions of the sources and, therefore, their a priori knowledge is not required by the proposed method; a quite precious degree of freedom in practice. Even more, unlike all the existing methods, the proposed technique does not need to assume the same angular distribution across all the sources.

In order to properly assess the performance of the new estimator, we also conduct a complete theoretical study of its statistical properties (i.e., its bias and variance). Furthermore, we derive an explicit expression for the CRLB of the underlying estimation problem. This fundamental lower bound, which reflects the best achievable performance ever [References], will be used as an overall benchmark against which we gauge the accuracy of the new estimator. Computer simulations will show that the proposed estimator outperforms ESB and RGC especially at low SNR values and/or low DOA separations. The new CRLBs will also reveal that the noncircularity of the signals becomes more informative about the angular parameters when the sources have different angular distributions and when the angular spreads increase.

The rest of this paper is organized as follows. In Section II, we introduce the system model and some of the basic assumptions that will be adopted throughout the article. In section III, we derive the new algorithm and in section IV we show how the estimation of the central DOAs can be decoupled from that of the angular spreads. In section V, we derive the statistical bias and variance of the new estimator. In section VI, we derive an explicit expression for the CRB of the underlying estimation problem. Computer simulations are presented in Section VII and concluding remarks are drawn out in Section VIII.

We list beforehand some of the common notations adopted throughout this paper. Matrices and vectors are represented by bold upper- and lower-case characters, respectively. Vectors are by default in column orientation. Moreover, we consider the following standard notations:

[TABLE]

[TABLE]

II System Model

Consider an array consisting of $L$ identical sensors (i.e., with the same gain, phase, and sensitivity pattern) that is immersed in the far-filed of $K$ scattered ID sources with the same central frequency $\omega_{0}$. Assume that the root mean square (rms) delay spread is small compared to the inverse bandwidth of the transmitted signals so that the narrowband assumption remains valid in the presence of scattering [References-References]. Under these mild conditions, the signal received by the $l$th sensor, $l=1,2,\ldots,L$, can be modeled as follows [References-References]:

[TABLE]

in which $n$ stands for the $n$th snapshot. Moreover, $w_{l}(n)$ is an additive zero-mean circularly symmetric Gaussian-distributed noise. The noise components are assumed to be temporally and spatially white, i.e., uncorrelated between snapshots and receiving antenna branches, respectively. Furthermore, $s_{k}(\theta,\bar{\bm{\psi}}_{k},n)$ is the data-modulated angular distribution (with respect to $\theta$) of the signal received from the $k$th source; parameterized here by the vector $\bar{\bm{\psi}}_{k}=[\bar{\Theta}_{k},\bar{\sigma}_{k}]^{T}$.

For any planar configuration of the receiving antenna array, $a_{l}(\theta)$ can be written as:

[TABLE]

For mathematical convenience, we gather all the unknown central DOAs and angular spreads in the following parameter vectors:

[TABLE]

Our goal in the remainder of this paper is to jointly estimate the angular parameters, $\bar{\bm{\Theta}}$ and $\bar{\bm{\sigma}}$, of the $K$ noncircular sources given the set of received signals, $x_{l}(n)$, $l=1,2,\ldots,L$. To that end, we stack the received data over the $L$ sensors at each snapshot $n$ in a single vector:

[TABLE]

From (1), $\mathbf{x}(n)$ is explicitly given by:

[TABLE]

where

[TABLE]

are the array noise and response vectors, respectively. For ID sources, the components impinging from different scatterers are uncorrelated thereby yielding:

[TABLE]

where $\sigma_{s_{k}}^{2}$ is the average power of the $k$th source. Since the sources are also assumed to radiate noncircular signals, we adopt the definition of noncircularity in [References, References]. Moreover, we exploit the property of signals’ correlation in the real sense [References, property $3.1$] to prove from (8) that $p^{\prime}_{kk}(\theta,\theta^{\prime};\bar{\bm{\psi}}_{k})$ can be written as:

[TABLE]

Here, $0\leq\bar{\gamma}_{k}\leq 1$ and $\bar{\varphi}_{k}$ are the noncircularity rate and phase of the $k$th source, respectively. As emphasized in Section II, all existing works on angular parameters estimation of ID sources assume the sources to be circular. As such, none of them makes use of the unconjugated kernels in (9) since they are identically zero in this case. In this paper, however, we consider the case of noncircular sources with maximum noncircularity rate (i.e., $\bar{\gamma}_{k}=1$), known in the open literature as strictly second-order noncircular or rectilinear signals. Examples of such signals include unfiltered BPSK-, OQPSK-, PAM-, ASK-, AM- and MSK-modulated signals [References]. Their unconjugated angular auto-correlation kernels are obtained from (10) as:

[TABLE]

Now, since the sources’ signals are uncorrelated from the noise components, the conjugated and unconjugated covariance matrices of $\mathbf{x}(n)$ defined, respectively, as $\mathbf{R}_{\mathbf{x}\mathbf{x}}=\mathbb{E}\{\mathbf{x}(n)\mathbf{x}(n)^{H}\}$ and $\mathbf{R}^{\prime}_{\mathbf{x}\mathbf{x}}=\mathbb{E}\{\mathbf{x}(n)\mathbf{x}(n)^{T}\}$ are explicitly given by:

[TABLE]

where $\sigma_{w}^{2}$ is the unknown noise variance. Note here that the unconjugated covariance matrix of the circular noise vector is identically zero and, therefore, it vanishes in (13) contrarily to (II). By further assuming the ID sources to be mutually uncorrelated, it follows that:

[TABLE]

where $\delta_{kk^{\prime}}$ is the Kronecker delta function defined as $\delta_{kk^{\prime}}=1$ for $k=k^{\prime}$ and [math] otherwise. Now, plugging (8) and (11) in (14) and (15), respectively, leads to:

[TABLE]

Consequently, (II) and (13) simplify to:

[TABLE]

III Angular Parameters Estimation in Presence of Noncircular Signals

In order to exploit the additional information contained in the unconjugated covariance matrix of noncircular signals, we define the following extended received vector:

[TABLE]

whose extended covariance matrix is given by:

[TABLE]

On the one hand, using the explicit expressions of $\mathbf{R}_{\mathbf{x}\mathbf{x}}$ and $\mathbf{R}^{\prime}_{\mathbf{x}\mathbf{x}}$ established, respectively, in (18) and (19) and resorting to some algebraic manipulations, it can be shown that:

[TABLE]

where $\widetilde{\mathbf{a}}(\theta,\bar{\varphi}_{k})$ is the extended array response vector defined as:

[TABLE]

We also define the extended (normalized) covariance matrix of the noise-free signal pertaining to the $k$th source as:

[TABLE]

Hence, the extended covariance matrix in (III) is simply given by:

[TABLE]

Next, we consider the following eigendecomposition of the extended covariance matrix in (29):

[TABLE]

where $\widetilde{\mathbf{U}}_{s}$ and $\widetilde{\mathbf{U}}_{w}$ denote the eigenvector matrices associated to the signal and noise subspaces, respectively. Moreover, $\bm{\Sigma}$ is a diagonal matrix containing the eigenvalues of the overall extended noise-free covariance matrix involved in (29), i.e.:

[TABLE]

Traditional subspace-based methods which are all designed for circular ID sources rely on the fact that the columns of each $k$th noise-free covariance matrix, $\mathbf{R}^{(k)}_{ss}(\bar{\bm{\psi}}_{k})$, are orthogonal to those of the pseudo-noise subspace, i.e.:

[TABLE]

in which $r$ is the effective dimension of the pseudosignal subspace [References]. In principle, the same orthogonality property in (32) holds for noncircular ID sources:

[TABLE]

and can be used, as well, to estimate the associated angular parameters. However, similar to all subspace methods, the estimation performance is critically affected if the effective dimension $r$ is not appropriately selected. Besides, the optimal choice of $r$ depends on the value of the angular spread which is itself considered as an unknown parameter in our work. To sidestep this problem, we will rather capitalize on the inverse of the extended covariance matrix as recently done in [References]:

[TABLE]

To that end, let $\bm{\psi}$ and $\varphi$ be the two generic variables that run over all the possible values of $\bar{\bm{\psi}}_{k}$ and $\bar{\varphi}_{k}$, respectively. Then, right-multiplying (34) by $\widetilde{\mathbf{R}}^{(k)}_{ss}(\bm{\psi},\varphi)$ yields:

[TABLE]

At relatively high SNR levels, the signal eigenvalues in $\bm{\Sigma}$ are relatively large and, therefore, the diagonal elements of $\bm{\Sigma}^{-1}$ are almost equal to zero. Consequently, the first term in the right-hand side of (35) does not vary appreciably with $\bm{\psi}$ and $\varphi$. The second term in (35) is thus dominant. Owing to (33), however, it is identically zero when $\bm{\psi}=\bar{\bm{\psi}}_{k}$ and $\varphi=\bar{\varphi}_{k}$ (for $k=1,2,\ldots,K$). Therefore, at favorable SNR conditions, the quantity $\|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-1}~{}\!\widetilde{\mathbf{R}}^{(k)}_{ss}(\bm{\psi},\varphi)\|_{\textrm{Fro}}$ attains its minimum at $(\bar{\bm{\psi}}_{k},~{}\bar{\varphi}_{k})$ for each $k=1,2,\ldots,K$. Based on this observation, the angular parameters can be estimated jointly with the sources’ noncircularity phases by resolving the following $K$ optimization problems:

[TABLE]

where $\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$ is the sample-mean estimate of the actual extended covariance matrix, $\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$, i.e.:

[TABLE]

in which $N$ stands for the number of snapshots. Further, if the sources have the same scatterers’ angular distribution [i.e., $\widetilde{\mathbf{R}}^{(k)}_{ss}(\bm{\psi},\varphi)=\widetilde{\mathbf{R}}_{ss}(\bm{\psi},\varphi)$, $\forall k$], then all the angular parameters can be estimated jointly by finding the location of the $K$ smallest values of the common cost function:

[TABLE]

where

[TABLE]

Note here that the cost function in (39) to be minimized requires a three-dimensional (3-D) search over the central DOA, $\Theta$, the angular spread, $\sigma$, and the noncircularity phase, $\varphi$. In the following, we will try to reduce the complexity of the proposed method by reducing the dimensionality of the cost function in (39).

Actually, using (27) in (40), it can be shown that:

[TABLE]

where $\mathbf{R}_{ss}(\bm{\psi})$ and $\mathbf{R}^{\prime}_{ss}(\bm{\psi})$ are, respectively, the normalized conjugated and unconjugated noise-free auto-covariance matrices of the sources which are explicitly given by:

[TABLE]

Assuming small angular spreads, we prove in Appendix A that $\mathbf{R}_{ss}(\bar{\bm{\psi}}_{k})$ can be written for each $k$th source as:

[TABLE]

with $\bm{\Phi}(\bar{\Theta}_{k})~{}\triangleq~{}\textrm{diag}\{\mathbf{a}(\bar{\Theta}_{k})\}$ and $\mathbf{T}(\bar{\bm{\psi}}_{k})$ is a real-valued $(L\times L)$ symmetric matrix whose $(p,l)$th entry is given by:

[TABLE]

and $f^{\prime}_{p-1}(\theta)$ stands for the first derivative of $f_{p-1}(\theta)$ with respect to $\theta$.

In the same way, we also show that the normalized unconjugated noise-free covariance matrix, $\mathbf{R}^{\prime}_{ss}(\bar{\bm{\psi}}_{k})$, of noncircular ID sources can be be expressed as:

[TABLE]

where $\mathbf{T}^{\prime}(\bar{\bm{\psi}}_{k})$ is also a real-valued $(L\times L)$ symmetric matrix whose $(p,l)$th entry is given by:

[TABLE]

Injecting (47) and (49) back into (44) with the generic $\bm{\psi}$ and $\Theta$ being substituted for $\bar{\bm{\psi}}_{k}$ and $\bar{\Theta}_{k}$, respectively, and resorting to some straightforward manipulations, it can be shown that:

[TABLE]

in which $\widetilde{\bm{\Phi}}(\Theta,\varphi)~{}=~{}\textrm{diag}\{\widetilde{\mathbf{a}}(\Theta,\varphi)\}$ and

[TABLE]

For mathematical convenience, we also introduce the following notations:

[TABLE]

Then, plugging (51) back into (39), we prove after tedious manipulations (cf. Appendix B), that the angular parameters, $\{\bar{\bm{\psi}}_{k}\}_{k=1}^{K}$, can now be estimated by minimizing the following compressed cost function (i.e., that depends on $\bm{\psi}$ only):

[TABLE]

Note here that the first version of our proposed method defined by the cost function in (60) is applicable for a general class of angular distributions (symmetric distributions with small angular spreads) and any planar array configuration. However, it requires the a priori knowledge of the angular distributions to calculate the matrices $\mathbf{A}$ and $\mathbf{B}$ from (56) and (57), respectively. Furthermore, finding the $K$ minima of (60) with respect to \bm{\psi}=\big{[}\Theta,~{}\sigma\big{]}^{T} still requires a two-dimensional (2-D) search over $\Theta$ and $\sigma$ and needs the angular distribution to be identical for all the sources to estimate jointly the angular parameters. In the following, we will build upon some properties of the matrices $\mathbf{T}(\bm{\psi})$ and ${\mathbf{T}^{\prime}}(\bm{\psi})$ in order to decouple the estimation of the central DOAs from that of the angular spreads [References]. These properties are valid for any symmetric source’s angular distribution with small angular spreads. Hence, the estimator can be implemented by two successive one-dimensional (1-D) parameter searches, thereby resulting in tremendous computational savings. Moreover, we will exploit these properties to establish unstructured models for $\mathbf{T}(\bm{\psi})$ and ${\mathbf{T}^{\prime}}(\bm{\psi})$ that are totally oblivious to the symmetric sources’ angular distributions. Therefore, we will obtain a new version of the proposed estimator that does not require the a priori knowledge of the sources’ angular distributions.

IV Robust Version of the Proposed Estimator

To begin with, for any array configuration, recall that $\mathbf{T}(\bm{\psi})$ is a real-valued symmetric matrix whose expression is given by (III). Moreover, we prove in the following that if $f^{\prime}_{p-1}(\bar{\Theta}_{k})$ is expressed as follows111(61) means that the antenna array must be an equally-spaced linear array.:

[TABLE]

where $g(\bar{\Theta}_{k})$ is a transformation of the central DOA $\bar{\Theta}_{k}$, then $\mathbf{T}(\bar{\bm{\psi}}_{k})$ is a symmetric Toeplitz matrix. In fact, injecting (61) in (III), we show that $[\mathbf{T}]_{pl}(\bar{\bm{\psi}}_{k})$ can be written as follows:

[TABLE]

From (62), we can simply verify that:

[TABLE]

Consequently, $\mathbf{T}(\bar{\bm{\psi}}_{k})$ is a symmetric Toeplitz matrix and, therefore, it can be fully constructed from its first column vector denoted here as $\mathbf{t}_{1}$, i.e.:

[TABLE]

Moreover, for any symmetric angular distribution, we prove in Appendix C that if its angular spread verifies the following condition:

[TABLE]

then the elements, $\{\mathbf{t}_{1}(l)\}_{l=1}^{L}$, of the vector, $\mathbf{t}_{1}$, satisfy the following property:

[TABLE]

(65) is a nonrestrictive condition for propagation environments characterized by small angular spreads, e.g., macro-cell environments [References-References]. Actually, (66) can be rewritten in the more succinct form222Note here that the notation $\mathbf{v}_{1}\leq\mathbf{v}_{2}$ for any tow $N-$dimensional vectors $\mathbf{x}=[x_{1},x_{2},\ldots,x_{N}]^{T}$ and $\mathbf{y}=[y_{1},y_{2},\ldots,y_{N}]^{T}$ means that $\mathbf{x}_{n}\leq y_{n}$ for $n=1,2,\ldots,N$.:

[TABLE]

where $\mathbf{J}_{n}$ is from now on a ($n\times n$) matrix given by:

[TABLE]

and $\mathbf{e}_{n}$ is a $n-$dimensional vector given by $\mathbf{e}_{n}=[1,0,\ldots,0]^{T}$.

For any array configuration, recall also that ${\mathbf{T}^{\prime}}(\bm{\psi})$ is a real-valued symmetric matrix whose expression is given by (III). Moreover, if $f^{\prime}_{p-1}(\bar{\Theta}_{k})$ satisfies (61), we show that $[\mathbf{T}^{\prime}]_{pl}(\bar{\bm{\psi}}_{k})$ can be written as:

[TABLE]

From (75), we can see that ${\mathbf{T}^{\prime}}(\bm{\psi})$ is a Hankel matrix. Therefore, it can be constructed from its first and last column vectors denoted, respectively, as $\mathbf{t}^{\prime}_{1}$ and $\mathbf{t}^{\prime}_{L}$ as follows:

[TABLE]

Moreover, for any symmetric source’s angular distribution, we also prove in Appendix C that if \sigma<1/\big{(}2\sqrt{2}\pi(L-1)g(\bar{\Theta}_{k})\big{)}, then the elements of $\mathbf{t}^{\prime}_{1}$ and $\mathbf{t}^{\prime}_{L}$ satisfy the following properties:

[TABLE]

Furthermore, we verify that the first column vector of ${\mathbf{T}^{\prime}}(\bm{\psi})$ is identical to the first column vector of $\mathbf{T}(\bm{\psi})$, i.e., we have the following relation:

[TABLE]

In order to exploit the interesting properties stated above in (67), (77) and (78), we consider an auxiliary vector $\mathbf{z}=[z_{1},\ldots,z_{L-1},z_{L},\ldots,z_{2L-2}]^{T}$ whose elements are all in $[0,1]$ and sorted in decreasing order:

[TABLE]

or equivalently:

[TABLE]

Then, we construct the following two auxiliary matrices:

[TABLE]

which also verify the constraints in (67) and (77), respectively. Therefore, bearing in mind the expressions of the matrices $\mathbf{A}$ and $\mathbf{B}$ in (56) and (57), respectively, it follows that instead of minimizing the 2-D criterion in (60), one can start by solving the following 1-D constrained optimization problem in order to find the central DOAs:

[TABLE]

where

[TABLE]

The optimization task in (83) can be solved efficiently via the well-known sequential quadratic programming (SQP) algorithm which is a rapidly converging descent method for nonlinearly-constrained optimization problems [References]. Interestingly enough, the estimator in (83) is also totally oblivious to the sources’ angular distributions provided that the latter be symmetric. In fact, the auxiliary matrices $\mathbf{Z}$ and $\mathbf{Z}^{\prime}$ involved in (80) were built for any symmetric angular distribution upon some general properties shared by $\mathbf{T}$ and $\mathbf{T}^{\prime}$, respectively, and not their true expressions as required in (60). Moreover, this estimator is applicable in the more challenging scenario where the sources have different angular distributions. These are actually quite precious degrees of freedom in practice since the angular distribution may vary from one environment to another and/or from source to source in real-world scenarios. After acquiring the central DOAs, \widehat{\bar{\mathbf{\Theta}}}\!~{}=~{}\!\big{[}\widehat{\bar{\Theta}}_{1},\widehat{\bar{\Theta}}_{2},\ldots,\widehat{\bar{\Theta}}_{K}\big{]}, as in (83), the angular spread pertaining to each $k$th source is estimated as follows:

[TABLE]

where $f_{c}(.)$ is the compressed cost function already established in (60). Our estimator actually reduces to (83) and (85), that is after going through three different derivation stages, from (39) to (60), ultimately leading to our final robust solution.

V Statistical Properties

In order to assess the theoretical performance limits of the proposed estimator, we will express its mean square error (MSE) analytically, based on the minimization of the original cost function, f\big{(}\bm{\psi},\varphi~{}\!\big{|}\!~{}\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}\big{)}, given in (39) instead of the compressed one in (60) due to the presence of the nonlinear modulus operator in it. We will also use $\bm{\alpha}_{k}=[\bar{\bm{\psi}}^{T}_{k},\bar{\varphi}_{k}]^{T}=[\bar{\Theta}_{k},\bar{\sigma}_{k},\bar{\varphi}_{k}]^{T}$ to denote the entire unknown parameter vector pertaining to each $k$th noncircular ID source. We further use:

[TABLE]

[TABLE]

to denote, respectively, the bias vector and covariance matrix of the estimate $\widehat{\bm{\alpha}}_{k}$. To begin with, it is easy to show that the mean square error (MSE), defined as \textrm{MSE}(\widehat{\bm{\alpha}}_{k})\triangleq\mathbb{E}\big{\{}(\widehat{\bm{\alpha}}_{k}-\bm{\alpha}_{k})(\widehat{\bm{\alpha}}_{k}-\bm{\alpha}_{k})^{T}\big{\}}, is given by:

[TABLE]

Moreover, similar to [References], let $\breve{\bm{\alpha}}_{k}$ denote the asymptotic estimate (obtained when the number of snapshots $N\rightarrow\infty$), and define $\Delta\bm{\alpha}_{k}~{}\triangleq~{}\breve{\bm{\alpha}}_{k}~{}-~{}\bm{\alpha}_{k}$ and $\Delta\breve{\bm{\alpha}}_{k}~{}\triangleq~{}\widehat{\bm{\alpha}}_{k}-\breve{\bm{\alpha}}_{k}$. Then, it immediately follows from (86) that $\mathbf{b}(\widehat{\bm{\alpha}}_{k})$ decomposes as the sum of the asymptotic bias and the residual bias stemming from the finite-sample effects:

[TABLE]

Furthermore, using some relatively straightforward algebraic manipulations, it can be shown that:

[TABLE]

Plugging (89) and (90) back into (88), it follows that:

[TABLE]

In order to establish an analytical expression for $\textrm{MSE}(\hat{\bm{\alpha}}_{k})$, we will derive hereafter the four expectations involved in (91) separately. To do so, we use $\bm{\alpha}$ and $\mathbf{R}$ as generic variables for $\bm{\alpha}_{k}$ and $\widehat{\mathbf{R}}^{-2}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$, respectively. We also denote the gradient vector and Hessian matrix of the scalar-valued objective function, $f(\bm{\alpha},\mathbf{R})$, as follows:

[TABLE]

From (39), the $i^{th}$ element of the vector $\mathbf{f}(\bm{\alpha}|\mathbf{R})$ is given by:

[TABLE]

where

[TABLE]

Furthermore, the entries of the Hessian matrix, $\mathbf{F}(\bm{\alpha}|\mathbf{R})$, are obtained as follows:

[TABLE]

where

[TABLE]

V-A Derivation of*  $\mathbb{E}\{\Delta\bm{\alpha}_{k}\}$  and  $\mathbb{E}\{\Delta\bm{\alpha}_{k}\Delta\bm{\alpha}_{k}^{T}\}$:*

First, using the properties of the complex Wishart distribution [References, p. $273$], it can be shown that the asymptotic sample-mean estimate of the extended covariance matrix, $\breve{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}=\lim_{N\to\infty}\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$, is a consistent estimate of the Hermitian matrix $\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$. This means that as $N\rightarrow\infty$, we have:

[TABLE]

To derive the asymptotic bias, $\Delta\bm{\alpha}_{k}$, we use as in [References, References] the first-order Taylor series expansion of $\mathbf{f}(\bm{\alpha},\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$ around the actual parameter vector $\bm{\alpha}_{k}$:

[TABLE]

By noticing that the asymptotic estimate, $\breve{\bm{\alpha}}_{k}$, also minimizes $f(\bm{\alpha}|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$, it follows that $\mathbf{f}(\breve{\bm{\alpha}}_{k}|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})=\mathbf{0}$. Therefore, by evaluating (96) at $\bm{\alpha}=\breve{\bm{\alpha}}_{k}$, it follows that:

[TABLE]

from which $\Delta\bm{\alpha}_{k}\triangleq(\breve{\bm{\alpha}}_{k}-\bm{\alpha}_{k})$ is obtained as:

[TABLE]

Consequently, the approximate expression for the asymptotic bias, $\Delta\bm{\alpha}_{k}$, and \mathbb{E}\big{\{}\Delta\bm{\alpha}_{k}\Delta\bm{\alpha}^{T}_{k}\big{\}} are obtained as follows:

[TABLE]

and

[TABLE]

V-B Derivation of*  $\mathbb{E}\{\Delta\breve{\bm{\alpha}}_{k}\}$  and  $\mathbb{E}\{\Delta\breve{\bm{\alpha}}_{k}\Delta\breve{\bm{\alpha}}_{k}^{T}\}$:*

After tedious algebraic manipulations, we also show in Appendix D that $\Delta\breve{\bm{\alpha}}_{k}$ is expressed as follows:

[TABLE]

where $\mathbf{v}(\breve{\bm{\alpha}}_{k}|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2},\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$ is a 3-dimensional vector whose $i^{th}$ element is explicitly given by

[TABLE]

Recall here that $\widetilde{\mathbf{R}}^{[i]}_{ss}$ was already defined in (93) and we further define $\Delta\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}$ as follows:

[TABLE]

Then, by exploiting the fact that $\partial\textrm{tr}\{\mathbf{A}\mathbf{B}\}/\partial\mathbf{A}=\mathbf{B}^{T}$ for any two matrices $\mathbf{A}$ and $\mathbf{B}$, it follows that:

[TABLE]

where the last equality follows from the identity $\textrm{tr}\{\mathbf{A}\mathbf{B}\}=\textrm{vec}^{T}\{\mathbf{A}^{T}\}\textrm{vec}\{\mathbf{B}\}$. Consequently, the vector \mathbf{v}\big{(}\breve{\bm{\alpha}}_{k}|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2},\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}\big{)} is expressed as follows:

[TABLE]

where the matrix $\mathbf{G}_{ss}$ is given by

[TABLE]

Plugging (103) back into (99), one obtains:

[TABLE]

whose expectation yields the required residual bias as follows:

[TABLE]

Furthermore, in presence of noncircular signals, it can be shown that $\mathbb{E}\{\Delta\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}\}$ is accurately approximated by333See [References] and [References] for more details about the proof in the case of circular sources that we generalize here to the noncircular case using the appropriate extended covariance matrices.:

[TABLE]

where $\widetilde{\mathbf{e}}_{n}$ is an eigenvector associated to the $n^{th}$ eigenvalue, $\lambda_{n}$, of the extended covariance matrix $\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$. From (105), it also immediately follows that:

[TABLE]

where

[TABLE]

The entries of $\mathbf{H}$ are also evaluated using the following accurate approximation:

[TABLE]

in which the weighting coefficients, $\omega_{nn^{\prime}}$, are simply given by:

[TABLE]

VI New CRLB for Noncircular Gaussian Distributed Signals Generated from ID Sources

In this section, we assume that the transmitted signals $\{\bm{s}(t)\}_{t=1,2,\ldots,N}$ are zero-mean Gaussian distributed and generated from noncircular ID sources. We also assume that the noncircularity rate of the signals is $0\leq\gamma\leq 1$. Now recall from (24) that the extended covariance matrix of the received signals is given by:

[TABLE]

Moreover, using (45) and (46) in (18) and (19), respectively, it follows that:

[TABLE]

Then, using (47) and (49) in (114) and (115), respectively, leads to:

[TABLE]

Recall also that the explicit expressions of $\bm{\Phi}(\bar{\Theta}_{k})$, $\mathbf{T}(\bar{\bm{\psi}}_{k})$ and ${\mathbf{T}^{\prime}}(\bar{\bm{\psi}}_{k})$ were already given in Section III. Our goal in this section is to find the CRLB of the unknown parameters of interest (i.e., namely the angular parameters) which are gathered in the following vector:

[TABLE]

The unknown nuisance parameters which are the noise variance, $\sigma_{w}$, the sources’ powers, $\bm{\beta}~{}\triangleq~{}\left[\sigma_{s_{1}}^{2},\ldots,\sigma_{s_{K}}^{2}\right]^{T}$, and their noncircularity phases, $\bm{\varphi}~{}\triangleq~{}\left[\varphi_{1},\ldots,\varphi_{K}\right]^{T}$, are also gathered in the vector:

[TABLE]

We will also group all the parameters in (119) an (118) in a single vector:

[TABLE]

The CRLB of the entire unknown parameter vector, $\bm{\upsilon}$, is defined as follows [References]:

[TABLE]

where $\mathbf{I}(\bm{\upsilon})$ is the so-called Fisher information matrix (FIM). Since the extended snapshot vectors, $\{\widetilde{\mathbf{x}}(t)\}_{t=1}^{N}$, defined in (21) are mutually independent, then according to [References] the $(i,j)$th entry of the FIM associated to the underlying estimation problem is given by:

[TABLE]

where $\upsilon_{i}$ is the $i$th element of the whole parameter vector given in (120). Using (116) and (117), we show in Appendix E that the CRLB for the angular parameters alone in presence of uncorrelated ID noncircular sources is explicitly given by:

[TABLE]

where the expressions of $\mathbf{I}_{\bm{\eta},\bm{\eta}}$, $\mathbf{I}_{\bm{\xi},\bm{\eta}}$, and $\mathbf{I}_{\bm{\xi},\bm{\xi}}$ are provided in Appendix E.

VII Simulation Results

In this section, we assess the performance of the newly proposed method and gauge it against the most recent state-of-the-art techniques that are geared toward multiple ID sources, namely ESB [References] and RGC [References]. Although the latter were derived specifically for ID circular sources, they can be applied to the noncircular case as well after completely ignoring the non-zero unconjugated covariance matrix. All the methods will be also gauged against the CRLB. In all simulations, we consider complex Gaussian transmitted signals and a uniform linear array of $6$ sensors separated by half a wavelength.

VII-A Assessment of the new estimator

In this subsection, the root mean-square error (RMSE) of each estimator is computed empirically by means of $2000$ Monte-Carlo runs. We first consider in Fig. 1 two uncorrelated ID noncircular sources with the same noncircularity rate ($\gamma_{1}=\gamma_{2}=1$) and noncircularity phases $\varphi_{1}=\frac{\pi}{3}$ and $\varphi_{2}=\frac{\pi}{4}$. Both ID sources have a Gaussian angular distribution (i.e., GID) and are located at central DOAs $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$ and $\bar{\Theta}_{2}=30{\textrm{\textdegree}}$ with respective angular spreads $\bar{\sigma}_{1}=1.5{\textrm{\textdegree}}$ and $\bar{\sigma}_{2}=3{\textrm{\textdegree}}$. The SNR is fixed to $5$ dB while the number of snapshots used to estimate the sample covariance matrix is increased from $100$ to $1000$ in steps of $100$. Figs. 1(a) and 1(b) depict the empirical RMSEs of all tested methods.

Clearly, our estimator is statistically more efficient and outperforms ESB and RGC both in terms of central DOAs and angular spreads estimation accuracy. Moreover, the performance improvements of the proposed method over ESB and RGC hold almost the same irrespectively of $N$. Therefore, we will hereafter fix $N=1000$.

Figs. 2(a) and 2(b) depict the empirical RMSEs of all tested methods versus the SNR. The analytical RMSE of the new estimator established in Section V is also plotted there. These figures show a very good agreement between the empirical and analytical RMSEs of the proposed estimator, thereby corroborating our analytical performance analysis of Section V. It also suggests that the proposed estimator outperforms ESB and RGC, both in terms of central DOAs and angular spreads estimation capabilities, especially under the adverse conditions of low SNR levels.

In Fig. 3, we consider two uncorrelated ID noncircular sources with different angular distributions. More specifically, the first source is uniformly distributed (UID) with central DOA, $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$, and angular spread $\bar{\sigma}_{1}=1.5{\textrm{\textdegree}}$ while the second is GID distributed with central DOA, $\bar{\Theta}_{2}=30{\textrm{\textdegree}}$, and angular spread $\bar{\sigma}_{2}=3{\textrm{\textdegree}}$.

To apply ESB in this setup, however, we assume that both sources are GID. In fact, in contrast to the proposed method and RGC, ESB was specifically derived in the case where all the sources have the same angular distribution. By comparing Figs. 2 and 3 (i.e., sources truly having the same distribution), we observe that ESB suffers from severe performance degradation. It even becomes less accurate than RGC at low SNR levels, that is in stark contrast to what was earlier reported in Fig. 2. The proposed estimator, however, keeps its superiority in terms of estimation accuracy thereby making it more attractive in practice where the sources are more likely to have different angular distributions.

Next, we examine the impact of the sources’ separation on the performance of the three estimators. To that end, we reconsider the case of noncircular ID sources with the same angular distribution (GID). The first source is kept fixed at $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$ with angular spread, $\bar{\sigma}_{1}=1.5{\textrm{\textdegree}}$, while the second (with $\bar{\sigma}_{2}=3^{\circ}$) is shifted from $18{\textrm{\textdegree}}$ to $30{\textrm{\textdegree}}$ with $2^{\circ}$. The results are plotted in Fig. 4 at $5$ dB SNR and suggest that all estimators expectedly improve their accuracy as the DOA separation increases. Yet, the proposed approach significantly outperforms ESB and RGC for small DOA separations, a more challenging scenario in practice.

Finally, we consider in Fig. 5 an even more challenging scenario where two uncorrelated GID noncircular sources with the same noncircularity rate ($\gamma_{1}=\gamma_{2}=1$) and noncircularity phases $\varphi_{1}=\frac{\pi}{3}$ and $\varphi_{2}=\frac{\pi}{4}$ are located at central DOAs $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$ and $\bar{\Theta}_{2}=15{\textrm{\textdegree}}$ with respective angular spreads $\bar{\sigma}_{1}=2{\textrm{\textdegree}}$ and $\bar{\sigma}_{2}=4{\textrm{\textdegree}}$. The number of snapshots is fixed to $N=100$. Figs. 5(a) and 5(b) show that the performance of the three methods is satisfactory, especially at high SNR values. However, our new estimator still outperforms the two other methods both in terms of central DOAs and angular spreads estimation performance.

VII-B *Assessment of the new CRLBs: *

In this subsection, we illustrate the newly derived CRLBs (i.e., $\textrm{CRLB}^{\textrm{noncir}}$) in different scenarios. We first consider two equipowered ID sources with identical noncircularity rate, $\gamma=1$, and noncircularity phases $\varphi_{1}=\pi/3$ and $\varphi_{2}=\pi/4$. The sources are located at central DOAs $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$ and $\bar{\Theta}_{2}=30{\textrm{\textdegree}}$ with respective angular spreads $\bar{\sigma}_{1}=3{\textrm{\textdegree}}$ and $\bar{\sigma}_{2}=5{\textrm{\textdegree}}$. Figs. 6(a) and 6(b) show both $\log(\textrm{CRLB}^{\textrm{noncir}})$ and $\log(\textrm{CRLB}^{\textrm{cir}})$ of $\bar{\Theta}_{1}$ and $\bar{\sigma}_{1}$, respectively, when the sources have: $i$) the same Gaussian angular distribution, and $ii$) different angular distributions (the first source is UID and the second source is GID).

We see from Fig. 6 that the CRLBs for noncircular ID sources are lower than their counterparts derived assuming ID circular sources, especially at low SNR values. This illustrates the performance gain that is achieved by exploiting the non-circularity feature of the sources in the estimation process. Moreover, $\textrm{CRLB}^{\textrm{cir}}$ converges faster to $\textrm{CRLB}^{\textrm{noncir}}$, at high SNR, when the sources have the same angular distribution (GID-GID in our case). Therefore, at high SNRs, the noncircularity of the signals is more informative about the angular parameters when the sources have different distributions. Next, we examine the impact of the angular spread on the estimation of the angular parameters, by fixing $\bar{\sigma}_{2}$ and varying $\bar{\sigma}_{1}$. Fig. 7 depicts $\log(\textrm{CRLB}^{\textrm{noncir}})$ and $\log(\textrm{CRLB}^{\textrm{cir}})$ as a function of the SNR for three different values of $\bar{\sigma}_{1}$. Moreover, we consider in Fig. 7(a) the case of point (or non-distributed) sources which corresponds to $\bar{\sigma}_{1}=\bar{\sigma}_{2}=0{\textrm{\textdegree}}$.

As intuitively expected, $\textrm{CRLB}^{\textrm{noncir}}$ and $\textrm{CRLB}^{\textrm{cir}}$ increase with the angular spread and so does the difference between them. This reveals that as the angular spread increases, there is more room for the noncircularity of the signals to improve the estimation performance. In fact, the signals become more dispersed and thus the unconjugated covariance matrix becomes more informative about the angular parameters. In Figs. 8 and 9, we study the effect of the signals’ noncircularity parameters on CRLB${}^{\textrm{noncir}}$ under different sources’ separations, $\Delta\Theta$, in terms of central DOAs. The first source is UID and fixed at $\bar{\Theta}_{1}=10{\textrm{\textdegree}}$ whereas the second source is GID and its central DOA, $\bar{\Theta}_{2}$, is varied from $18{\textrm{\textdegree}}$ to $30{\textrm{\textdegree}}$.

We observe from Figs. 8-(a) and 8-(b) that $\textrm{CRLB}^{\textrm{noncir}}$ of the two angular parameters decrease as the noncircularity steps rate increases. Moreover, the gap between $\textrm{CRLB}^{\textrm{noncir}}$ and $\textrm{CRLB}^{\textrm{cir}}$ increases as the DOA separation $\Delta\Theta$ decreases. In fact, the ratio between the two CRLBs tends to zero at low DOA separations (for $\Delta\Theta=8^{\circ}$). More specifically, at low DOA separations, $\textrm{CRLB}^{\textrm{noncir}}$ becomes very small compared to $\textrm{CRLB}^{\textrm{cir}}$, meaning that huge performance gains can be achieved in this challenging scenario by exploiting the additional information carried by the unconjugated covariance matrix. Fig. 9 also reveals that $\textrm{CRLB}^{\textrm{noncir}}$ is more sensitive to the noncircularity phase separation at small DOA separations.

VIII Conclusion

In this paper, we developed a new method for the estimation of the angular parameters in the presence of noncircular ID sources. The new estimator decouples the estimation of the central DOAs from that of the angular spreads by means of two consecutive 1-D searches, thereby resulting in tremendous computational savings as compared to the brute-force 2D grid search solution. It is also oblivious to the sources’ angular distribution or any mismatch thereof. This estimator is particulary interessant for symmetric sources’ distributions with small angular spreads.

The proposed estimator outperforms most recent state-of-the-art techniques, especially for small DOA separations and/or low SNR levels. Its performance was also assessed analytically and the obtained results were corroborated by Monte-Carlo simulations. In order to benchmark the new estimator, we also derived for the first time an explicit expression for the stochastic CRLBs of the underlying estimation problem. The analysis of the new CRLB unambiguously shows that the noncircularity of the signals brings valuable additional information about the angular parameters especially when the sources have different angular distributions and/or when the angular spreads increase. Besides, the noncircular CRLBs decrease as the noncircularity rate increases. And, they are much smaller than the circular CRLBs at small DOA separations. In which case they also become more sensitive to the noncircularity phase separation.

Appendix A—Proof of (III)

From (45), the $(p,l)$th element of $\mathbf{R}_{ss}(\bar{\bm{\psi}}_{k})$ has the following expression:

[TABLE]

where $a_{p}(\theta)=e^{j2\pi f_{p-1}(\theta)}$. Otherwise, we denote by $\widetilde{\theta}$ the deviation of the direction $\theta$ from the central DOA $\bar{\Theta}_{k}$ as follows:

[TABLE]

For small angular spreads, $\widetilde{\theta}$ tends to zero. We can therefore use the following approximation:

[TABLE]

where $f^{\prime}_{p-1}(\theta)$ stands for the first derivative of $f_{p-1}(\theta)$ with respect to $\theta$. Hence, we obtain the following expression for \big{[}\mathbf{R}_{ss}\big{]}_{pl}(\bar{\bm{\psi}}_{k}):

[TABLE]

\big{[}\mathbf{R}_{ss}\big{]}_{pl}(\bar{\bm{\psi}}_{k}) can be written equivalently as follows:

[TABLE]

where \big{[}\mathbf{T}\big{]}_{pl}(\bar{\bm{\psi}}_{k}) is given by:

[TABLE]

From (129), the complex conjugate of $[\mathbf{T}]_{pl}(\bar{\bm{\psi}}_{k})$ is given by:

[TABLE]

Moreover, by assuming that the angular distribution is symmetric with respect to the central DOA $\bar{\Theta}_{k}$, we have the following relation:

[TABLE]

From (131), it follows that:

[TABLE]

From (130) and (132), we have that:

[TABLE]

This proves that $\mathbf{T}(\bar{\bm{\psi}}_{k})$ is a real-valued $(L\times L)$ symmetric matrix. Otherwise, $\exp\{j2\pi\left(f^{\prime}_{p-1}(\bar{\Theta}_{k})-f^{\prime}_{l-1}(\bar{\Theta}_{k})\right)\widetilde{\theta}\}$ can be written as:

[TABLE]

Therefore, $[\mathbf{T}]_{pl}(\bar{\bm{\psi}}_{k})$ can be expressed as follows:

[TABLE]

Since $\mathbf{T}(\bar{\bm{\psi}}_{k})$ is a real-valued matrix, then we can deduce that:

[TABLE]

Consequently, $[\mathbf{T}]_{pl}(\bar{\bm{\psi}}_{k})$ can be reduced to:

[TABLE]

thereby leading to the result given in (III). Moreover, since, we have:

[TABLE]

then we can conclude that $\forall p,l$, we have:

[TABLE]

Appendix B—Proof of (154)

Substituting (51) in (39) and using the identity $\textrm{tr}\{\mathbf{A}\mathbf{B}\mathbf{C}\}=\textrm{tr}\{\mathbf{C}\mathbf{A}\mathbf{B}\}$ for any square matrices, $\mathbf{A},\mathbf{B}$ and $\mathbf{C}$, along with the fact that $\widetilde{\bm{\Phi}}(\bar{\Theta},\varphi)^{H}\widetilde{\bm{\Phi}}(\bar{\Theta},\varphi)=\mathbf{I}_{2L}$, we obtain:

[TABLE]

On the other hand, by recalling the expression of the extended array response vector, $\widetilde{\mathbf{a}}(\bar{\Theta},\varphi)$, in (27), it follows that $\widetilde{\bm{\Phi}}(\bar{\Theta},\varphi)=\textrm{diag}\{\widetilde{\mathbf{a}}(\bar{\Theta},\varphi)\}$ is given by:

[TABLE]

Furthermore, by recalling the expression of $\widetilde{\mathbf{T}}(\bm{\psi})$ in (55), it is easy to show that:

[TABLE]

where

[TABLE]

Similar to (30), the estimated extended covariance matrix, $\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}$, is eigendecomposed as follows:

[TABLE]

from which it can be shown that:

[TABLE]

Moreover, as shown in [References], $\widehat{\widetilde{\mathbf{U}}}_{s}$ and $\widehat{\widetilde{\mathbf{U}}}_{w}$ can be partitioned as follows:

[TABLE]

where $\mathbf{D}_{s}$ and $\mathbf{D}_{w}$ are some diagonal matrices whose complex diagonal entries are of unit modulus. Injecting (147) and (148) back into (146), we show after some algebraic manipulations that $\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}$ has the following block diagonal structure:

[TABLE]

where

[TABLE]

Substituting (141), (144), and (151) back into (137), and resorting to some algebraic manipulations yields the following result:

[TABLE]

in which the complex numbers $z_{1}(\bm{\psi})$ and $z_{2}(\bm{\psi})$ are explicitly given by:

[TABLE]

In order to reduce the dimensionality of the optimization problem at hand, we begin by minimizing the underlying cost function with respect to the unknown noncircularity phase $\varphi$. To that end, we use $z_{2}(\bm{\psi})=|z_{2}(\bm{\psi})|{\exp\{j\angle z_{2}(\bm{\psi})}\}$ and rewrite (154) as follows:

[TABLE]

From (158), it is clear (for a fixed $\bm{\psi}$) that the function $f(.)$ attains its minimum (with respect to $\varphi$) at the point:

[TABLE]

Substituting (159) back into (158) and recalling (155) and (156), we obtain the following cost function that depends on $\bm{\psi}$ only:

[TABLE]

which is equivalent to the result given in (60).

Appendix C—Proof of (66) and (77)

$\mathbf{T}(\bm{\psi})$ is a symmetric Toeplitz matrix constructed from its first column vector $\mathbf{t}_{1}$ as follows:

[TABLE]

where the $l$th element of the vector $\mathbf{t}_{1}$ is given from (62) by:

[TABLE]

For small angular spreads, we use a second-order Taylor-series development of $\cos\left(2\pi(l-1)g(\bar{\Theta}_{k})\!(\theta\!-\!\bar{\Theta}_{k})\right)$ to obtain the following equality:

[TABLE]

Then, $\mathbf{t}_{1}(l)$ can be approximated as follows:

[TABLE]

From (163), we clearly see that the elements, $\{\mathbf{t}_{1}(l)\}_{l=1}^{L}$, of the vector, $\mathbf{t}_{1}$, satisfy the following property:

[TABLE]

Moreover, if \sigma<1/\big{(}\sqrt{2}\pi(l-1)g(\bar{\Theta}_{k})\big{)}~{}\forall~{}\!l=1,\ldots,L or equivalently \sigma<1/\big{(}\sqrt{2}\pi(L-1)g(\bar{\Theta}_{k})\big{)}, then $\mathbf{t}_{1}(L)\geq 0$.

In the same way, ${\mathbf{T}^{\prime}}(\bm{\psi})$ is a Hankel matrix defined from its first and last column vectors $\mathbf{t}^{\prime}_{1}$ and $\mathbf{t}^{\prime}_{L}$ as follows:

[TABLE]

where the $l$th elements of the vectors $\mathbf{t}^{\prime}_{1}$ and $\mathbf{t}^{\prime}_{L}$ are given, respectively, by:

[TABLE]

For small angular spreads, $\mathbf{t}^{\prime}_{1}(l)$ and $\mathbf{t}^{\prime}_{L}(l)$ can be approximated as follows:

[TABLE]

From (168) and (169), we see clearly that the elements of $\mathbf{t}^{\prime}_{1}$ and $\mathbf{t}^{\prime}_{L}$ satisfy the following properties:

[TABLE]

Moreover, if \sigma<1/\big{(}2\sqrt{2}\pi(L-1)g(\bar{\Theta}_{k})\big{)}, then $\mathbf{t^{\prime}}_{L}(L)\geq 0$.

Appendix D—Proof of (99)

For any small vector and matrix perturbations, $\delta\mathbf{x}$ and $\delta\mathbf{X}$, and scalar-valued function $g(\mathbf{x},\mathbf{X})$, we have the following Taylor series expansion [References, References] around $\mathbf{x}_{0}$ and $\mathbf{X}_{0}$:

[TABLE]

The result in (Appendix D—Proof of (99)) is applied with $\mathbf{x}=\bm{\alpha}$ and $\mathbf{X}=\mathbf{R}$ to the functions:

[TABLE]

in order to obtain their Taylor series expansions around the point $(\mathbf{x}_{0},\mathbf{X}_{0})=(\breve{\bm{\alpha}}_{k},\breve{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$. In this way, for $k=1,2,\ldots K$, the underlying perturbations are given by:

[TABLE]

By doing so, we obtain for $i=1,2,3$:

[TABLE]

where

[TABLE]

Moreover, notice from (92) and (94) that $f_{\alpha_{i}}(\bm{\alpha}|\mathbf{R})$ and $\left[\partial f_{\alpha_{i}}(\bm{\alpha}|\mathbf{R})/\partial\bm{\alpha}\right]^{T}$ are, respectively, the $i^{th}$ element of the gradient vector, $\mathbf{f}(\bm{\alpha}|\mathbf{R})$, and the $i^{th}$ row of the Hessian matrix $\mathbf{F}(\bm{\alpha}|\mathbf{R})$. Therefore, by further defining the vector $\mathbf{v}=[v_{1},v_{2},v_{3}]^{T}$, the results of (Appendix D—Proof of (99)) for $i=1,2,3$ are rewritten in the following matrix/vector form:

[TABLE]

Evaluating the expansion in (178) at (\bm{\alpha},\mathbf{R})=\big{(}\widehat{\bm{\alpha}}_{k}, \widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}\big{)} and using $\Delta\breve{\bm{\alpha}}_{k}\triangleq\widehat{\bm{\alpha}}_{k}-\breve{\bm{\alpha}}_{k}$ leads to:

[TABLE]

The finite-sample and asymptotic estimates, $\widehat{\bm{\alpha}}_{k}$ and $\breve{\bm{\alpha}}_{k}$, are obtained by minimizing $f(\bm{\alpha}|\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$ and $f(\bm{\alpha}|\breve{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$, respectively. Therefore, the gradient of the latter objective function is identically zero at $\widehat{\bm{\alpha}}_{k}$ and $\breve{\bm{\alpha}}_{k}$, i.e.: ​​​​​​

[TABLE]

Exploiting (181) back into (Appendix D—Proof of (99)) and resolving for $\Delta\breve{\bm{\alpha}}_{k}$, one obtains:

[TABLE]

in which owing to (95) we also replaced $\breve{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}$ by $\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2}$. To find the explicit expression of $\mathbf{v}(\breve{\bm{\alpha}}_{k}|\mathbf{R}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2},\widehat{\mathbf{R}}_{\widetilde{\mathbf{x}}\widetilde{\mathbf{x}}}^{-2})$, involved we further denote:

[TABLE]

Then, using (92) and (174) in (177), it follows that:

[TABLE]

where $\widetilde{\mathbf{R}}^{[i]}_{ss}$ is given by (93).

Appendix E—Derivation of $\textrm{CRLB}(\bm{\eta})$

We have the following parameter vector:

[TABLE]

Therefore, the associated FIM can be written as:

[TABLE]

whose $ij$th entry is expressed as:

[TABLE]

with

[TABLE]

In (192), $\upsilon_{i}$ is the $i$th element of $\bm{\upsilon}$ and the involved partial derivatives of $\mathbf{R}_{\mathbf{x}\mathbf{x}}$ are given by:

[TABLE]

Furthermore, it can be shown that the partial derivatives of $\mathbf{R}^{\prime}_{\mathbf{x}\mathbf{x}}$ are given by:

[TABLE]

Recall that our goal is to find the CRLB of the angular parameters, $\bm{\eta}$, denoted as $\textrm{CRLB}(\bm{\eta})$. Therefore, we are interested in the $\bm{\eta}$-block of $\mathbf{I}^{-1}(\bm{\upsilon})$ only. From (187), the whole FIM, $\mathbf{I}(\bm{\upsilon})$, is a block matrix with $\mathbf{I}_{\bm{\eta},\bm{\eta}}$ being its first diagonal block. Thus, we use the block matrices inversion Lemma [References] to obtain the following expression for $\textrm{CRLB}(\bm{\eta})$:

[TABLE]