A Novel Direction-of-Arrival Estimation via Phase Retrieval with Unknown Sensor Gain-and-Phase Errors

In signal array processing, high-resolution direction-of-arrival (DOA) estimation algorithms work well on the assumption that the system models are perfect. However, in practicality, there are imperfect system models in which sensor gain-and-phase errors are considered. In this paper, we propose a novel framework that can effectively solve direction-of-arrival estimation tasks in the presence of sensor gain-and-phase errors. In contrast to existing approaches based on phase retrieval, our method eliminates gain errors by using the compensated covariance matrix. Meanwhile, we propose a data preprocessing method by taking only one column of the compensated covariance matrix without losing any magnitude information. Additionally, the phase retrieval problem is formed by the proposed data preprocessing method. Furthermore, the phase retrieval problem is solved by the recently proposed sparse feasible point pursuit algorithm, and DOA estimates are obtained. To prevent the model from ambiguities, we employ the known DOA to place reference sources. Numerical results show that the proposed scheme achieves better performance compared to state-of-the-art approaches.


Introduction
Direction-of-arrival (DOA) estimation of plane waves using sensor arrays is a key problem which is frequently encountered in many applications, such as radar, sonar, biomedical engineering, and astronomy. Various high-resolution algorithms for DOA estimation, such as Multiple Signals Classification (MUSIC) [1], Estimation of Signal Parameters via Rotational Invariance Technique (ESPRIT) [2], and the Propagator method [3] have been proposed. However, this class of approaches is based on the assumption of the ideal array manifold. In practical applications, there always exist some array perturbations, such as sensor gain-and-phase errors, which will cause severe degradation in performance [4][5][6].

Related Work
Over the last four decades, a number of approaches have been proposed to calibrate the gain-and-phase errors. The approaches developed in [7,8] have excellent performance; however, these approaches require a multidimensional search so that they are complicated and time-consuming, and they require calibration signals with known directions. Paulraj and Kailath [9] first investigated the problem of DOA estimation for uniform linear arrays (ULA) in the presence of unknown gain-and-phase errors, and proposed an approach to correct the array perturbations by estimating these errors, but could not provide a unique solution.
The approaches in [10,11] can simultaneously estimate the DOAs of signals and array parameters. Unfortunately, these approaches suffer from suboptimal convergence because of the joint iteration between DOA estimation and array parameter estimation, and they are based on the assumption that the array errors are small.
Wylie et al. [12] also considered the problem in the presence of phase errors, and developed an approach that can obtain a unique solution with some heuristic constraints of phase errors. The primary contribution of this work is a simple, robust (and hence, practically useful) linear least-square (LS) -based algorithm for self-calibration of linear-equispaced (LES) arrays with constant but unknown phase errors that do not require any a priori knowledge of the DOA of any of the received signals, and is consequently useful in situations where a cooperating source may not exist, or is temporarily unavailable. LS estimation of the unknown DOA and array phase response is accomplished via the imposition of a novel constraint, and the resultant estimator bias is shown to be inversely proportional to the number of array elements.
In [13], Li and Er modified the tedious approach in [9] and proposed a class of simplified calibration algorithms based on different diagonal lines of the covariance matrix. They concluded that more diagonals may not result in better performance. Moreover, they presented the theoretical performance of the gain-and-phase-errors estimation. The simplified approach in [13] is superior to the original in [9]. A novel gain-and-phase calibration approach using independent component analysis was proposed in [14], and able to achieve better performance than conventional approaches, but only applicable to no-Gaussian sources.
Zhao et al. [15] proposed a sparse-based DOA estimation approach under sensor gain-and-phase uncertainties, and the approach utilized the sparsity of DOAs and uncertainty matrix and designed an optimization problem for a joint estimation method. Then, an iterative two-step process was proposed to solve the optimization problem, but Zhao's approach is based on the assumption that the array errors are small.
However, these calibration approaches are sensitive to phase errors, and the DOA estimation level drops as the phase errors increase [7][8][9][10][11][12][13][14][15]. This serves as motivation for the introduction of a DOA estimation approach that is robust against phase errors.
The approach in [16,17] can perform independently of the phase errors, but are limited to the configuration of the arrays. These approaches cannot be applied to linear arrays. Moreover, they are not applicable to cases where only one signal is observed.

