Underdetermined DOA Estimation of Quasi-Stationary Signals Using a Partly-Calibrated Array

Abstract

Quasi-stationary signals have been widely found in practical applications, which have time-varying second-order statistics while staying static within local time frames. In this paper, we develop a robust direction-of-arrival (DOA) estimation algorithm for quasi-stationary signals based on the Khatri–Rao (KR) subspace approach. A partly-calibrated array is considered, in which some of the sensors have an inaccurate knowledge of the gain and phase. In detail, we first develop a closed-form solution to estimate the unknown sensor gains and phases. The array is then calibrated using the estimated sensor gains and phases which enables the improved DOA estimation. To reduce the computational complexity, we also proposed a reduced-dimensional method for DOA estimation. The exploitation of the KR subspace approach enables the proposed method to achieve a larger number of degrees-of-freedom, i.e., more sources than sensors can be estimated. The unique identification condition for the proposed method is also derived. Simulation results demonstrate the effectiveness of the proposed underdetermined DOA estimation algorithm for quasi-stationary signals.

1. Introduction

Direction-of-arrival (DOA) estimation is an important array signal processing technique which has found broad applications in radar, sonar, navigation and wireless communication [1]. In the past decades, a number of subspace-based DOA estimation algorithms have been developed including maximum likelihood (ML) [2,3,4], Capon beamforming [5,6,7], multiple signal classification (MUSIC) [8,9,10], and estimation of parameters via the rotational invariance technique (ESPRIT) [11,12]. Among them, MUSIC is the most popular one for passive localization. However, the recently developed time-reversal MUSIC algorithm can also implement active detection by utilizing the multistatic data matrix [13,14,15,16]. Moreover, [17] provides a learning-by-example approach based on the support vector machine (SVM), which exploits a multi-scaling procedure to enhance the angular resolution of the detection process in the region of incidence of the incoming waves. In addition, the sparse signal reconstruction-based DOA estimation algorithm [18] also shows good performance even with a single snapshot.

The methods mentioned above provide high accuracy DOA estimation when the array is well calibrated. However, their performances degrade when some of the array sensors are mis-calibrated. For instance, an array may suffer imperfect calibration [19] and/or mutual coupling [20] which may vary over time and environment, yielding perturbations in the sensor gain and phase characteristics. In this paper, we mainly focus on the calibration problem. Generally, the calibration methods are classified into three categories, i.e., pilot calibration, self-calibration and auto-calibration. The pilot calibration methods [21,22,23] utilize sources with known parameters, such as direction and location, to estimate the array uncertainties analytically by exploiting the mathematical model of the array response [24]. These kinds of algorithms have highly accurate estimation performance. However, they also highly depend on the accurate knowledge of the pilot sources, which may not always be available in practice. In addition, the pilot calibration may be expensive and time-consuming because of the off-line operation. The self-calibration methods [19,25] have aroused much interest because they do not require the pilot sources. The self-calibration methods can estimate the gain-phase uncertainties and DOA simultaneously, which makes them more attractive. The auto-calibration methods [24,26] do not use any external source for calibration. Instead, the array elements are transceivers transmitting to each other and these signals are used to remove the array uncertainties.

In practice, the response of some sensor elements is poorly known or even unknown when the array sensors are incompletely calibrated. This kind of array is known as partly-calibrated array [27,28,29]. A number of calibration algorithms have been proposed to ensure robust DOA estimation using partly-calibrated arrays [19,27,28,29]. Among them, an iterative method is proposed in [19] to perform array calibration and DOA estimation simultaneously. This approach, however, is suboptimal and works properly only when the array perturbations are small. An iterative algorithm is developed in [27] for DOA and sensor gain-phase joint estimation, where a good initialization is required to achieve the convergence to the global minimum. In [28], an ESPRIT-like algorithm is developed to achieve the joint estimation of DOA and unknown gains and phases without any iteration. In [29], a subspace-based DOA estimation approach is proposed, which is applicable to the partly-calibrated arrays composed of multiple well-calibrated subarrays of arbitrary known geometry. Moreover, the displacements of the sensor elements also lead to DOA estimation performance degradation. Utilizing the stochastic collocation method, a method is proposed in [30] to calculate the cumulative density functions (cdfs) of the DOA estimates due to a random displacement of one sensor element in an array. However, it does not achieve a robust DOA estimation result.

Limited by the degrees-of-freedom (DOFs), the number of sources is usually assumed to be less than the number of sensors, as in [27,28,29]. Hence, increasing the DOFs has attracted continuing interest. In the past few years, several algorithms have been developed to deal with the underdetermined DOA estimation problem, i.e., the number of sources is larger than the number of sensors, either by exploiting the high-order statistics [31] or the sparse array structures [32,33,34,35]. In [36], a Khatri–Rao (KR) subspace approach is proposed by exploiting the statistical properties of the quasi-stationary signals, which represent a class of commonly encountered nonstationary signals including speech and audio signals. The statistics of quasi-stationary signals are locally static over a short time period, but exhibit difference from one local time frame to another [36,37]. By exploiting such statistical properties, the KR subspace approach achieves $2(N-1)$ DOFs with an N-element uniform linear array (ULA).

