Cumulant-Based DOA Estimation of Noncircular Signals against Unknown Mutual Coupling

Abstract

To effectively find the direction of non-circular signals received by a uniform linear array (ULA) in the presence of non-negligible perturbations between array elements, i.e., mutual coupling, in colored noise, a direction of arrival (DOA) estimation approach in the context of high order statistics is proposed in this correspondence. Exploiting the non-circularity hidden behind a certain class of wireless communication signals to build up an augmented cumulant matrix, and carrying out a reformulation of the distorted steering vector to extract the angular information from the unknown mutual coupling, by exploiting the characteristic of mutual coupling, i.e., a limited operating range and an inverse relation of coupling effects to interspace, we develop a MUSIC-like estimator based on the rank-reduction (RARE) technique to directly determine directions of incident signals without mutual coupling compensation. Besides, we provide a solution to the problem of coherency between signals and mutual coupling between sensors co-existing, by selecting a middle sub-array to mitigate the undesirable effects and exploiting the rotation-invariant property to blindly separate the coherent signals into different groups to enhance the degrees of freedom. Compared with the existing robust DOA methods to the unknown mutual coupling under the framework of fourth-order cumulants (FOC), our work takes advantage of the larger virtual array and is able to resolve more signals due to greater degrees of freedom. Additionally, as the effective aperture is virtually extended, the developed estimator can achieve better performance under scenarios with high degree of mutual coupling between two sensors. Simulation results demonstrate the validity and efficiency of the proposed method.

1. Introduction

Direction of arrival (DOA) estimation, an important research area of sensor array signal processing, has attracted a large amount of attention because of its wide applications to electromagnetic, acoustic, seismic sensing, etc [1,2,3,4,5,6]. However, non-negligible interaction between sensors, i.e., mutual coupling, makes prevailing DOA estimation algorithms invalid since extra unknowns are introduced. Weiss and Friedlander [7] first discuss the structure of the mutual coupling matrix (MCM) in uniform linear and circular arrays and estimated the DOAs and mutual coupling coefficients in an iterative way. This work is followed by Sellone et al. whose approach is more robust to the detrimental effects and does not need a preliminary estimate [8]. Both schemes suffer from high computational complexity, which is inevitable given their reliance on multidimensional search processes. To circumvent the adverse effects of mutual coupling, the technique in [9,10] uses a subset of auxiliary sensors to make the “middle sub-array” mutual coupling free, and prevalent super-resolution algorithms can be directly applied the middle sub-array to resolve the DOAs. The main merit of these kind methods is the lower computational complexity due to not requiring any iterations but the degrees of freedom (DOFs) is reduced since only the middle sub-array is utilized. In a departure from middle sub-array methods [9,10], a recently proposed subspace-based method [11] takes advantage of the whole array by reparameterizing the actual array response and, hence improves the estimation accuracy. Recently, some novel sparsity inducing algorithms are developed to have robustness to the electromagnetic nuisance. Liu and Zhou first deal with the issue via the perspective of sparse Bayesian learning [12]. ${\ell }_1$-SVD is applied in [13] in conjunction with transforming the steering vector with mutual coupling, which is favorable for sparse recovery. Chen et al. investigate the problem of off-grid DOA estimation with coupling effects, and adopt the ${\ell }_p$-norm-based technique [14] and the relevance vector machine [15], respectively, to handle different cases, reducing amount of computation while preserving satisfactory estimation performance. Taking advantage of the uncorrelation between signals, Wang et al. address the issue from the perspective of group sparsity reconstruction of a long vector with mitigated noise components, and bring about an even better estimation accuracy and resolution [16].

The aforementioned work provides the solutions following the second-order statistics (SOS) and is unable to work properly in spatial colored noise if the noise covariance matrix is not available in advance. Since high-order statistics are insensitive to colored noise and can inherently enhance the DOFs [17,18,19], some study has been devoted to the suppression of unknown mutual coupling and colored noise using fourth-order cumulants (FOC) [20,21,22,23,24]. However, there are many additional pseudo-estimates in [20] which may affect the estimation performance, especially for the case where the true DOAs do not lie on the pseudo-peaks, while only the middle sub-array is utilized in [21] reducing the array aperture. Though promising results are shown in [22] and its follow up work [23,24], the array aperture in the former is fixed to a part of the physical aperture while the extended aperture is presented in the latter but still constraint to the middle sub-array, which gives rise to incomplete analysis and mining of statistical information implicated in the higher-order moments.

Recently, non-circularity embedded into wireless communication signals has been exploited to enhance DOFs for DOA estimation in the context of second-order statistics [25,26,27,28,29,30,31,32,33,34,35,36,37]. However, little work has been done on leveraging non-circularity in high-order moments to combat the unknown mutual coupling. To fill this gap, in this correspondence, a new FOC-based approach is proposed for the direction finding of non-circular signals in the presence of the deleterious effects between sensors. Making use of the characteristic of mutual coupling, such as a limited operating range and an inverse relation of coupling effects to interspace, to parameterize the steering vector, as well as the non-circularity of the observations to form a virtual array, we develop a MUSIC-like estimator by means of the rank-reduction (RARE) property to resolve the DOA estimates without any calibration process. Additionally, we introduce a solution to the problem of coherent signal estimation in the presence of mutual coupling between sensors, mitigating the detrimental effects via a middle sub-array and enhancing the DOFs through blind separation of the coherent signals into different groups. Compared with the existing FOC-based methods [21,22,23,24], in our proposed method, the DOFs as well as the effective aperture are further improved and hence, better accuracy and resolution can be obtained. Additionally, we discuss the identifiability of DOA estimation which has not been done in previous work, and the theoretical analysis shows that the work in this submission has an essential advantage in identifying more sources compared with [21,22,23,24] if the array has a relatively strong coupling between sensors.

The rest of the paper is organized as follows. In Section 2, the array model is discussed, when non-circular data observed by a sensor array with mutual coupling is constructed. In Section 3, we first develop a MUSIC-like method to directly determine the DOA estimates in conjunction with the non-circularity, but also isolate the mutual coupling effects, then discuss an even more challenging scenario where non-circular signals from one user in the uplink channel are coherent to each other, provide a solution to separate the coherent signals into different groups with the help of the middle array that is coupling-free, and then perform a DOA estimation in each coherent group via the forward–backward spatial smoothing (FBSS) technique. The estimation performance in terms of identifiability and computational complexity is analyzed in Section 4. Simulation results in Section 6 show the significant improvement of the proposed method. Finally, some concluding remarks are given in Section 7.

Throughout this submission, the following notations will be adopted: the operators ${(·)}^T$, ${(·)}^{*}$, ${(·)}^H$, ${(·)}^{+}$, $|·|$, $E[·]$, $det\{·\}$, ⊗, ∘, and ${\parallel ·\parallel }_2$ denote the operation of transpose, conjugate, conjugate transpose, pseudo-inverse, modulus, expectation, determinant, Kronecker product, Khatri–Rao (KR) product, and Euclidean (${\ell }_2$) norm, respectively. The symbol $\mathrm{diag}\{z_1,\cdots ,z_N\}$ represents a diagonal matrix with diagonal entries $z_1,\cdots ,z_N$, $\mathrm{blkdiag}\{{\mathbf{Z}}_1,{\mathbf{Z}}_2\}$ symbolizes a block diagonal matrix with diagonal entries ${\mathbf{Z}}_1$ and ${\mathbf{Z}}_2$, $\mathrm{Toeplitz}\left\{·\right\}$ represents a symmetric Toeplitz matrix constructed by a vector in the brace, and $\mathrm{cum}\{z_1,z_2,z_3,z_4\}$ symbolizes a cumulant calculated from data $z_1$, $z_2$, $z_3$, and $z_4$. The symbol ${\mathbf{I}}_K$ stands for an identity matrix of size $K\times K$. The symbol $\mathbf{Z}(a:b,c:d)$ refers to a constructed sub-matrix by the entries from a to b-th row and c to d-th column of $\mathbf{Z}$, and the symbol $\mathbf{Z}(a,b)$ denotes the entry in the a-th row and b-th column of $\mathbf{Z}$.

2. Problem Formulation

2.1. Strictly Second-Order Non-Circular Signals