Motivation and Contribution
Recently, the study on phase retrieval has attracted attention in signal processing [18][19][20]. The problem of phase retrieval refers to recovering a signal from its magnitude measurements, and appears in various fields, such as electron microscopy [21,22], crystallography [23], and optical imaging [24,25], where it is difficult to measure the phase information [26]. To remove the effects of phase errors, Kim et al. first applied the theory of phase retrieval in DOA estimation and proposed a new DOA estimation approach from magnitude-only measurements [27]. Kim's approach does not make use of sensor phase information, so it is robust to phase errors. However, the gain errors were not eliminated in the approach; hence, the DOA estimating performance of Kim's approach becomes degraded as the gain errors increase.
In this paper, we develop a direction-of-arrival (DOA) estimation approach for uniform linear array (ULA) in the presence of gain-and-phase errors based on phase retrieval. The proposed approach can solve the potential impracticality of Kim's approach when gain-and-phase errors coexist. Firstly, we start off by estimating the gain errors, and propose a data preprocessing method for the compensated covariance matrix by taking one column of the matrix, though we include all the magnitude information. In addition, the phase retrieval problem is formed by the proposed data preprocessing method. The estimates of ambiguities are gotten rid of by introducing reference sources.
By combining this with the phase retrieval technique in [28], the proposed DOA estimation approach outperforms Kim's approach in the presence of gain-and-phase errors.

Organizations
The remainder of this paper is arranged as follows. In Section 2, we present the direction-of-arrival estimation system model in the presence of gain-and-phase errors. In Section 3, we propose a new DOA estimation framework based on phase retrieval. Section 4 shows the simulation results of DOA estimation performance for the proposed approach, and conclusions are given in Section 5.

Data Model
Consider a ULA consisting of M sensors with element spacing d. Let λ denote the wavelength of the carrier signal. Assume there are Q far-field narrowband source signals impinging on the ULA. The arriving direction of the qth signal is specified by θ q , q = 1, 2, ..., Q. Then, the steering vector of the ULA is a(θ q ) = [1, e j2πd sin θ q /λ , · · · , e j2π(M−1)d sin θ q /λ ] T .
The signal received by the ULA array is given by the M × 1 vector where A = [a(θ 1 ), a(θ 2 ), · · · , a(θ Q )] is the M × Q steering matrix of the ULA with respect to the directions of the incoming source signals of Q, s(t) is the Q × 1 source vector representing Q source waveforms arriving at the ULA, and n(t) is the M × 1 additive Gaussian noise vector. In this paper, the superscripts T and H represent the transpose and conjugate transpose operations, respectively. For the (2) model, three assumptions are introduced, which are considered to hold throughout the paper.
Assumption 1: The sources are zero-mean, stationary, and unrelated. Assumption 2: The source s(t) is independent of the additive noise n(t). Assumption 3: n(t) is assumed to be stationary, zero-mean, and spatially white Gaussian. Denote the gain-and-phase error of the mth sensor with ρ m and ϕ m , respectively. The received signal in the presence of gain-and-phase errors is where Φ and G are diagonal matrices, and their mth diagonal elements are e jϕ m and ρ m , respectively. Without loss of generality, we assume that ρ 1 = 1.

Gain Errors Estimation from Covariance Matrix
can be expressed as where , σ 2 q denotes the power of the qth source, σ 2 n denotes the power of the noise, and I M is an M × M identity matrix.
The eigendecomposition of R is given as where the eigenvalues {v m } M m=1 are arranged in descending order and u m is the corresponding eigenvector. Then, σ 2 n can be estimated aŝ Define the mth diagonal element R of R(m, m). From (4), we have Recall that ρ 1 = 1, where the gain errors can be estimated aŝ

Gain Errors Calibration
As the gain errors and noise power have been estimated, we are going to find a way to calibrate the gain errors. Since the information on arriving angle is included in the covariance matrix, it can be a feasible way to estimate the arriving angles from the covariance matrix, meaning that we need to eliminate the gain errors in this matrix.
As the noise power is estimated from (6), we can remove the noise in the covariance matrix first: Using the estimated gain errors from (8), we can eliminate the gain errors by compensating for the covariance matrix: whereĜ = diag(1,ρ 2 , . . . ,ρ M ). For the compensated covariance matrix R c , we can find that the impact of gain errors and noise can be removed if they are correctly estimated. The next object is to eliminate the phase errors in the compensated covariance matrix, R c .