In this paper, we propose a DOA estimation algorithm based on the KR subspace approach using the partly-calibrated array. The proposed method is a kind of self-calibration method. Due to the existence of the unknown sensor gain and phase perturbations, the conventional KR-MUSIC and KR-ESPRIT methods [36] generally suffer significant performance degradation. To address this issue, we develop a KR subspace-based robust DOA estimation algorithm by simultaneously estimating the unknown gain and phase uncertainties. By utilizing an ESPRIT-like algorithm, the DOA as well as the unknown gain and phase uncertainties are jointly estimated with closed-form expressions. As such, unlike other joint DOA estimation and calibration methods requiring iterative solution or spectral search, the proposed method can be effectively implemented to achieve a high accuracy. Different to the method in [28], the proposed method achieves higher DOFs by exploiting the KR subspace. After the array is fully calibrated with the estimated gain and phase errors, it enables the underdetermined DOA estimation problem to be solved with more sources than sensors. In order to reduce the computational complexity, we further propose a reduced-dimensional (RD) method by designing a RD matrix. Simulation results demonstrate the effectiveness of the proposed DOA estimation algorithm in the presence of several sensor gain and phase errors even with more sources than sensors.

The remainder of this paper is organized as follows. In Section 2, we first describe the signal model with a partly-calibrated array. Using the quasi-stationary signals, the KR subspace is then presented. In Section 3, we propose an underdetermined DOA estimation algorithm with a partly-calibrated array for quasi-stationary signals. Furthermore, a RD method is also provided. In Section 4, we discuss the condition for unique identification and computational complexity of the proposed method. Simulation results are given in Section 5 and Section 6 concludes this paper.

Notations used in this paper are as follows. Lower-case (upper-case) bold characters are used to denote vectors (matrices). ${\left(·\right)}^T$, ${\left(·\right)}^{*}$ and ${\left(·\right)}^H$ denote the transpose, conjugate and Hermitian operations of a matrix or vector, respectively. $\mathrm{E}[·]$ stands for the expectation operator. ${∥·∥}_F$ denotes the Frobenius norm. $\mathrm{trace}[\mathbf{A}]$ and $\mathrm{rank}[\mathbf{A}]$ denote the trace and rank of matrix $\mathbf{A}$, respectively. ${\mathbb{C}}^{M\times N}$ denotes an $M\times N$ complex matrix or vector (when $N=1$). $\mathrm{vec}(\mathbf{x})$ is a vectorization operator that turns a matrix $\mathbf{x}$ into a vector by concatenating all columns. In addition, $\mathrm{diag}(\mathbf{x})$ denotes a diagonal matrix with the elements of $\mathbf{x}$ constituting the diagonal entries. ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix, ${\mathbf{1}}_N$ and ${\mathbf{0}}_N$ stand for the $N\times 1$ all-one and all-zero vector, respectively. ⊗, ⊙ and ∘ denote the Kronecker product, Khatri–Rao product and Schur–Hadamart product, respectively. $\angle (x)$ stands for the phase of variable x.

2. Signal Model

In this section, we first describe the partly-calibrated array signal model, then present the KR subspace utilizing the quasi-stationary signals.

2.1. Partly-Calibrated Array Model

Consider a partly-calibrated ULA with N omnidirectional antennas. We assume that there are K uncorrelated far-field sources impinging on the array with distinct DOAs $\mathit{\theta }={[{\theta }_1,\cdots ,{\theta }_K]}^T$. For the DOA estimation problem, it is a common scenario to use the first sensor as the reference. Thus, without loss of generality, it is practical to assume that the first $N_c$ sensors are well-calibrated, whereas the last $N-N_c$ sensors are uncalibrated with unknown, direction-independent gain and phase uncertainties [27,28]. It is worth noting that the algorithm remains effective when the calibrated sensors are located in arbitrary positions, provided that there is at least one pair of consecutive calibrated sensors [28]. The complex-valued baseband array received signal vector at time index t can be modeled as

$$\mathbf{x}(t)=\overline{\mathbf{A}}\mathbf{s}(t)+\mathbf{n}(t),$$

(1)