For a signal $s\left(t\right)$, if its elliptic covariance is $E\left[s^2\left(t\right)\right]=\rho e^{j\varphi }{E\left[\right|s\left(t\right)|}^2]=\rho e^{j\varphi }{\sigma }_s^2\ne 0$, where $\varphi$ is the deterministic non-circularity phase, and $\rho$ is non-circularity rate satisfying $0<\rho \le 1$, then we say $s\left(t\right)$ is non-circular. The special case $\rho =1$ is common in wireless communications, and the received baseband signals satisfying such a condition is referred to as strictly second-order non-circular, like amplitude modulation (AM), binary phase shift keying (BPSK), amplitude shift keying (ASK), offset quadrature phase shift keying (OQPSK), and pulse amplitude modulation (PAM), etc. [25,26,27,28,30,38,39]. To aid the understanding of such a particular class of signals, we take the baseband signal of BPSK or ASK modulation as an example. In the framework of DOA estimation, the analog received BPSK or ASK modulated signals are bandpass filtered and down converted to baseband, the in-phase and quadrature components are matched-filtered, sampled, and paired to obtain complex signal $s\left(t\right)$ that can be factorized into $s\left(t\right)=e^{j\phi }{s}_0\left(t\right)$ where $s_0\left(t\right)$ is a real-valued symbol and $\phi$ is an arbitrary phase shift, due to the initial phase of the transmitted carrier, that can be different for each signal but constant with time. Therefore, the baseband signal $s\left(t\right)$ at the receiver side is complex-valued. This model has been used by numerous authors who have studied DOA estimation of non-circular signals [25,26,27,28,30,38,39]. Examining the covariance and elliptic covariance of $s\left(t\right)$, one gets $E\left[s\left(t\right)s^{*}\left(t\right)\right]=E\left[e^{j\phi }{s}_0\left(t\right)e^{-j\phi }{s}_0\left(t\right)\right]=E\left[s_0^2\left(t\right)\right]={\sigma }_s^2$ and $E\left[s\left(t\right)s\left(t\right)\right]=E\left[e^{j\phi }{s}_0\left(t\right)e^{j\phi }{s}_0\left(t\right)\right]=e^{j2\phi }E\left[s_0^2\left(t\right)\right]=e^{j2\phi }{\sigma }_s^2\ne 0$. By the definition of non-circular signal, one can readily find that $\rho =1$ and $\varphi =2\phi$, i.e., the baseband signal of BPSK or ASK modulation is strictly non-circular, and the non-circularity phase is twice as much as the initial phase. For the case of the received baseband signals with unknown modulation, one can resort to two kinds of classification approaches to identify the modulation of the signals. One class is the likelihood-based approach [40,41,42] and the other is the so-called feature-based approach [43,44,45,46,47,48]. If the signals are recognized as non-circular, such as BPSK, ASK, OQPSK, or PAM, the proposed method in the sequel is applicable, otherwise it fails to work. The research on modulation classification is beyond the scope of this manuscript, so we assume that the modulation of the received baseband signals is known a priori or can be recognized by using existing classification approaches if not specified, and the signals to be dealt with are strictly second-order non-circular.

2.2. Array Model for Non-Circular Signals

Consider narrowband non-circular signals $s_{0,i}\left(t\right)$ in the far filed impinge on an M-elements uniform linear array (ULA) from N different angles ${\theta }_i$, $i=1,2,\cdots ,N$, and the mutual coupling between sensors cannot be neglected, then at the time index t the corresponding array observation can be expressed as

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& =\sum _{i=1}^N\mathbf{C}\mathbf{a}\left({\theta }_i\right)e^{j\frac{{\varphi }_i}{2}}{s}_{0,i}\left(t\right)+\mathbf{n}\left(t\right)=\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)+\mathbf{n}\left(t\right)\hfill \end{array}$$

(1)

where $\mathbf{a}\left(\theta \right)={\left[\begin{array}{c}1,\beta \left(\theta \right),\cdots ,{\beta }^{M-1}\left(\theta \right)\end{array}\right]}^T\in {\mathbb{C}}^M$ is the steering vector, $\beta \left(\theta \right)=e^{j\frac{2\pi d}{\lambda }sin\theta }$, $\lambda$, and d are the carrier wavelength and the spacing between adjacent sensors, respectively, $\mathbf{C}$ denotes the MCM quantifying the degrees of electromagnetic coupling among elements, $\mathbf{A}=\left[\begin{array}{c}\mathbf{a}\left({\theta }_1\right),\cdots ,\mathbf{a}\left({\theta }_N\right)\end{array}\right]$ is the array manifold, $\mathsf{\Psi }=\mathrm{diag}\left\{e^{j\frac{{\varphi }_1}{2}},\cdots ,e^{j\frac{{\varphi }_N}{2}}\right\}$, ${\mathbf{s}}_0\left(t\right)={\left[\begin{array}{c}{s}_{0,1}\left(t\right),\cdots ,s_{0,N}\left(t\right)\end{array}\right]}^T\in {\mathbb{R}}^N$, and $\mathbf{n}\left(t\right)$ is the noise following Gaussian distribution but spatially colored. In addition, it is assumed that $\mathbf{A}$ is unambiguous, i.e., the steering vectors ${\left\{\mathbf{a}\left({\theta }_i\right)\right\}}_{i=1}^N$ are linearly independent for any set of distinct ${\left\{{\theta }_i\right\}}_{i=1}^N$.

As described in [9,10,11,49,50,51], one can sufficiently model the perturbed ULA by taking the electromagnetic characteristic of mutual coupling, such as a limited operating range and an inverse relation of coupling effects to interspace, into consideration. To be specific, If the mutual coupling acting on two sensors covers P interelement spacing at most such that the resultant MCM has finite non-zero coefficients, then the resultant MCM is banded symmetric Toeplitz, which can be formulated as

$$\begin{array}{cc}\hfill \mathbf{C}& =\mathrm{Toeplitz}\left\{\left[1,c_1,\cdots ,c_{P-1},{\mathbf{0}}_{1\times (M-P)}\right]\right\}\hfill \end{array}$$

(2)

where $0<|c_1|,|c_2|,\cdots ,|c_{P-1}|<c_0=1$ are the mutual coupling coefficients.

Since the incident signals follow strictly second-order non-circularity, they are deterministic and non-Gaussian, and the FOC matrices from the observation blocks can be given by

$$\begin{array}{cc}\hfill {\mathbf{C}}_{x1}& =\mathrm{cum}\left\{{\mathbf{x}}^{*}\left(t\right),{\mathbf{x}}^T\left(t\right),{\mathbf{x}}^{*}\left(t\right),{\mathbf{x}}^T\left(t\right)\right\}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathbf{C}}_{x2}& =\mathrm{cum}\left\{\mathbf{x}\left(t\right),{\mathbf{x}}^T\left(t\right),\mathbf{x}\left(t\right),{\mathbf{x}}^T\left(t\right)\right\}\hfill \end{array}$$

(4)

whose entries in the $\left[\left(k_1-1\right)M+k_2\right]$-th row and the $\left[\left(l_1-1\right)M+l_2\right]$-th column, $k_1,k_2,l_1,l_2=1,2,\cdots ,M$, are defined as