A Simple Data Preprocessing Method
For the compensated covariance matrix R c , we propose a new simple data preprocessing method in order to simplify the complexity of the problem. Define the mnth element of R c as r mn ; r mn can then be expressed as (assume d = λ/2) Let B= ∑ Q q=1 σ 2 q e j(π(m−n) sin θ q ) , and Equation (11) can then be written as Based on the Euler Theorem, r mn = B(cos(ϕ m − ϕ n ) + j sin(ϕ m − ϕ n )). Therefore, taking the magnitude squared of r mn , we obtain It is easy to see that the magnitude information of r mn is independent of phase errors, but depends on the phase information related to the arriving angle θ q . This motivates the introduction of a DOA estimation approach based on the magnitude information of r mn . Taking the first column of R c , we obtain an M × 1 vector Then, taking the magnitude squared of r, we obtain where y = |r| 2 and the absolute value operator is applied element-wise to the vector r. The magnitude squared of the nth column of R c can be expressed as where y n = |r n | 2 and r n denotes the nth column of R c . Based on the following equation: we have lemma 1 about the magnitude squared of the elements in the compensated covariance matrix R c .

Lemma 1.
For any element in y n , n = 2...M, there must be an equivalent element in y. Therefore, we can say that y contains all the magnitude information of R c .
Proof. The elements in y are given as ∑ Note that y contains all the magnitude information of R c , and we can use y = |r| 2 , only a column of R c , to form the phase retrieval problem in order to reduce the computational cost. Therefore, we only use one column of R c without missing the magnitude information. The simple data preprocessing method for the compensated covariance matrix helps us to reduce the redundancy of the matrix and form the phase retrieval problem easily.

A New Way to Eliminate the Phase Errors
r can be written as where p = [σ 2 1 , σ 2 2 , . . . , σ 2 Q ] T , and Φ can be regarded as a new phase error matrix: Φ = diag([1, e j(ϕ 2 −ϕ 1 ) , . . . , e j(ϕ M −ϕ 1 ) ]). In order to convert the DOA estimating problem into a nonlinear optimization in a sparse framework, we assume the source's possible DOAs comply with a grid of S points, with S Q. Consider the estimation error for the gain errors and noise power, where (17) can be revised as wherep is defined as an S × 1 sparse vector that only has Q non-zero elements representing source powers, andÃ is defined as a M × S dictionary matrix with a sth column: where θ s refers to a possible DOA. Our goal is to determine the Q non-zero elements ofp from the knowledge of r andÃ. The support of the sparse vectorp has a one-to-one relation with the Q arriving angles θ q , q = 1, 2, ..., Q. Based on Φ Ãp 2 = Ãp 2 , we can obtain Therefore, we have converted the DOA estimation problem into a phase retrieval problem [13]. The rows ofÃ can be regarded as the known measurement vectors. We can now propose estimating the support ofp fromÃ and y, and we can find that the phase errors are totally removed in (21) by taking the magnitude squared of the elements in r.

Phase Ambiguities Solution
However, the phase ambiguity problem associated with phase retrieval solutions should be considered in the DOA estimation. It is significant to resolve the phase ambiguity problem to get a unique solution.
Consider the scenario of two sources, assume the compensated covariance matrix is approximated correctly, and let β q = π(m − 1) sin θ q , where y m can be written as: Assume the source power is known, and the arriving angle can be estimated from Equation (22). To simplify the phase ambiguity problem, let σ 2

Equation (22) then becomes
The solution ∆β to (23) is subject to phase ambiguities, which consist of phase shift ambiguity and phase mirroring ambiguity. The phase shift ambiguity refers to the case that there exists an arbitrary The phase mirroring ambiguity originates from the fact that ∆β and −∆β can be the two solutions to (23). The phase ambiguity can be removed by placing known reference sources, where the DOA of a source can be estimated unambiguously if the source is restricted to only one side of the reference source. For example, if the DOA of the reference source is set at θ re f = 0 rad, the estimating DOA can be estimated unambiguously if it is restricted to the region [0, π/2]. As depicted in Figure 1, the angle ambiguities problem can be solved after a reference target is placed at the cost of a small loss in the visible region. For the multiple-sources DOA estimation problem, the ambiguity can be removed by setting multiple reference sources.