where $\overline{\mathbf{A}}=\mathsf{\Gamma }(\mathit{\gamma })\mathbf{A}$ is the actual array manifold matrix, $\mathsf{\Gamma }(\mathit{\gamma })=\mathrm{diag}(\mathit{\gamma })$ with $\mathit{\gamma }={[{\mathbf{1}}_{N_c}^T,{\rho }_1{e}^{j{\phi }_1},\cdots ,{\rho }_{N-N_c}{e}^{j{\phi }_{N-N_c}}]}^T$ corresponding to the unknown sensor gain-phase vector. It is assumed that the sensor gains are nonzero, i.e., ${\rho }_n\ne 0,\forall n=1,\cdots ,N-N_c$. In addition, $\mathbf{A}=[\mathbf{a}({\theta }_1),\cdots ,\mathbf{a}({\theta }_K)]$ is the assumed array manifold matrix with $\mathbf{a}({\theta }_k)={[1,e^{-j2\pi d/\lambda sin{\theta }_k},\cdots ,e^{-j2\pi d/\lambda (N-1)sin{\theta }_k}]}^T\in {\mathbb{C}}^N$, where d and $\lambda$ denote the inter-element spacing and the signal wavelength, respectively. $\mathbf{s}(t)={[s_1(t),\cdots ,s_K(t)]}^T\in {\mathbb{C}}^K$ denotes the zero-mean quasi-stationary signal waveforms, and $\mathbf{n}(t)\sim \mathcal{CN}(\mathbf{0},{\sigma }_n^2\mathbf{I})\in {\mathbb{C}}^N$ is the additive white Gaussian noise vector with zero-mean and covariance matrix ${\sigma }_n^2{\mathbf{I}}_N$. We further assume that the source signals $s_k(t),k=1,\cdots ,K$ are independent of each other. In addition, the source signals and noise are also assumed to be uncorrelated.

2.2. Khatri–Rao Product Subspace

Quasi-stationary signals have a wide distribution in nature, such as speech and audio signals. Their second-order statistics remain static at a local time frame, however, they vary with different time frames. As shown in Figure 1, we illustrate the waveform of a quasi-stationary signal. Similar to [36], we model a quasi-stationary source signal $s_k(t)$ to have M non-overlapped time frames and the length of each frame is assumed to be L. In this case,

$$\mathrm{E}\left[{\left|s_k(t)\right|}^2\right]={\sigma }_{mk}^2,\forall t\in \left[\left(m-1\right)L,mL-1\right],m=1,\cdots M,$$

(2)

where ${\sigma }_{mk}^2$ is the signal power of $s_k(t)$ at the mth time frame. For the mth time frame, the sample covariance matrix is written as

$${\mathbf{R}}_m=\mathrm{E}\left[\mathbf{x}(t){\mathbf{x}}^H(t)\right]=\overline{\mathbf{A}}{\mathbf{D}}_m{\overline{\mathbf{A}}}^H+{\sigma }_n^2{\mathbf{I}}_N,$$

(3)

where ${\mathbf{D}}_m=\mathsf{\Gamma }\left([{\sigma }_{m1}^2,\cdots ,{\sigma }_{mK}^2]\right)$. In practice, due to the finite number of snapshots, the covariance matrix is evaluated as ${\widehat{\mathbf{R}}}_m=(1/L){\sum }_{t=(m-1)L}^{mL-1}\mathbf{x}(t){\mathbf{x}}^H(t)$.

We vectorize ${\widehat{\mathbf{R}}}_m$ to an $N^2\times 1$ vector as:

$${\mathbf{y}}_m=\mathrm{vec}({\widehat{\mathbf{R}}}_m)=\breve{\mathbf{A}}{\mathbf{d}}_m+{\sigma }_n^2\mathbf{i},$$

(4)

where $\breve{\mathbf{A}}={\overline{\mathbf{A}}}^{*}\odot \overline{\mathbf{A}}$ is the generalized manifold matrix, ${\mathbf{d}}_m={[{\sigma }_{m1}^2,\cdots ,{\sigma }_{mK}^2]}^T$ and $\mathbf{i}=\mathrm{vec}({\mathbf{I}}_N)$. We can see that the expression of (4) is similar to the signal model in (1). To exploit the second-order properties of the quasi-stationary signals, we stack ${\mathbf{y}}_m$ over all M time frames to yield the following matrix

$$\mathbf{Y}=[{\mathbf{y}}_1,\cdots ,{\mathbf{y}}_M]=\breve{\mathbf{A}}\mathsf{\Psi }+{\sigma }_n^2\mathbf{i}{\mathbf{1}}_M^T,$$

(5)

where $\mathsf{\Psi }=[{\mathbf{d}}_1,\cdots ,{\mathbf{d}}_M]$. Thus, the DOF of the array is extended.

As shown in [36], the noise covariance term ${\sigma }_n^2\mathbf{i}{\mathbf{1}}_M^T$ can be eliminated by utilizing the orthogonal complement matrix ${\mathbf{P}}_{{\mathbf{1}}_M}^{\perp }={\mathbf{I}}_M-(1/M){\mathbf{1}}_M{\mathbf{1}}_M^T$, as

$$\tilde{\mathbf{Y}}=\mathbf{Y}{\mathbf{P}}_{{\mathbf{1}}_M}^{\perp }=\breve{\mathbf{A}}\mathsf{\Psi }{\mathbf{P}}_{{\mathbf{1}}_M}^{\perp },$$

(6)

the singular value decomposition (SVD) of which is expressed as

$$\tilde{\mathbf{Y}}=\left[{\mathbf{E}}_s{\mathbf{E}}_n\right]\left[\begin{array}{cc}{\mathsf{\Sigma }}_s& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_s^H\\ {\mathbf{V}}_n^H\end{array}\right],$$