$$\begin{array}{cc}\hfill & {\mathbf{C}}_{x1}\left(\left(k_1-1\right)M+k_2,\left(l_1-1\right)M+l_2\right)\hfill \\ \hfill & =\mathrm{cum}\{x_{k_1}^{*}\left(t\right),x_{l_1}\left(t\right),x_{k_2}^{*}\left(t\right),x_{l_2}\left(t\right)\}\hfill \\ \hfill & =E\left[x_{k_1}^{*}\left(t\right)x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)x_{l_2}\left(t\right)\right]-E\left[x_{k_1}^{*}\left(t\right)x_{k_2}^{*}\left(t\right)\right]\hfill \\ \hfill & \times E\left[x_{l_1}\left(t\right)x_{l_2}\left(t\right)\right]-E\left[x_{k_1}^{*}\left(t\right)x_{l_1}\left(t\right)\right]E\left[x_{k_2}^{*}\left(t\right)x_{l_2}\left(t\right)\right]\hfill \\ \hfill & -E\left[x_{k_1}^{*}\left(t\right)x_{l_2}\left(t\right)\right]E\left[x_{k_2}^{*}\left(t\right)x_{l_1}\left(t\right)\right]\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\mathbf{C}}_{x2}\left(\left(k_1-1\right)M+k_2,\left(l_1-1\right)M+l_2\right)\hfill \\ \hfill & =\mathrm{cum}\{x_{k_1}\left(t\right),x_{l_1}\left(t\right),x_{k_2}\left(t\right),x_{l_2}\left(t\right)\}\hfill \\ \hfill & =E\left[x_{k_1}\left(t\right)x_{k_2}\left(t\right)x_{l_1}\left(t\right)x_{l_2}\left(t\right)\right]-E\left[x_{k_1}\left(t\right)x_{k_2}\left(t\right)\right]\hfill \\ \hfill & \times E\left[x_{l_1}\left(t\right)x_{l_2}\left(t\right)\right]-E\left[x_{k_1}\left(t\right)x_{l_1}\left(t\right)\right]E\left[x_{k_2}\left(t\right)x_{l_2}\left(t\right)\right]\hfill \\ \hfill & -E\left[x_{k_1}\left(t\right)x_{l_2}\left(t\right)\right]E\left[x_{k_2}\left(t\right)x_{l_1}\left(t\right)\right]\hfill \end{array}$$

(6)

where $x_m\left(t\right)$ is the m-th entry of $\mathbf{x}\left(t\right)$. Substituting Equation (1) into Equations (3) and (4), and utilizing the properties of high-order statistics [CP1]–[CP5] in [18], one can proceed to

$$\begin{array}{cc}\hfill {\mathbf{C}}_{x1}& =\mathrm{cum}\left\{{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^{*},{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^T,{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^{*},\right.\hfill \\ \hfill & \left.{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^T\right\}+\mathrm{cum}\left\{{\mathbf{n}}^{*}\left(t\right),{\mathbf{n}}^T\left(t\right),{\left(\mathbf{n}\left(t\right)\right)}^{*},{\left(\mathbf{n}\left(t\right)\right)}^T\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{p=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_p\right)\right)}^{*}{e}^{-j\frac{{\varphi }_p}{2}}{s}_{0,p}^{*}\left(t\right),\sum _{m=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_m\right)\right)}^T\right.\hfill \\ \hfill & \times e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{q=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_q\right)\right)}^{*}{e}^{-j\frac{{\varphi }_q}{2}}{s}_{0,q}^{*}\left(t\right),\hfill \\ \hfill & \left.\sum _{n=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_n\right)\right)}^T{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{p=1}^N\sum _{q=1}^N\sum _{m=1}^N\sum _{n=1}^N{\left(e^{j\frac{{\varphi }_p+{\varphi }_q}{2}}\left(\mathbf{C}\mathbf{a}\left({\theta }_p\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_q\right)\right)\right)}^{*}\left(e^{j\frac{{\varphi }_m+{\varphi }_n}{2}}\left(\mathbf{C}\mathbf{a}\left({\theta }_m\right)\right)\right.\hfill \\ \hfill & {\left.\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_n\right)\right)\right)}^T\mathrm{cum}\{s_{0,p}^{*}\left(t\right),s_{0,m}\left(t\right),s_{0,q}^{*}\left(t\right),s_{0,n}\left(t\right)\}\hfill \\ \hfill & =\sum _{i=1}^N{\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)}^{*}{\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)}^T\hfill \\ \hfill & \mathrm{cum}\{s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right)\}\hfill \\ \hfill & =\sum _{i=1}^N{\gamma }_i{\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)}^{*}{\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)}^T\hfill \\ \hfill & ={\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right)}^{*}{\mathbf{C}}_s{\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right)}^T\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\mathbf{C}}_{x2}& =\mathrm{cum}\left\{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right),{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^T,\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right),\right.\hfill \\ \hfill & \left.{\left(\mathbf{C}\mathbf{A}\mathsf{\Psi }{\mathbf{s}}_0\left(t\right)\right)}^T\right\}+\mathrm{cum}\left\{\mathbf{n}\left(t\right),{\mathbf{n}}^T\left(t\right),\left(\mathbf{n}\left(t\right)\right),{\left(\mathbf{n}\left(t\right)\right)}^T\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{p=1}^N\mathbf{C}\mathbf{a}\left({\theta }_p\right)e^{j\frac{{\varphi }_p}{2}}{s}_{0,p}\left(t\right),\sum _{m=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_m\right)\right)}^T\right.\hfill \\ \hfill & \times e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{q=1}^N\mathbf{C}\mathbf{a}\left({\theta }_q\right)e^{j\frac{{\varphi }_q}{2}}{s}_{0,q}\left(t\right),\hfill \\ \hfill & \left.\sum _{n=1}^N{\left(\mathbf{C}\mathbf{a}\left({\theta }_n\right)\right)}^T{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{p=1}^N\sum _{q=1}^N\sum _{m=1}^N\sum _{n=1}^N\left(e^{j\frac{{\varphi }_p+{\varphi }_q}{2}}\left(\mathbf{C}\mathbf{a}\left({\theta }_p\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_q\right)\right)\right)\hfill \\ \hfill & \times {\left(e^{j\frac{{\varphi }_m+{\varphi }_n}{2}}\left(\mathbf{C}\mathbf{a}\left({\theta }_m\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_n\right)\right)\right)}^T\mathrm{cum}\{s_{0,p}\left(t\right),\hfill \\ \hfill & s_{0,m}\left(t\right),s_{0,q}\left(t\right),s_{0,n}\left(t\right)\}\hfill \\ \hfill & =\sum _{i=1}^N{e}^{j2{\varphi }_i}\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right.\hfill \\ \hfill & {\left.\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)}^T\mathrm{cum}\{s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right)\}\hfill \\ \hfill & =\sum _{i=1}^N{e}^{j2{\varphi }_i}{\gamma }_i\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right){\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)}^T\hfill \\ \hfill & =\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right){\mathsf{\Psi }}^2{\mathbf{C}}_s{\mathsf{\Psi }}^2{\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right)}^T\hfill \end{array}$$

(8)

where ${\gamma }_i\triangleq cum\left\{s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right)\right\}$ and ${\mathbf{C}}_s\triangleq \mathrm{diag}\left\{{\gamma }_1,\cdots ,{\gamma }_N\right\}\in {\mathbb{R}}^{N\times N}$.

3. Proposed Non-Circular FOC-Based Estimator

To increase the effective aperture as well as the DOFs, one can design an augmented cumulant matrix, of size $2M^2\times 2M^2$, as

$$\begin{array}{cc}\hfill {\mathbf{C}}_x& =\left[\begin{array}{cc}{\mathbf{C}}_{x1}& {\mathbf{C}}_{x2}^{*}\\ {\mathbf{C}}_{x2}& {\mathbf{C}}_{x1}^{*}\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}{\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right)}^{*}{\mathsf{\Psi }}^{-2}\\ \left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right){\mathsf{\Psi }}^2\end{array}\right]{\mathbf{C}}_s{\left[\begin{array}{c}{\left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right)}^{*}{\mathsf{\Psi }}^{-2}\\ \left(\left(\mathbf{C}\mathbf{A}\right)\circ \left(\mathbf{C}\mathbf{A}\right)\right){\mathsf{\Psi }}^2\end{array}\right]}^H.\hfill \end{array}$$

(9)

DOA Estimation Without Mutual Coupling Compensation

The augmented matrix ${\mathbf{C}}_x$ has the following factorization through the singular value decomposition (SVD)

$$\begin{array}{cc}\hfill {\mathbf{C}}_x& =\mathbf{U}\mathsf{\Sigma }{\mathbf{V}}^H\hfill \end{array}$$

(10)

where $\mathsf{\Sigma }=\mathrm{diag}\left\{{\lambda }_1,\cdots ,{\lambda }_{2M^2}\right\}$ is composed of $2M^2$ singular values that can be sorted in descending order as ${\lambda }_1\ge \cdots \ge {\lambda }_N>{\lambda }_{N+1}=\cdots ={\lambda }_{2M^2}=0$. The matrix ${\mathbf{U}}_s\triangleq \mathbf{U}(:,1:N)$ collects the singular vectors corresponding to the N largest singular values, while ${\mathbf{U}}_n\triangleq \mathbf{U}(:,N+1:2M^2)$ contains the rest of the singular vectors corresponding to the $2M^2-N$ zero singular values. It is known that ${\mathbf{U}}_s$ spans the signal subspace, while its orthogonal space, namely the noise subspace, is spanned by ${\mathbf{U}}_n$, thus one has