DOA Estimation Using Sparse Least-Square Feasible Point Pursuit (LS-FFP) Approach
Based on (21), the DOA estimation problem can be considered as the optimization problem miñ p y − Ãp 2 2 2 . (24) Let b m denote the mth row ofÃ, where (24) can be recast as where (22) can be recognized as the least-squares (LS) formulation for phase retrieval.
To solve (25), we use the sparse Least-Square Feasible Point Pursuit (LS-FFP) algorithm, a phase retrieval algorithm that exploits sparsity [28]. Firstly, (25) is recast in the following equivalent form: where w = [w 1 , . . . , w M ] T , and the equality constraints in (26) are rewritten as It is clear that (28) is a non-convex constraint. In order to transform the optimization problem into a convex one, following the principles in [29], (28) can be replaced by where z is an S × 1 random vector, and h m ≥ 0 is a slack variable. We thus obtained the following convex, quadratically constrained quadratic programming (QCQP) by: where l 1 -norm p 1 is an ideal description of sparsity, and µ balances the objective function and the sparsity ofp. As indicated in Section 3.5, reference sources should be added to resolve the ambiguity problem. In this paper, the sources are random variables distributed according to the standard Gaussian distribution, so the average power of each source is 1. For the scenario of the single-source DOA estimation, a reference source is set at θ r = 0, and the estimating source is coming from [0, π/2] rad.
The M × S dictionary matrix isÃ = [a(0), a(∆), · · · , a(π/2 − ∆), a(π/2)], and S = π/2∆ + 1. Thus, we modified (31) to the following optimization problem by adding two constraints, miñ p,w,h where p s denotes the power of the sth source, and p 1 = 1 describes the reference source with direction-of-arrival (DOA) θ r = 0 rad. Using the convex optimization toolbox convex (CVX) [30], we were able to solve the convex quadratically constrained quadratic programming (QCQP) (32) in MATLAB. LS-FPP is an iterative algorithm starting with a random initial S × 1 vector z. In the kth iteration, we solve (32) to obtainp k , then setting z k+1 =p k . Since the optimal value of the cost function in each iteration step is non-increasing [28], the iteration stops until the difference of the optimal value between two adjacent iterations is less than a given threshold. Asp is estimated, the DOA estimates can be obtained by seeking the peaks ofp.
Consequently, the proposed approach for DOA estimation is summarized as follows: Step 1: Estimate the covariance matrix with L snapshotŝ Step 2: Estimate the noise power and gain errors by (6) and (8).
Step 3: Calculate the compensated covariance matrix R c as (10).
Step 4: Take the first column of R c and construct the data model as (18).
Step 5: Use the sparse LS-FPP approach to solve the convex QCQP (32) and obtainp.
Step 6: Obtain DOA estimates by seeking the peaks ofp.

Computational Complexity
Computational complexity of the proposed DOA estimation framework lies in three parts: gain errors estimation, compensated covariance matrix construction, and DOA estimation using the Sparse LS-FFP approach. For the gain errors estimation part, the computational complexity is O(M 2 L + M 3 ), which comes from the covariance matrix calculation and eigenvalue decomposition of the covariance matrix. For the compensated covariance matrix construction part, the computational complexity is O(M 3 ), which comes from Equation (10). For the part where the DOA estimation was done using the Sparse LS-FFP approach, the computational complexity mainly comes from solving the optimization problem (32), which is O((S + 3M) 3.5 ). Assuming the maximum iteration is I max , the complexity of the part where DOA estimation is done using the Sparse LS-FFP approach is O(I max (S + 3M) 3.5 ).
The computational complexity of Kim's approach is O(I max (S + M) 2 ). For the approach in [15], it cost O(I max (S + M) 3 ). Figure 2 shows the computational complexity comparison of the proposed approach and other existing approaches. The maximum iteration is assumed to be I max = 10 and S = 30. It can be seen that the proposed approach requires more computational complexity to improve performance.