(7)

where ${\mathbf{E}}_s$ and ${\mathbf{E}}_n$ are the left singular matrices corresponding to the nonzero and zero singular values, respectively, ${\mathbf{V}}_s$ and ${\mathbf{V}}_n$ are the respective right singular matrices, and ${\mathsf{\Sigma }}_s$ is a diagonal matrix whose diagonal entries contain the nonzero singular values. Based on the subspace theory, matrix ${\mathbf{E}}_s$ spans the same subspace as the generalized manifold matrix $\breve{\mathbf{A}}$ which is denoted as the signal subspace. As such, these two matrices can be related by

$${\mathbf{E}}_s=\breve{\mathbf{A}}\mathbf{T},$$

(8)

where $\mathbf{T}$ is a $K\times K$ nonsingular matrix.

3. The Proposed Method

Due to the unknown sensor gain and phase uncertainties, the conventional KR subspace approach-based DOA estimation algorithms [36,37,38] suffer performance degradation or even failed operation. In this section, we first estimate the unknown gains and phases, then, the DOAs of the sources are estimated utilizing the ESPRIT-like algorithm based on the KR subspace approach.

3.1. Joint Parameters Estimation

The generalized manifold matrix $\breve{\mathbf{A}}$ in (4) can be rewritten as

$$\begin{array}{cc}\hfill \breve{\mathbf{A}}& ={\left(\mathsf{\Gamma }(\mathit{\gamma })\mathbf{A}\right)}^{*}\odot \left(\mathsf{\Gamma }(\mathit{\gamma })\mathbf{A}\right)\hfill \\ & =\mathsf{\Gamma }({\mathit{\gamma }}^{*}\otimes \mathit{\gamma })({\mathbf{A}}^{*}\odot \mathbf{A})=\mathsf{\Gamma }(\tilde{\mathit{\gamma }})\tilde{\mathbf{A}},\hfill \end{array}$$

(9)

with $\tilde{\mathit{\gamma }}={\mathit{\gamma }}^{*}\otimes \mathit{\gamma }$ and $\tilde{\mathbf{A}}={\mathbf{A}}^{*}\odot \mathbf{A}$.

Similar to the conventional ESPRIT algorithm, the partly-calibrated ULA array is divided into two overlapping subarrays to exploit the rotational invariance property. Denote ${\mathit{\gamma }}_1$ and ${\mathit{\gamma }}_2$ as vectors containing the first and last $N-1$ elements of $\mathit{\gamma }$, and ${\mathbf{A}}_1$ and ${\mathbf{A}}_2$ matrices with the first and last $N-1$ rows of $\mathbf{A}$, respectively. In addition, we define ${\tilde{\mathit{\gamma }}}_1={\mathit{\gamma }}_1^{*}\otimes \mathit{\gamma }$, ${\tilde{\mathit{\gamma }}}_2={\mathit{\gamma }}_2^{*}\otimes \mathit{\gamma }$, ${\tilde{\mathbf{A}}}_1={\mathbf{A}}_1^{*}\odot \mathbf{A}$ and ${\tilde{\mathbf{A}}}_2={\mathbf{A}}_2^{*}\odot \mathbf{A}$. By selecting the first and last $N(N-1)$ rows of ${\mathbf{E}}_s$, we have the two submatrices ${\mathbf{E}}_{s1}$ and ${\mathbf{E}}_{s2}$ as

$$\begin{array}{c}{\mathbf{E}}_{s1}=\mathsf{\Gamma }({\tilde{\mathit{\gamma }}}_1){\tilde{\mathbf{A}}}_1\mathbf{T},\hfill \\ {\mathbf{E}}_{s2}=\mathsf{\Gamma }({\tilde{\mathit{\gamma }}}_2){\tilde{\mathbf{A}}}_2\mathbf{T}.\hfill \end{array}$$

(10)

Using the rotational invariance property ${\tilde{\mathbf{A}}}_2={\tilde{\mathbf{A}}}_1\mathsf{\Phi }$ with $\mathsf{\Phi }=\mathsf{\Gamma }\left({\left[e^{-j2\pi d/\lambda sin{\theta }_1},\cdots ,e^{-j2\pi d/\lambda sin{\theta }_K}\right]}^T\right)$, we can further obtain the following equation from the result of (10)

$$\mathsf{\Gamma }(\mathit{Y}){\mathbf{E}}_{s2}={\mathbf{E}}_{s1}\mathsf{\Pi },$$

(11)

where $\mathsf{\Pi }={\mathbf{T}}^{-1}\mathsf{\Phi }\mathbf{T}$, $\mathsf{\Gamma }(\mathit{Y})=\mathsf{\Gamma }({\tilde{\mathit{\gamma }}}_1)\mathsf{\Gamma }({\tilde{\mathit{\gamma }}}_2^{-1})=\mathsf{\Gamma }({\tilde{\mathit{\gamma }}}_1\circ {\tilde{\mathit{\gamma }}}_2^{-1})$ with