$$\begin{array}{cc}\hfill & {∥{\left[\begin{array}{c}{e}^{-j{\varphi }_i}{\left(\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\right)}^{*}\\ e^{j{\varphi }_i}\left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\otimes \left(\mathbf{C}\mathbf{a}\left({\theta }_i\right)\right)\end{array}\right]}^H{\mathbf{U}}_n∥}_2^2=0,\mathrm{for}i=1,2,\cdots ,N.\hfill \end{array}$$

(11)

However, due to the unknown $\mathbf{C}$ one cannot directly obtain the DOA estimates of the desired non-circular signals via 1-D search for making Equation (11) hold. For the purpose of decoupling the angular information from the electromagnetic impact between sensors of the ULA, we reparameterize the distorted steering vector as [11]

$$\begin{array}{cc}\hfill \mathbf{C}\mathbf{a}\left(\theta \right)& =\mathbf{T}\left(\theta \right)\mathit{\alpha }\hfill \end{array}$$

(12)

where

$$\begin{array}{cc}\hfill \mathbf{T}\left(\theta \right)& =\mathrm{blkdiag}\left\{{\mathbf{T}}_1,{\mathbf{T}}_2,{\mathbf{T}}_3\right\}\in {\mathbb{C}}^{M\times (2P-1)}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \mathit{\alpha }& ={\left[\begin{array}{c}{\mu }_1,\cdots ,{\mu }_{P-1},\tau \left(\theta \right),{\alpha }_1,\cdots ,{\alpha }_{P-1}\end{array}\right]}^T\in {\mathbb{C}}^{2P-1}\hfill \end{array}$$

(14)

with

$$\begin{array}{cc}\hfill {\mathbf{T}}_1& =\mathrm{diag}\left\{1,\beta \left(\theta \right),\cdots ,{\beta }^{P-2}\left(\theta \right)\right\}\in {\mathbb{C}}^{(P-1)\times P-1)}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\mathbf{T}}_2& ={\left[{\beta }^{P-1}\left(\theta \right),\cdots ,{\beta }^{M-P}\left(\theta \right)\right]}^T\in {\mathbb{C}}^{M-2P+2}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathbf{T}}_3& =\mathrm{diag}\left\{{\beta }^{M-P+1}\left(\theta \right),\cdots ,{\beta }^{M-1}\left(\theta \right)\right\}\in {\mathbb{C}}^{(P-1)\times P-1)}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\mu }_l& =1+\sum _{k=1}^{P-1}{c}_k{\beta }^k\left(\theta \right)+\sum _{k=1}^{l-1}{c}_k{\beta }^{-k}\left(\theta \right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\alpha }_l& =1+\sum _{k=1}^{P-1}{c}_k{\beta }^{-k}\left(\theta \right)+\sum _{k=1}^{P-1-l}{c}_k{\beta }^k\left(\theta \right)\hfill \end{array}$$

(19)

and

$$\begin{array}{cc}\hfill \tau \left(\theta \right)& =1+\sum _{k=1}^{P-1}{c}_k\left({\beta }^k\left(\theta \right)+{\beta }^{-k}\left(\theta \right)\right).\hfill \end{array}$$

(20)

It should be noted that $\tau \left(\theta \right)$ in Equation (14) is generally not equal to zero, otherwise “blind angles” cannot be identified [11]. Substituting Equation (12) back to Equation (11), one has

$$\begin{array}{cc}\hfill & {∥{\left[\begin{array}{c}{e}^{-j{\varphi }_i}{\left(\left(\mathbf{T}\left({\theta }_i\right)\alpha \right)\otimes \left(\mathbf{T}\left({\theta }_i\right)\alpha \right)\right)}^{*}\\ e^{j{\varphi }_i}\left(\mathbf{T}\left({\theta }_i\right)\alpha \right)\otimes \left(\mathbf{T}\left({\theta }_i\right)\alpha \right)\end{array}\right]}^H{\mathbf{U}}_n∥}_2^2\hfill \\ \hfill & ={∥{\left[\begin{array}{c}{e}^{-j{\varphi }_i}{\left(\left(\mathbf{T}\left({\theta }_i\right)\otimes \mathbf{T}\left({\theta }_i\right)\right)\left(\alpha \otimes \alpha \right)\right)}^{*}\\ e^{j{\varphi }_i}\left(\mathbf{T}\left({\theta }_i\right)\otimes \mathbf{T}\left({\theta }_i\right)\right)\left(\alpha \otimes \alpha \right)\end{array}\right]}^H{\mathbf{U}}_n∥}_2^2\hfill \\ \hfill & ={∥{\upsilon }^H{\overline{\mathbf{T}}}^H\left({\theta }_i\right){\mathbf{U}}_n∥}_2^2\hfill \\ \hfill & ={\upsilon }^H{\overline{\mathbf{T}}}^H\left({\theta }_i\right){\mathbf{U}}_n{\mathbf{U}}_n^H\overline{\mathbf{T}}\left({\theta }_i\right)\upsilon \hfill \\ \hfill & ={\upsilon }^H\mathbf{M}\left({\theta }_i\right)\upsilon =0\hfill \end{array}$$

(21)

where $\mathit{\upsilon }\triangleq {\left[e^{-j{\varphi }_i}{\left(\mathit{\alpha }\otimes \mathit{\alpha }\right)}^H,e^{j{\varphi }_i}{\left(\mathit{\alpha }\otimes \mathit{\alpha }\right)}^T\right]}^T\in {\mathbb{C}}^{2{(2P-1)}^2}$, $\overline{\mathbf{T}}\left({\theta }_i\right)\triangleq \mathrm{blkdiag}\left\{{\left(\mathbf{T}\left({\theta }_i\right)\otimes \mathbf{T}\left({\theta }_i\right)\right)}^{*},\mathbf{T}\left({\theta }_i\right)\otimes \mathbf{T}\left({\theta }_i\right)\right\}\in {\mathbb{C}}^{2M^2\times 2{(2P-1)}^2}$, and $\mathbf{M}\left({\theta }_i\right)\triangleq {\overline{\mathbf{T}}}^H\left({\theta }_i\right){\mathbf{U}}_n{\mathbf{U}}_n^H\overline{\mathbf{T}}\left({\theta }_i\right)\in {\mathbb{C}}^{2{(2P-1)}^2\times 2{(2P-1)}^2}$.

It is clear that ${\overline{\mathbf{T}}}^H\left(\theta \right){\mathbf{U}}_n$ is of size $2{(2P-1)}^2\times (2M^2-N)$. If $2{(2P-1)}^2\le 2M^2-N$, generally speaking, then the matrix ${\overline{\mathbf{T}}}^H\left(\theta \right){\mathbf{U}}_n$ has a full row rank and $\mathbf{M}\left(\theta \right)$ has a full rank accordingly. Nevertheless, for some special cases that $\theta$ is exactly the same with any one of the N true angles, i.e., $\theta ={\theta }_i,i=1,\cdots ,N$, the expression in Equation (21) becomes zero. Since $\mathit{\upsilon }\ne \mathbf{0}$, Equation (21) is valid only when $\mathbf{M}\left(\theta \right)$ is singular that it takes as a zero-value determinant. As a result, the DOA estimation is now dependent on the determinant of $\mathbf{M}\left(\theta \right)$. Since $\mathbf{M}\left(\theta \right)$ does not contain any mutual coupling coefficients, $\det \left\{\mathbf{M}\left(\theta \right)\right\}$ is insensitive to the nuisances. Therefore, one can estimate the bearings as

$$\begin{array}{cc}\hfill \widehat{\theta }& =arg\underset{\theta }{min}det\left\{\mathbf{M}\left(\theta \right)\right\}.\hfill \end{array}$$

(22)

4. Performance Analysis

In this section, the identifiability as well as the computational complexity of the proposed method is discussed, with a comparison of the existing work in the context of FOC, e.g., Liu’s [21] and Liao’s [24] approaches.