Numerical Simulations
In this section, numerical simulations are run to demonstrate the performance of the proposed DOA estimation approach. All the simulations are based on a ULA with M = 32 sensors, and the sensor spacing d = 0.5λ. The number of snapshot is L = 100. Let ∆ = π/60 rad, so a grid of S = 31 points is uniformly spaced over the [0, π/2] rad. A reference source of the known DOA is set at θ re f = 0 rad, and the estimating source is located at θ = 2π/15 rad. We also evaluated the performance of the proposed approach for the multiple sources estimation, where two known reference sources were set at θ re f 1 = 0 rad and θ re f 2 = π/15 rad, and the estimating sources are located at θ 1 = 2π/15 rad and θ 2 = π/5 rad. The input signal-to-noise ratio (SNR) of the qth signal is defined as SNR = −10log 10 σ 2 n . The gain-and-phase errors were generated by the following formulas, respectively: where β m and η m are uniform random variables in interval [−0.5, 0.5], σ ρ , and σ ϕ are the standard deviations of gain-and-phase errors, respectively. To compare with Kim's approach, we used the probability of error to evaluate the performance. The probability of error is defined as the ratio of incorrectly estimated DOAs to 100 simulations. Moreover, we compare this with a state-of-the-art sparse DOA estimation approach proposed by Zhao et al. [15]. Firstly, we compared the probability of error of the DOA estimation versus the standard deviation of gain errors (σ ρ ) for the existing approaches and the proposed approach of one or two sources. Here, the SNR was set at 10 dB, and the standard deviation of phase errors (σ ϕ ) was assumed to be 0.3π rad. As shown in Figure 3, the performance of Kim's approach noticeably began to deteriorate as σ ρ increased. For the proposed approach, it has a probability of error lower than 0.1 for σ ρ ≤ 0.75. Therefore, the proposed approach performs more robustly than Kim's approach for a wide range of σ ρ . The proposed approach performs slightly better than Zhao's approach; however, Kim's approach is better when the gain errors do not exist. In the next two experiments, the standard deviation of gain errors is set at 0.5, and Kim's approach is not considered in the next experiments because it is not suitable when the gain errors exist, as we can know from Figure 3. Then, we evaluated the probability of error of the DOA estimation versus SNR for the proposed approach of one or two sources. The standard deviation of phase errors (σ ϕ ) was assumed to be 0.3π rad, and the standard deviation of gain errors (σ ρ ) was set at 0.5. From Figure 4, it was shown that the probability of error of the proposed approach is lower than 0.1 for all SNRs. The proposed approach was shown to provide better performance than Zhao's approach. It is remarkable that the proposed approach is robust to noise. To demonstrate the effect of phase errors on DOA estimation, we evaluated the probability of error of the DOA estimation for the proposed approach with respect to the standard deviation of phase errors from 0 rad to 2π rad. Here, the standard deviation of gain errors (σ ρ ) was set at 0.5, and the SNR was 10 dB. Figure 5 shows that the probability of error of the proposed approach is lower than 0.1 as the standard deviation of phase errors varies. However, the performance of Zhao's approach becomes degraded when σ ϕ > 0.3π rad. Zhao's approach was based on the assumption that the array errors are small [15], and it is indicated that the phase errors have no effect on the proposed approach, which is a key benefit based on phase retrieval. In order to explore the results better, we also evaluated the root mean square error (RMSE) on the DOA estimates. We compared the RMSE of the DOA estimates versus the standard deviation of gain errors for the existing approaches and the proposed approach. The SNR was set at 10 dB, and the standard deviation of phase errors was assumed to be 0.3π rad. We also compared with another existing techniques in [13]. As shown in Figure 6, the proposed approach performs better than other existing approaches. Then, we compared the RMSE of the DOA estimates versus the standard deviation of gain errors for the existing approaches and the proposed approach. The SNR was set at 10 dB, and the standard deviation of gain errors was set at 0.15. Shown in Figure 7, we can see that the proposed approach and Kim's approach perform independently of phase errors. Kim's approach even performs slightly better than the proposed approach, because the standard deviation of gain errors is relatively small in this experiment. However, Zhao's approach and Li's approach perform badly when the standard deviation of phase errors rise to 0.5π rad.

Conclusions
In this paper, we presented a new DOA estimation approach for ULA in the presence of gain-and-phase errors based on phase retrieval. In contrast to the existing DOA estimation approach using phase retrieval, the proposed approach can provide more robust DOA-estimating performance as the gain errors increase. Meanwhile, because the effect of phase errors is eliminated by taking the magnitude squared of the elements in the compensated covariance matrix, the proposed approach performs independently of phase errors.

Conflicts of Interest:
The authors declare no conflict of interest.