$$\begin{array}{cc}\hfill \mathit{Y}={\tilde{\mathit{\gamma }}}_1\circ {\tilde{\mathit{\gamma }}}_2^{-1}& =({\mathit{\gamma }}_1^{*}\otimes \mathit{\gamma })\circ {({\mathit{\gamma }}_2^{*}\otimes \mathit{\gamma })}^{-1}\hfill \\ & =\left({\mathit{\gamma }}_1^{*}\circ {({\mathit{\gamma }}_2^{*})}^{-1}\right)\otimes {\mathbf{1}}_N.\hfill \end{array}$$

(12)

More specifically, $\mathit{Y}$ can be expressed as

$$\begin{array}{c}\mathit{Y}=[{\mathbf{1}}_{(N_c-1)}^T,{\left({\rho }_1{e}^{j{\phi }_1}\right)}^{-1},{\rho }_1{e}^{j{\phi }_1}{({\rho }_2{e}^{j{\phi }_2})}^{-1},\cdots ,{\rho }_{N-N_c-1}{e}^{j{\phi }_{N-N_c-1}}{({\rho }_{N-N_c}{e}^{j{\phi }_{N-N_c}})}^{-1}{]}^H\otimes {\mathbf{1}}_N.\end{array}$$

(13)

In the underlying problem with uncalibrated sensors, both $\mathsf{\Gamma }(\mathit{Y})$ and $\mathsf{\Pi }$ in (11) are unknown. As such, the conventional ESPRIT algorithm [39] cannot be applied directly. Here, $\mathsf{\Gamma }(\mathit{Y})$ and $\mathsf{\Pi }$ can be estimated by solving the following optimization problem

$$\begin{array}{cc}\hfill \underset{\mathit{Y},\mathsf{\Pi }}{min}& {∥\mathsf{\Gamma }(\mathit{Y}){\mathbf{E}}_{s2}-{\mathbf{E}}_{s1}\mathsf{\Pi }∥}_F^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{W}\mathit{Y}={\mathbf{1}}_{(N_c-1)N},\hfill \end{array}$$

(14)

where $\mathbf{W}$ is defined as

$$\mathbf{W}=\left[\begin{array}{cc}{\mathbf{I}}_{(N_c-1)N}& {\mathbf{0}}_{(N_c-1)N\times (N(N-1)-(N_c-1)N)}\end{array}\right].$$

The equality constraint in (14) is introduced to ensure that the first $N_c$ sensors take the calibrated gain and phase values. When the calibrated sensors are in arbitrary positions, the matrix $\mathbf{W}$ should be re-arranged to ensure that the corresponding sensors are well-calibrated.

With some simplifications as shown in the Appendix A, the original optimization problem (14) can be further reformulated as

$$\begin{array}{cc}\hfill \underset{\mathit{Y}}{min}& {\mathit{Y}}^H\left({({\mathbf{E}}_{s2}{\mathbf{E}}_{s2}^H)}^T\circ \mathbf{P}\right)\mathit{Y}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{W}\mathit{Y}={\mathbf{1}}_{(N_c-1)N},\hfill \end{array}$$

(15)

which can be effectively solved using the Lagrange multipliers method. The matrix $\mathbf{P}$ is expressed as $\mathbf{P}={\mathbf{I}}_{N(N-1)}-{\mathbf{E}}_{s1}{({\mathbf{E}}_{s1}^H{\mathbf{E}}_{s1})}^{-1}{\mathbf{E}}_{s1}^H$. The detailed derivation of the solution is given in Appendix B. The solution of the optimization problem (15) is expressed as

$$\widehat{\mathit{Y}}={\mathbf{Q}}^{-1}{\mathbf{W}}^T{(\mathbf{W}{\mathbf{Q}}^{-1}{\mathbf{W}}^T)}^{-1}{\mathbf{1}}_{(N_c-1)N},$$

(16)

where $\mathbf{Q}={({\mathbf{E}}_{s2}{\mathbf{E}}_{s2}^H)}^T\circ \mathbf{P}$. Note that $\mathbf{Q}$ is required to be nonsingular to estimate the $\widehat{\mathit{Y}}$. In the case of infinite snapshots, $\mathbf{Q}$ is a nonsingular matrix. To ensure that $\mathbf{Q}$ is also a nonsingular matrix with finite snapshots in practice, diagonal loading [40,41] is a possible method to handle this problem. In addition, from the extensive experiments with finite snapshots that we have made, the matrix $\mathbf{Q}$ is always nonsingular. Thus, it is not necessary to use the diagonal loading method for $\mathbf{Q}$ in general [28].

Combining with (13), the ith entry of $\widehat{\mathit{\gamma }}$ can be estimated as

$${\widehat{\mathit{\gamma }}}_i={\left(\prod _{n=N_c}^{i-1}\mathcal{T}(n)\right)}^{-1},i=N_c+1,\cdots ,N,$$

(17)