4.1. Identifiability of DOA Estimation

As $2{(2P-1)}^2\le 2M^2-N$ makes $\mathbf{M}\left(\theta \right)$ be of full rank, which is a necessary condition of the proposed algorithm, one can deduce that $N\le 2M^2-2{(2P-1)}^2$, which implies that the maximum DOFs of the developed estimator are $2M^2-2{(2P-1)}^2$. By contrast, existing methods cannot estimate as many sources as our technique. Liao’s method is based on ESPRIT and as a result, the upper bound of estimation identifiability depends on the dimension of signal subspace. It is known that the signal subspace produced by Liao’s method has a size of $M(M-2P+1)\times N$, and there are also $(M-2P+2)(M-2P+1)-2(M-2P+1)=(M-2P)(M-2P+1)$ linearly dependent rows due to the property of the middle sub-array. Therefore, the row dimension of the signal subspace is $2P(M-2P+1)$, in other words, Liao’s method can estimate up to $2P(M-2P+1)$ signals. In Liu’s method, only the middle sub-array utilized, which means the extended aperture by FOC is $2(M-2P+2)-1$ elements. Therefore, Liu’s method can handle at most $2(M-2P+1)$ DOAs. It is evident that the proposed method has the largest estimation capacity among the three algorithms, followed by Liao’s method and then Liu’s method. Consequently, Liu’s method suffers from the reduced effective aperture. Besides, if the coupling between sensors is prominent, i.e., P is large, or equivalently M is not sufficiently large, then the size of the middle sub-array will be very small. For instance, if $P=\frac{M}{2}$ and $N=4$, Liu’s method is not valid. This is because the effective array aperture [18] is 3 elements in this case, and at most 2 DOAs can be estimated.

4.2. Computational Complexity

Next, we compare the major computations of the proposed estimator, Liao’s and Liu’s methods, which involve in the fourth-order statistical matrices construction, SVD, and spectral search. For Liao’s approach, the main calculation cost is to form $M-2P+2$ FOC matrices of the same size $M\times M$ and to perform the SVD of the augmented $M(M-2P+2)\times M$ cumulant matrix, which requires $\mathcal{O}\left(9(M-2P+2)M^{2}L+M^3(M-2P+2)\right)$ flops where L is the number of snapshots. Here, a flop stands for a complex-valued floating point multiplication operation. Liu’s technique constructs one ${(M-2P+2)}^2\times {(M-2P+2)}^2$ cumulant matrix, implements its SVD, and one-dimensional MUSIC spectral search, respectively, and it thus needs $\mathcal{O}\left(9{(M-2P+2)}^{4}L+{(M-2P+2)}^6+\frac{180}{\delta }({(M-2P+2)}^2+1)({(M-2P+2)}^2-N)\right)$ flops in total where $\delta$ is the sampling grid spacing. By comparison, the proposed scheme establishes two $M^2\times M^2$ cumulant matrices, carries out SVDs of the augmented $2M^2\times 2M^2$ matrix, and executes a one-dimensional spectral search on determinants of a $2{\left(2P-1\right)}^2$ matrix. The resulting flops taken in our work are in order of $\mathcal{O}\left(18M^{4}L+8M^6+\frac{1440}{\delta }{\left(2P-1\right)}^6\right)$. As in general L is relatively large, say several thousands, the proposed method has the highest computational complexity among the three algorithms, followed by Liao’s approach and then Liu’s technique.

5. Solution to the Case of Coherent Signals

It should be noted that the aforementioned method is applicable to incident signals being independent, and the coherency between signals may lead to severe performance degradation or even total failure due to the rank deficiency. This motivates us to provide a solution to the knotty problem, i.e., multiple groups of coherent signals received by a tightly coupled ULA in colored noise, and another direction finding algorithm based on FOC is proposed in this section, exploiting rotation-invariant property to blindly separate the coherent signals into different groups, and then resolving the DOAs in each coherent group via rank restoration techniques such as FBSS.

First, in order to mitigate the unknown mutual coupling, one can extract the observations of the middle sub-array by $\tilde{\mathbf{x}}\left(t\right)=\mathbf{W}\mathbf{x}\left(t\right)$ where the selection matrix $\mathbf{W}=\left[\begin{array}{c}{\mathbf{0}}_{(M-2P+2)\times (P-1)},{\mathbf{I}}_{M-2P+2},{\mathbf{0}}_{(M-2P+2)\times (P-1)}\end{array}\right]$.

Consider that there are K groups of signals from K users sensed by the sensor array deployed at the base station where the k-th group has $P_k$ signals that are coherent to each other, but the signals from different groups are independent. Referring to [52,53,54], the array observation with a shrinking size can be rewritten as

$$\begin{array}{cc}\hfill \overline{\mathbf{x}}\left(t\right)& =\overline{\mathbf{A}}\tilde{\mathsf{\Gamma }}\mathbf{s}\left(t\right)+\overline{\mathbf{n}}\left(t\right)\hfill \end{array}$$

(23)

where $\overline{\mathbf{A}}=\mathbf{J}\mathbf{A}$, $\tilde{\mathsf{\Gamma }}=\mathbf{D}\mathsf{\Gamma }$ with $\mathbf{D}=\mathrm{diag}\{\mu \left({\theta }_1\right),\mu \left({\theta }_2\right),\cdots ,\mu \left({\theta }_N\right)\},\mu \left({\theta }_i\right)=c_{P-1}+\cdots +{\beta }^{P-1}\left({\theta }_i\right)+\cdots +c_{P-1}{\beta }^{2P-2}\left({\theta }_i\right)$, and $\mathsf{\Gamma }=\mathrm{blkdiag}\{{\mathit{\alpha }}_1,\cdots ,{\mathit{\alpha }}_K\}$ with ${\mathit{\alpha }}_k={\left[\begin{array}{c}{\alpha }_{k1},\cdots ,{\alpha }_{k{P}_k}\end{array}\right]}^T$ collecting fading coefficients.

Let $\tilde{\mathbf{A}}=\overline{\mathbf{A}}\tilde{\mathsf{\Gamma }}=\left[{\mathbf{b}}_1,\cdots ,{\mathbf{b}}_K\right]$, being the so-called generalized array manifold, where ${\mathbf{b}}_k=\left[\mathbf{a}\left({\theta }_{k1}\right),\cdots ,\mathbf{a}\left({\theta }_{kPk}\right)\right]{\tilde{\mathit{\alpha }}}_k$ is the generalized steering vector, one can get four cumulant matrices as follows

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_1& \triangleq \mathrm{cum}\left\{{\overline{\mathbf{x}}}_1\left(t\right),{\overline{\mathbf{x}}}_1^{*}\left(t\right),\overline{\mathbf{x}}\left(t\right),{\overline{\mathbf{x}}}^H\left(t\right)\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{m=1}^K\tilde{\mathbf{A}}(1,m)e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{p=1}^K{\tilde{\mathbf{A}}}^{*}(1,p)e^{-j\frac{{\varphi }_p}{2}}{s}_{0,p}^{*}\left(t\right),\sum _{n=1}^K{\mathbf{b}}_n{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right),\sum _{r=1}^K{\mathbf{b}}_r^H{e}^{-j\frac{{\varphi }_r}{2}}{s}_{0,r}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{m=1}^K\sum _{n=1}^K\sum _{p=1}^K\sum _{r=1}^K{e}^{j\frac{{\varphi }_m-{\varphi }_p+{\varphi }_n-{\varphi }_r}{2}}\tilde{\mathbf{A}}(1,m){\tilde{\mathbf{A}}}^{*}(1,p){\mathbf{b}}_n{\mathbf{b}}_r^H\mathrm{cum}\left\{s_{0,m}\left(t\right),s_{0,p}^{*}\left(t\right),s_{0,n}\left(t\right),s_{0,r}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^H\mathrm{cum}\left\{s_{0,i}\left(t\right),s_{0,i}^{*}\left(t\right),s_{0,i}\left(t\right),s_{0,i}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{\gamma }_i{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^H\hfill \\ \hfill & =\tilde{\mathbf{A}}{\tilde{\mathbf{C}}}_s{\tilde{\mathbf{A}}}^H\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_2& \triangleq \mathrm{cum}\left\{{\overline{\mathbf{x}}}_2\left(t\right),{\overline{\mathbf{x}}}_1^{*}\left(t\right),\overline{\mathbf{x}\left(t\right)},{\overline{\mathbf{x}}}^H\left(t\right)\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{m=1}^K\tilde{\mathbf{A}}(2,m)e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{p=1}^K{\tilde{\mathbf{A}}}^{*}(1,p)e^{-j\frac{{\varphi }_p}{2}}{s}_{0,p}^{*}\left(t\right),\sum _{n=1}^K{\mathbf{b}}_n{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right),\sum _{r=1}^K{\mathbf{b}}_r^H{e}^{-j\frac{{\varphi }_r}{2}}{s}_{0,r}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{m=1}^K\sum _{n=1}^K\sum _{p=1}^K\sum _{r=1}^K{e}^{j\frac{{\varphi }_m-{\varphi }_p+{\varphi }_n-{\varphi }_r}{2}}\tilde{\mathbf{A}}(2,m){\tilde{\mathbf{A}}}^{*}(1,p){\mathbf{b}}_n{\mathbf{b}}_r^H\mathrm{cum}\left\{s_{0,m}\left(t\right),s_{0,p}^{*}\left(t\right),s_{0,n}\left(t\right),s_{0,r}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K\frac{\tilde{\mathbf{A}}(2,m)}{\tilde{\mathbf{A}}(1,m)}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^H\mathrm{cum}\left\{s_{0,i}\left(t\right),s_{0,i}^{*}\left(t\right),s_{0,i}\left(t\right),s_{0,i}^{*}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{\gamma }_i\frac{\tilde{\mathbf{A}}(2,m)}{\tilde{\mathbf{A}}(1,m)}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^H\hfill \\ \hfill & =\tilde{\mathbf{A}}\tilde{\mathbf{D}}{\tilde{\mathbf{C}}}_s{\tilde{\mathbf{A}}}^H\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_3& \triangleq \mathrm{cum}\left\{{\overline{\mathbf{x}}}_1\left(t\right),{\overline{\mathbf{x}}}_1^{*}\left(t\right),\overline{\mathbf{x}}\left(t\right),{\overline{\mathbf{x}}}^T\left(t\right)\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{m=1}^K\tilde{\mathbf{A}}(1,m)e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{p=1}^K{\tilde{\mathbf{A}}}^{*}(1,p)e^{-j\frac{{\varphi }_p}{2}}{s}_{0,p}\left(t\right),\sum _{n=1}^K{\mathbf{b}}_n{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right),\sum _{r=1}^K{\mathbf{b}}_r^T{e}^{j\frac{{\varphi }_r}{2}}{s}_{0,r}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{m=1}^K\sum _{n=1}^K\sum _{p=1}^K\sum _{r=1}^K{e}^{j\frac{{\varphi }_m-{\varphi }_p+{\varphi }_n+{\varphi }_r}{2}}\tilde{\mathbf{A}}(1,m){\tilde{\mathbf{A}}}^{*}(1,p){\mathbf{b}}_n{\mathbf{b}}_r^T\mathrm{cum}\left\{s_{0,m}\left(t\right),s_{0,p}\left(t\right),s_{0,n}\left(t\right),s_{0,r}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{e}^{j{\varphi }_i}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^T\mathrm{cum}\left\{s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{\gamma }_i{e}^{j{\varphi }_i}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^T\hfill \\ \hfill & =\tilde{\mathbf{A}}\mathsf{\Psi }{\tilde{\mathbf{C}}}_s\mathsf{\Psi }{\tilde{\mathbf{A}}}^T\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_4& \triangleq \mathrm{cum}\left\{{\overline{\mathbf{x}}}_2\left(t\right),{\overline{\mathbf{x}}}_1^{*}\left(t\right),\overline{\mathbf{x}}\left(t\right),{\overline{\mathbf{x}}}^T\left(t\right)\right\}\hfill \\ \hfill & =\mathrm{cum}\left\{\sum _{m=1}^K\tilde{\mathbf{A}}(2,m)e^{j\frac{{\varphi }_m}{2}}{s}_{0,m}\left(t\right),\sum _{p=1}^K{\tilde{\mathbf{A}}}^{*}(1,p)e^{-j\frac{{\varphi }_p}{2}}{s}_{0,p}\left(t\right),\sum _{n=1}^K{\mathbf{b}}_n{e}^{j\frac{{\varphi }_n}{2}}{s}_{0,n}\left(t\right),\sum _{r=1}^K{\mathbf{b}}_r^T{e}^{j\frac{{\varphi }_r}{2}}{s}_{0,r}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{m=1}^K\sum _{n=1}^K\sum _{p=1}^K\sum _{r=1}^K{e}^{j\frac{{\varphi }_m-{\varphi }_p+{\varphi }_n+{\varphi }_r}{2}}\tilde{\mathbf{A}}(2,m){\tilde{\mathbf{A}}}^{*}(1,p){\mathbf{b}}_n{\mathbf{b}}_r^T\mathrm{cum}\left\{s_{0,m}\left(t\right),s_{0,p}\left(t\right),s_{0,n}\left(t\right),s_{0,r}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{e}^{j{\varphi }_i}\frac{\tilde{\mathbf{A}}(2,m)}{\tilde{\mathbf{A}}(1,m)}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^T\mathrm{cum}\left\{s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right),s_{0,i}\left(t\right)\right\}\hfill \\ \hfill & =\sum _{i=1}^K{\gamma }_i{e}^{j{\varphi }_i}\frac{\tilde{\mathbf{A}}(2,m)}{\tilde{\mathbf{A}}(1,m)}{\left|\tilde{\mathbf{A}}(1,i)\right|}^2{\mathbf{b}}_i{\mathbf{b}}_i^T\hfill \\ \hfill & =\tilde{\mathbf{A}}\mathsf{\Psi }\tilde{\mathbf{D}}{\tilde{\mathbf{C}}}_s\mathsf{\Psi }{\tilde{\mathbf{A}}}^T\hfill \end{array}$$

(27)

where ${\overline{\mathbf{x}}}_1\left(t\right)$ and ${\overline{\mathbf{x}}}_2\left(t\right)$ are the first and second columns of $\overline{\mathbf{x}}\left(t\right)$, respectively.

Then, one can take the eigen-decomposition to ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}$ and ${\mathsf{\Xi }}_4{\mathsf{\Xi }}_3^{+}$ respectively to work out the scaled generalized steering vectors. To see how it works, one can perform the following operations:

$$\begin{array}{cc}\hfill {\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}\tilde{\mathbf{A}}& =\tilde{\mathbf{A}}\tilde{\mathbf{D}}{\tilde{\mathbf{C}}}_s{\tilde{\mathbf{A}}}^H{\left(\tilde{\mathbf{A}}{\tilde{\mathbf{C}}}_s{\tilde{\mathbf{A}}}^H\right)}^{+}\tilde{\mathbf{A}}\hfill \\ \hfill & =\tilde{\mathbf{A}}\tilde{\mathbf{D}}{\tilde{\mathbf{C}}}_s{\tilde{\mathbf{A}}}^H{\left({\tilde{\mathbf{A}}}^H\right)}^{+}{\tilde{\mathbf{C}}}_s^{-1}{\tilde{\mathbf{A}}}^{+}\tilde{\mathbf{A}}\hfill \\ \hfill & =\tilde{\mathbf{A}}\tilde{\mathbf{D}}{\tilde{\mathbf{C}}}_s{\left({\tilde{\mathbf{A}}}^{+}\tilde{\mathbf{A}}\right)}^H{\tilde{\mathbf{C}}}_s^{-1}\hfill \\ \hfill & =\tilde{\mathbf{A}}\tilde{\mathbf{D}}\hfill \end{array}$$

(28)

where in the third equality we have used the identity ${\tilde{\mathbf{A}}}^{+}\tilde{\mathbf{A}}={\mathbf{I}}_K$ since $\tilde{\mathbf{A}}$ is of full column rank. Through eigen-decomposition, one has ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}={\tilde{\mathbf{E}}}_s{\tilde{\mathsf{\Sigma }}}_s{\tilde{\mathbf{E}}}_s^H$, where ${\tilde{\mathbf{E}}}_s$ contains the eigenvectors of ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}$ while ${\tilde{\mathsf{\Sigma }}}_s$ is a diagonal matrix comprising eigenvalues. Equivalently, one can deduce that ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}{\tilde{\mathbf{E}}}_s={\tilde{\mathbf{E}}}_s{\tilde{\mathsf{\Sigma }}}_s$. As a result, it is natural to reveal ${\tilde{\mathbf{E}}}_s=\tilde{\mathbf{A}}\mathbf{\Lambda }$, where $\mathbf{\Lambda }$ is an arbitrary diagonal matrix with nonzero entries, and ${\tilde{\mathsf{\Sigma }}}_s=\tilde{\mathbf{D}}$. Following the same principle, one can also obtain the eigenvectors of ${\mathsf{\Xi }}_4{\mathsf{\Xi }}_3^{+}$ as the scaled generalized steering vectors. Denoting ${\tilde{\mathbf{u}}}_k$ and ${\breve{\mathbf{u}}}_k$ are the the eigenvectors of ${\mathsf{\Xi }}_2{\mathsf{\Xi }}_1^{+}$ and ${\mathsf{\Xi }}_4{\mathsf{\Xi }}_3^{+}$, respectively, one can get a more robust estimate of the scaled generalized steering vectors by averaging them, i.e.,

$$\begin{array}{cc}\hfill {\tilde{\mathbf{b}}}_k& \triangleq \beta {\mathbf{b}}_k=\frac{1}{2}\left({\tilde{\mathbf{u}}}_k+{\breve{\mathbf{u}}}_k\right).\hfill \end{array}$$

(29)

Then one can perform FBSS to ${\tilde{\mathbf{b}}}_k$ to recover the rank deficiency

$$\begin{array}{cc}\hfill {\mathbf{B}}_k^{fb}& =\frac{1}{2q}\sum _{i=1}^q\left[{\mathbf{W}}_i{\tilde{\mathbf{b}}}_k{\tilde{\mathbf{b}}}_k^H{\mathbf{W}}_i^H+\mathbf{J}{\left({\mathbf{W}}_i{\tilde{\mathbf{b}}}_k{\tilde{\mathbf{b}}}_k^H{\mathbf{W}}_i^H\right)}^{*}\mathbf{J}\right]\hfill \end{array}$$

(30)

where ${\mathbf{W}}_i=\left[{\mathbf{0}}_{m\times (i-1)},{\mathbf{I}}_m,{\mathbf{0}}_{m\times (q-i)}\right]$ is the selection matrix for the i-th sub-array, with $m=M-2P+3-q$ being is the number of sensors in each sub-array, and $\mathbf{J}\in {\mathbb{R}}^{m\times m}$ is the exchange matrix. For the k-th group of coherent signals, one now can apply ESPRIT to ${\mathbf{B}}_k^{fb}$ to determine the DOA estimates.

Remark 1.

The parameter q plays a significant role in rank restoration since $rank\left({\mathit{B}}_k^{fb}\right)$ after FBSS becomes $1+2q$, and it also restricts the DOFs and effective aperture for estimation by $m=M-2P+3-q$. We consider two extreme cases to discuss the choices of q. If $q=⌈\frac{P_i}{2}⌉$, $i=1,2\cdots ,K$, i.e., q achieves its lower bound, then one barely restores the rank deficiency with the maximum number of sensors after FBSS, but this may cause not all coherent signals to be detected, especially at low signal-to-noise ratio (SNR) or for few snapshots. On the other hand, If $q=M-2P+2-P_i$, i.e., q achieves its upper bound, then one restores the rank deficiency with an excessive number of times of FBSS, but this may cause biased estimates due to only one dimensional noise subspace being available. To our empirical knowledge and simulation results, the proposed method performs well when q is selected appropriately between the two bounds. However, the optimum choice of q is still an open problem and no theoretical guidance for the selection has been provided. A plausible way to choose the “optimal” q is to check whether ${\mathit{B}}_k^{fb}$ induced by q has the largest discrepancy between the $P_i$ largest eigenvalues and the $m-P_i$ smallest eigenvalues, in which case ${\mathit{B}}_k^{fb}$ can be considered to have the largest SNR. Therefore, the optimal q can be obtained by finding ${\mathit{B}}_k^{fb}$ with the largest discrepancy between the signal and the noise subspaces. Although this method for choosing optimal q is completely ad-hoc, it seems to be very simple and effective from simulations.

Remark 2.

A cumulant requires many more computations than a covariance and, hence, the proposed FOC-based estimator of non-circular signals against mutual coupling works well on condition that the observation window is sufficiently long and the target is (quasi-) stationary. We attempted to implement the proposed method in a testbed composed of a Virtex-7 series FPGA and a TMS320C6x series DSP, but found that it is infeasible to complete the DOA estimation within a time less than the scale of milliseconds in the case of $M=8$ and $L=2000$ due to the cumbersome calculations of cumulants and the SVD of the augmented cumulant matrix ${\mathit{C}}_x$, of size $16\times 16$. As a result, for some applications, like maneuverable targets, the coherence time of such cumulants is quite short, and it is hard and even unlikely to realize the proposed algorithm by the prevailing hardware available in the market. On the other hand, for some applications that expect not much real-time quality, such as radio, hydrological, or meteorological environment monitoring where the incident signals have cyclostationarity, cumulant-based signal selective algorithms can be implemented for location estimation of far-field signals; the price to be paid is the need for the large number of computations and large data lengths for reliable estimation of the cumulants.

Remark 3.

For the case of wideband signals, one can divide the observations at each channel into some (possibly overlapping) segments, where for each segment, a number of frequency sub-bands are computed by the short-time Fourier transform (STFT). The idea of frequency-domain processing is to decouple the wideband model into a multitude of narrowband models, and then at each frequency subband one can apply the proposed method to obtain $\mathit{M}(f,\theta )$ constructed by Equation (21). In the frequency-domain approach, the final step is to combine $\mathit{M}(f,\theta )$ at various frequencies to obtain a DOA spectrum fusion, that is, $\widehat{\theta }=arg\underset{\theta }{min}{\sum }_{f\in {\mathcal{B}}_d}det\left\{\mathit{M}\left(\theta \right)\right\}$ where ${\mathcal{B}}_d\subset \left[0,\frac{1}{2}\right]$ is the discrete normalized frequency band of the received signals. The combination in the spectrum fusion equation follows the principle of incoherent signal subspace method (ISSM) in [55]. By contrast, our solution is inapplicable to the focusing-based approaches [56,57], which is another category of DOA estimation algorithms for wideband signals, since initial DOA estimates required therein cannot be obtained by conventional beamforming in the presence of unknown mutual coupling.

Remark 4.

The array perturbations, such as mutual coupling, gain-phase errors, and sensor location uncertainties, can be calibrated by placing signal sources at known positions, which is the so-called active calibration [58,59,60,61,62,63,64,65]. However, the process of measuring the array manifold including various perturbations can be time consuming and expensive, and it is inconvenient or even infeasible to make the calibration source available in some cases. In addition to the deployment of signals of opportunity, keeping the calibration effective is another essential issue. Plenty of factors contribute to the variation of array response over time: gradual changes in the behaviors of the sensor itself, the electronic circuitry, and the analog-to-digital converter (due to thermal effects, aging of components, etc.), changes in the electromagnetic environment (e.g., metal objects beside an antenna array cause a distortion of the beam pattern), and changes in the sensor locations (e.g., an antenna array mounted on the vibrating wing of an aircraft or a hydrophone array towed behind a ship). The active calibration scheme has its inherent shortcomings that are difficult to be overcome and, hence, in this paper we resort to the self-calibration solution that is not reliant on the array manifold measurement as well as the calibration sources whose positions are known in advance.

Remark 5.

The received data can be prewhitened with a square-root of the covariance matrix of the colored noise, and then the covariance matrix of colored noise is reduced to an identity matrix accordingly, i.e., the colored noise becomes “isotropic”. It should be noted that this operation relies on a prerequisite that the noise covariance matrix is available in advance. Generally speaking, noise in the internal circuitry of the antenna elements is slowly time-varying, and the prewhitening technique works well if the colored noise is wide-sense stationary over a period of time. However, this prerequisite may not hold in some scenarios. For instance, in wireless communication systems, the cell sites equipped with antenna arrays always receive signals from users, so the noise covariance matrix cannot be obtained separately; the reverberation noise in sonar has a fast time-variance due to the underwater dynamic media, so the noise covariance matrix at the next moment varies dramatically from that at the previous moment. As a result, the prewhitening technique is invalid under these adverse circumstances while fourth-order cumulants can handle such challenges.

6. Simulation Results and Discussion

In this section, numerous simulations are offered to assess the performance of the proposed non-circular FOC method that exploits the non-circular structure via actual steering vector reparameterization. Specifically, we compare our solution with its circular counterparts, Liu’s [21] and Liao’s approaches [24] that are based on the FOC. A ULA with half-wavelength spacing between adjacent elements is considered. Similar to the settings in [39], it is assumed that the BPSK modulated incident signals are statistically independent and have identical power. The noise is assumed to be spatially-colored Gaussian and the $(m,n)$-th entry of the covariance matrix is given by $\mathbf{R}(m,n)={\sigma }_n^2{0.85}^{|m-n|}{e}^{j\frac{\pi |m-n|}{16}}$ [66,67]. The signal-to-noise ratio (SNR) is defined as $10{log}_{10}({\sigma }_s^2/{\sigma }_n^2)$. The accuracy of the DOA estimate is measured from 1000 Monte Carlo runs in terms of the root mean square error (RMSE) which is defined as

$$\begin{array}{c}\hfill \mathrm{RMSE}=\sqrt{\frac{1}{1000N}\sum _{n=1}^{1000}\sum _{i=1}^N{({\widehat{\theta }}_i^{\left(n\right)}-{\theta }_i)}^2}\end{array}$$

(31)

where ${\widehat{\theta }}_i^{\left(n\right)}$ is the estimate of ${\theta }_i$ for the n-th trial, and N is the number of signals.

We first assumed that four independent sources from $\left[-{23}^{\circ },-6^{\circ },{14}^{\circ },{36}^{\circ }\right]$ impinging on a ten-element array with mutual coupling where $P=3$ with $c_1=-0.1545+0.4755j,c_2=0.1618-0.1176j$. From Figure 1, it can be seen that with the increase of the two variables, the SNR and the number of snapshots, the RMSE of the DOA estimates declines slowly for all three methods and then stabilizes for certain SNR or snapshot values due to the manifold error resulting from the unknown mutual coupling, whereas Liao’s method has a relatively large and constant error at approximately ${1.57}^{\circ }$. It is evident that our solution is superior to Liu’s method and Liao’s method, particularly at low SNRs as well as small snapshot sizes, since a larger array aperture has been utilized for in our algorithm. Between the other two approaches, one can tell that in this scenario, Liu’s method generally outperforms Liao’s method and achieves similar performance to our method in moderate conditions (i.e., SNR $\ge -7$ dB or the number of snapshots is larger than 3500) at a cost of computation. This demonstrates that when the ULA has sufficient sensors and the degree of coupling is not high (i.e., P is small), Liu’s method performs well. In contrast, Liao’s method has a clear advantage over Liu’s for SNR lower than $-10$ dB or for less than 2000 snapshots.

In the second scenario, we further appraise the performance of the algorithm in this submission under stronger mutual coupling. More precisely, we consider the case of $P=4$ with $c_1=-0.1545+0.4755j,c_2=0.1618-0.1176j,c_3=0.0211+0.0651j$. The corresponding RMSEs versus SNR and snapshot size are illustrated in Figure 2. It can be seen that the proposed approach performs the best over the entire range of SNR values and the moderate numbers of snapshots. Nevertheless, a performance degradation appears for less than 2000 snapshots mainly because Equation (12) is guaranteed to hold under strong coupling effects only for sufficiently large samples. Contrary to the performance in the first scenario, Liu’s method is strictly inferior to Liao’s method for all SNRs. As analyzed in Section 4.1, the size of the middle sub-array adopted by Liu’s method is $M-2P+2=4$, and the size of the noise subspace is $7\times 3$, which is much smaller than the $200\times 196$ noise subspace of the proposed method and the $24\times 4$ signal subspace of Liao’s method. Due to the fact that a consistent estimate of a small subspace requires a sufficiently large number of snapshots, Liu’s method has relatively large errors for a moderate number of snapshots, even at high SNRs.

Next, we change the number of sensors to $M=8$. The other settings are kept the same as in the second scenario. As discussed in Section 4.1, Liu’s approach fails to work in this example as the number of sensors in the middle sub-array becomes $M-2P+2=2$, and up to two DOAs can be identified, but actually $N=4$ signals should be estimated herein. However, the developed algorithm in this submission and Liao’s method can still work under this more challenging scenario. The resultant RMSEs versus SNR and the number of snapshots are depicted in Figure 3. One can observe that the proposed solution significantly outperforms Liao’s method. This can be mainly attributed to the fact that Liao’s method can only utilize the signal subspace with a size of $8\times 4$, whereas the proposed method exploits the noise subspace with a larger size of $128\times 124$ and achieves a lower RMSE for the DOA estimates. However, as in the second scenario, the performance of our method deteriorates quickly and becomes inferior to Liao’s method when the number of snapshots is low. It has been shown in the literature that ESPRIT-like algorithms have clear advantages over MUSIC-like algorithms in terms of estimation errors and robustness to array manifold perturbations in the limited snapshots situation [68]. In our simulations, it was found that at least 2000 snapshots is required to sufficiently characterize the properties of FOC, especially when mutual coupling is relatively strong. It is thus reasonable that Liao’s method (ESPRIT-like) offers better performance than the proposed method (MUSIC-like) when the number of snapshots is lower than 2000.

Then, we verify the proposed scheme for the more challenging scenario where coherency between signals and mutual coupling between sensors coexist. Consider that three groups of coherent signals, three paths per group, from $\left[-{47}^{\circ },-{24}^{\circ },-9^{\circ }\right]$, $\left[-{70}^{\circ },{30}^{\circ },{60}^{\circ }\right]$, and $\left[8^{\circ },{23}^{\circ },{45}^{\circ }\right]$ impinging on a ten-element coupled ULA with $P=2$ and $c_1=-0.1545+0.4755j$. The fading amplitudes of the coherent signals were $\left[1,0.9,0.8\right]$, $\left[1,0.7,0.6\right]$, and $\left[1,0.7,0.4\right]$, while the fading phases were $\left[{48.74}^{\circ },{121.15}^{\circ },{35.66}^{\circ }\right]$, $\left[{9.35}^{\circ },{251.47}^{\circ },{103.56}^{\circ }\right]$, and $\left[{130.21}^{\circ },{16.88}^{\circ },{319.69}^{\circ }\right]$, respectively. The number of times of FBSS is $q=2$, and the SNR and number of snapshots are fixed at 15 dB and 5000. Figure 4 illustrates the DOA estimates from 50 independent trials by the developed method in such a harsh case. The solid lines mark the true DOAs. It can be seen that each DOA is correctly determined, and our solution performs satisfactorily.

In the last experiment, the computational complexity of the three algorithms versus the number of snapshots is examined. The searching step-size in the MUSIC-like methods is set as ${0.1}^{\circ }$, and the number of snapshots varies from 1000 to $\mathrm{10,000}$. Suppose that the number of sensors in the ULA $M=10$, the coupling length $P=3$, and the number of independent signals $N=2$. It can be observed from Figure 5 that our method costs the highest computation, followed by Liao’s method, while Liu’s method is most computationally efficient, which supports the analysis in Section 4.2.

7. Conclusions

This paper has addressed the problem of DOA estimation of non-circular signals with ULAs under the coexistence of unknown mutual coupling and colored noise. A new MUSIC-like approach in the context of FOC is developed to offer superior performance over the existing techniques. The proposed algorithm makes efficient use of the extended array observations and enables us to deal with DOA estimation without a priori knowledge of mutual coupling between sensors. In addition, a solution to coherent signals with mutual coupling between sensors is introduced, by the middle sub-array to suppress the electromagnetic nuisance and blind separation of coherent signals to enhance DOFs. Compared with the previous work, our solution is able to handle more signals and achieve more accurate DOA estimates even with relatively stronger mutual coupling. A series of simulation results verify the validity and efficiency of the proposed method.