where $\mathcal{T}(n)=\frac{1}{N}{\sum }_{l=(n-1)N+1}^{nN}{\widehat{\mathit{Y}}}_l^{*}$, and ${\widehat{\mathit{Y}}}_l$ is the lth element of $\widehat{\mathit{Y}}$.

By substituting the estimated $\widehat{\breve{\mathit{\gamma }}}$ in (A8) into (A1), $\mathsf{\Pi }$ can be estimated as $\widehat{\mathsf{\Pi }}$. Since $\mathsf{\Pi }$ and $\mathsf{\Phi }$ are similar matrices, the eigenvalues of $\mathsf{\Pi }$ are the diagonal entries of $\mathsf{\Phi }$ and the columns of $\mathbf{T}$ are the eigenvectors of $\mathsf{\Pi }$. Thus, the kth signal DOA can be estimated by

$${\widehat{\theta }}_k=arcsin\left(-\frac{\lambda }{2\pi d}\angle ({\widehat{\upsilon }}_k)\right),k=1,\cdots ,K,$$

(18)

where ${\widehat{\upsilon }}_k$ is the kth eigenvalue of $\widehat{\mathsf{\Pi }}$.

3.2. The Proposed RD Method

Once the unknown sensor gain-phase vector $\widehat{\mathit{\gamma }}$ is estimated, the sensor array can be calibrated and the corresponding signal subspace can be expressed as

$${\widehat{\mathbf{E}}}_s=\mathsf{\Gamma }({\widehat{\tilde{\mathit{\gamma }}}}^{-1}){\mathbf{E}}_s=\tilde{\mathbf{A}}\mathbf{T},$$

(19)

where $\widehat{\tilde{\mathit{\gamma }}}={({\widehat{\mathit{\gamma }}}^{-1})}^{*}\otimes \widehat{\mathit{\gamma }}$.

Owing to the fact that there only exist $2N-1$ distinct rows in $\tilde{\mathbf{A}}$, the effective dimension of $\tilde{\mathbf{A}}$ can be reduced without any performance loss [36]. Let $\tilde{\mathbf{A}}=\mathbf{G}\mathbf{B}$, where $\mathbf{G}$ is an $N^2\times (2N-1)$ RD matrix and $\mathbf{B}=[{\mathbf{b}}_1,\cdots ,{\mathbf{b}}_K]$ with ${\mathbf{b}}_k=[e^{j2\pi (d/\lambda )(N-1)sin{\theta }_k},\cdots ,e^{j2\pi (d/\lambda )sin{\theta }_k},1,e^{-j2\pi (d/\lambda )sin{\theta }_k},\cdots ,$ $e^{-j2\pi (d/\lambda )(N-1)sin{\theta }_k}]\in {\mathbb{C}}^{2N-1}$.

With ${\widehat{\mathbf{E}}}_s=\tilde{\mathbf{A}}\mathbf{T}$ and $\tilde{\mathbf{A}}=\mathbf{G}\mathbf{B}$, the RD signal subspace can be expressed as

$${\widehat{\mathbf{E}}}_{ds}=\mathsf{\Theta }{\mathbf{G}}^T{\widehat{\mathbf{E}}}_s=\mathsf{\Theta }{\mathbf{G}}^T\tilde{\mathbf{A}}\mathbf{T}=\mathbf{B}\mathbf{T},$$

(20)

where $\mathsf{\Theta }={({\mathbf{G}}^T\mathbf{G})}^{-1}$.

To apply the ESPRIT algorithm, we construct two overlapping matrices by selecting the first $2N-2$ and last $2N-2$ rows of ${\widehat{\mathbf{E}}}_{ds}$, which are denoted as ${\widehat{\mathbf{E}}}_{ds1}$ and ${\widehat{\mathbf{E}}}_{ds2}$, respectively. We have

$$\begin{array}{c}{\widehat{\mathbf{E}}}_{ds1}={\mathbf{B}}_1\mathbf{T}\hfill \\ {\widehat{\mathbf{E}}}_{ds2}={\mathbf{B}}_2\mathbf{T},\hfill \end{array}$$

(21)

where ${\mathbf{B}}_1$ and ${\mathbf{B}}_2$ denote the first $2N-2$ rows and last $2N-2$ rows of $\mathbf{B}$, respectively. Obviously, ${\mathbf{B}}_1$ and ${\mathbf{B}}_2$ hold the rotational invariance property, i.e., ${\mathbf{B}}_2={\mathbf{B}}_1\mathsf{\Phi }$. Thus, the DOAs can be estimated by utilizing the ESPRIT algorithm.

Although the proposed method is implemented using a ULA, the proposed method is applicable on other kinds of arrays which hold rotational invariance property, such as uniformly circle arrays (UCAs).

4. Condition for Unique Identification and Computational Complexity

First, we discuss the condition under which the true DOA can be uniquely estimated, which is equivalent to finding the condition for ${\widehat{\mathbf{E}}}_{ds1}$ and ${\widehat{\mathbf{E}}}_{ds2}$ to have the full column rank.

Since $\tilde{\mathit{\gamma }}$ is an $N^2\times 1$ vector with nonzero entries, we have $\mathrm{rank}[\mathsf{\Gamma }(\tilde{\mathit{\gamma }})]=N^2$. Also note that $\breve{\mathbf{A}}=\mathsf{\Gamma }(\tilde{\mathit{\gamma }})\tilde{\mathbf{A}}=\mathsf{\Gamma }(\tilde{\mathit{\gamma }})\mathbf{G}\mathbf{B}$, $\mathrm{rank}[\mathbf{G}]=2N-1$ and $\mathrm{rank}[\mathbf{B}]=min\{2N-1,K\}$, the rank of $\breve{\mathbf{A}}$ is given as

$$\mathrm{rank}[\breve{\mathbf{A}}]=min\left\{\mathrm{rank}\left[\mathsf{\Gamma }(\tilde{\mathit{\gamma }})\right],\mathrm{rank}[\mathbf{G}],\mathrm{rank}[\mathbf{B}]\right\}=min\left\{N^2,2N-1,K\right\}.$$

(22)

For a multi-sensor array, $N>1$, thus $N^2>2N-1$, and $\mathrm{rank}[\breve{\mathbf{A}}]=min\{2N-1,K\}$. Then, the full column rank condition for $\breve{\mathbf{A}}$ is guaranteed as long as $K\le 2N-1$ is satisfied.

Since the DOAs are estimated by exploiting the ESPRIT algorithm, the KR subspaces ${\tilde{\mathbf{E}}}_{ds1}$ and ${\tilde{\mathbf{E}}}_{ds2}$ will lose one DOF. As a result, the proposed method can resolve up to $K\le 2(N-1)$ sources. Compared with the standard ESPRIT method which only resolves up to $N-1$ sources, the proposed method is very attractive when there are more sources than sensors.

Next, we discuss the computational complexity of the proposed method. The SVD operation requires $\mathcal{O}\{N^2{M}^2\}$. To estimate the unknown gain-phase error using (A8), $\mathcal{O}\{N^3{(N-1)}^2(N_c-1)+N^3(N-1)(N_c-1)(N+N_c-2)+N^3(N_c-1)\}$ is required. In this paper, the ESPRIT algorithm is applied to perform DOA estimation, the main computational load is the eigenvalue decomposition (EVD) process which requires $\mathcal{O}\{N^3{(N-1)}^3\}$. Utilizing the proposed RD method, the dimension of the signal subspace is reduced from $N(N-1)\times K$ to $(2N-2)\times K$. Thus, the computational load for EVD operation becomes $\mathcal{O}\{{(2N-2)}^3\}$. We summarize the computational complexities of the two methods in

Table 1. Comparison of the computational complexity of two methods.

Method	SVD Operation	Gain-Phase Error Estimation	EVD for DOA Estimation
Proposed method	$\mathcal{O}\left\{N^2{M}^2\right\}$	$\mathcal{O}\left\{N^3{(N-1)}^2(N_c-1)+N^3(N-1)(N_c-1)(N+N_c-2)+N^3(N_c-1)\right\}$	$\mathcal{O}\left\{N^3{(N-1)}^3\right\}$
Proposed RD method			$\mathcal{O}\left\{{(2N-2)}^3\right\}$

. It is clear that the computational complexity is significantly reduced by using the proposed RD method.

5. Simulation

In this section, simulation results are provided to demonstrate the effectiveness of the proposed DOA estimation algorithm based on the KR subspace approach using a partly-calibrated array.

In the first simulation, we consider a 10-element ($N=10$) partly-calibrated ULA with half-wavelength interelement spacing. The first five sensors are assumed to be well-calibrated, i.e., $N_c=5$, while the last five sensors are uncalibrated with direction-independent gain and phase uncertainties, $\mathit{\gamma }=[{\mathbf{1}}_5^T,0.8e^{j\pi /5},1.2e^{j\pi /5},1.53e^{-j\pi /5},0.75e^{j\pi /4},$ $1.36e^{-j\pi /3}{]}^T$. Assume that there exist four far-field sources with DOAs $-{25}^{\circ },{10}^{\circ },{20}^{\circ }$ and ${30}^{\circ }$, respectively. The frame period of each source is set as $L=1024$ and the number of frames is $M=50$. The signal-to-noise ratio (SNR) is defined as [36]

$$\mathrm{SNR}=\frac{\frac{1}{T}{\sum }_{t=1}^T\mathrm{E}\left[\parallel \overline{\mathbf{A}}{\mathbf{s}(t)\parallel }^2\right]}{\mathrm{E}\left[{\parallel \mathbf{n}(t)\parallel }^2\right]},$$

(23)

where $T=LM$. In the following simulation, we use 100 Monte Carlo trials to evaluate the performance of the proposed method.

In

Table 2. Performance of gain and phase estimation.

(a) Gain Estimation Results				(b) Phase Estimation Results (Radian)
Index	True Value	Mean	STD	Index	True Value	Mean	STD
${\widehat{\rho }}_1$	0.8000	0.8056	0.0049	${\widehat{\phi }}_1$	0.6283	0.6286	0.0055
${\widehat{\rho }}_2$	1.2500	1.2510	0.0085	${\widehat{\phi }}_2$	−1.0472	−1.0473	0.0061
${\widehat{\rho }}_3$	1.5300	1.5344	0.0118	${\widehat{\phi }}_3$	−0.6283	−0.6284	0.0074
${\widehat{\rho }}_4$	0.7500	0.7608	0.0077	${\widehat{\phi }}_4$	0.7854	0.7855	0.0092
${\widehat{\rho }}_5$	1.3600	1.3700	0.0164	${\widehat{\phi }}_5$	−1.0472	−1.0481	0.0099

, we present the performance of the proposed gain-phase uncertainties estimation algorithm. The mean value and standard deviation results are calculated at a fixed SNR of 10 dB. It is evident that the proposed gain-phase estimation method achieves accurate gain and phase estimation.

In Figure 2, the root mean square error (RMSE) of the proposed ESPRIT-like algorithm is compared with that of the KR-ESPRIT algorithm [36] at different SNRs. Here, the RMSE of DOA estimation is defined as

$$\mathrm{RMSE}=\sqrt{\mathrm{E}\left[\frac{1}{K}\sum _{k=1}^K{|{\widehat{\theta }}_k-{\theta }_k|}^2\right]},$$

(24)

where ${\widehat{\theta }}_k$ is the DOA estimation of the kth source. As shown in Figure 2, the KR-ESPRIT algorithm suffers significant performance degradation, while the proposed method shows accurate DOA estimation results. It also verifies that the proposed RD ESPRIT algorithm achieves almost the same DOA estimation without performance loss.

To evaluate the effect of frame length on estimation performance, we compare the RMSEs with different frame period lengths. In this simulation, the proposed RD method is used. Without loss of generality, the total number of snapshots used in this simulation is fixed as $T=ML=51,200$, and M varies with 256, 512 and 1024. From Figure 3, we can claim that with the increase of M, the performance degrades in the low SNR cases. In the high SNR cases, almost the same performance can be achieved with different frame periods. Note that the estimation of KR signal subspace ${\mathbf{E}}_s$ requires a sufficient number of frames. When L becomes small, the estimation of ${\mathbf{E}}_s$ will be distorted, which leads to a degraded estimation performance in the low SNR cases.

In the next simulation, we consider the underdetermined DOA estimation case. Consider a 5-element partly-calibrated ULA with $\mathit{\gamma }={[1,1,1.2e^{j\pi /4},1,0.86e^{-j\pi /6}]}^T$. Assume that there exist six sources ($K=6$) with DOAs $-{35}^{\circ },-{25}^{\circ },-{15}^{\circ },0^{\circ },{20}^{\circ }$, and ${30}^{\circ }$, respectively. Here, the SNR is set as 10 dB. In addition, the frame length is selected as $L=1024$, and the total number of snapshots is set as $T=51,200$.

The performance comparison using the proposed algorithm method and the proposed RD method is plotted in Figure 4. It demonstrates that the proposed methods work well even when the number of sources is greater than the number of sensors. In addition, the proposed RD method achieves almost the same estimation performance with reduced complexity.

The gain and phase uncertainties estimation with the underdetermined case is presented in

Table 3. Performance of gain and phase estimation with the underdetermined case.

(a) Gain Estimation Results				(b) Phase Estimation Results (Radian)
Index	True Value	Mean	STD	Index	True Value	Mean	STD
${\widehat{\rho }}_1$	1.2000	1.1990	0.0138	${\widehat{\phi }}_1$	0.7853	0.7836	0.0093
${\widehat{\rho }}_2$	0.8600	0.8749	0.0267	${\widehat{\phi }}_2$	−0.5236	−0.5269	0.0123

. The SNR used in this experiment is 10 dB. From Table 3, we can see that the gain-phase uncertainties estimation still has a solid performance.

The performance of the proposed DOA estimation method with the underdetermined case is evaluated by calculating the RMSEs versus SNRs with a different number of targets. As is shown in Figure 5, with the increased number of sources, the performance is heavily degraded in the low SNR case. In the high SNR case, the RMSE curves converge and the proposed method achieves good performance.

6. Conclusions

In this paper, we proposed a DOA estimation algorithm for quasi-stationary signals based on the KR subspace approach with a partly-calibrated array. We first developed a closed-form expression to estimate the unknown sensor gains and phases in the KR subspace. Then, the signal DOAs were estimated by applying the ESPRIT algorithm. In order to reduce the computational complexity, a RD method was developed. Finally, we derived the unique identification condition for the proposed method and the DOFs of the proposed method were also analyzed. We have analytically shown and verified through simulations that the proposed RD method can estimate the DOAs of more sources than sensors with a partly-calibrated array. In the future, it will be interesting and practical to study the correlated/coherent sources scenario based on the KR subspace method with a partly-calibrated array.