Two Novel Two-Stage Direction of Arrival Estimation Algorithms for Two-Dimensional Mixed Noncircular and Circular Sources

This paper addresses the two-dimensional (2D) direction-of-arrival (DOA) estimation problem with two novel methods for mixed noncircular and circular signals. The first proposed method is named the two-stage direction-of-arrival matrix (TSDOAM) method, and the other is called the two-stage rank reduction (TSRARE) method. The proposed methods utilize both the circularity and the direction-of-arrival differences between the noncircular and circular sources to estimate the 2D directions-of-arrival (DOAs). The maximum detectable 2D angle parameters of the TSDOAM and TSRARE methods are twice those of the existing methods. Moreover, the TSRARE method can detect more incident signals than the TSDOAM method due to the array aperture of two parallel uniform linear arrays (ULAs) being fully utilized. Simulation results show that compared to the existing methods for the small angle separation of 2D directions-of-arrival, the two proposed methods perform well in terms of the signal-to-noise ratio (SNR) and snapshots.

In order to improve the direction-of-arrival estimation performance, many effective noncircular algorithms for 2D directions-of-arrival have been presented in [31][32][33]. In [31], Liu et al. proposed an extended rank reduction (ERARE) method with noncircular information exploited for two-parallel uniform linear arrays (ULAs) which achieved an improved estimation accuracy compared to [15]. Based on [25], the conjugate information of the observed data was utilized to realize a better 2D direction-of-arrival estimation [32]. use of the conjugate information of the observed data to realize a better 2D direction-of-arrival estimation. A method that applied noncircular direction finding to the hexagonally-shaped electronically steerable parasitic antenna radiator (ESPAR) array was presented in [33], and the Cramér-Rao bound (CRB) was analyzed for comparison. However, the aforementioned algorithms cannot cope with the direction-of-arrival estimation problem for the mixed noncircular (e.g., binary phase shift keying, BPSK) and circular (e.g., quaternary phase shift keying, QPSK) signal scenario. Although Yin's method [15] and Xia's method [16] could be applicable to the above-mentioned mixed signal scenario, the distinguishable signals were less than the array elements.
In [34][35][36], direction-of-arrival estimation schemes for joint noncircular and circular signal estimation were proposed for 1D direction finding. In [34], Gao. et al. proposed a method that constructed a new data vector with the original data and the conjugate ones to form two estimators for noncircular and circular signal finding, respectively. However, the method in [34] cannot deal with the coincident directions-of-arrival of noncircular and circular signals, and its estimation performance degraded severely in small angle separation . In addition, the maximum number of detected signals was still limited. An improved method was presented in [35] to solve the above problems, which adopted the direction-of-arrival circularity difference rather than the direction-of-arrival difference between the noncircular and circular signals to estimate the direction-of-arrival. Nevertheless, the direction-of-arrival estimation performance dropped due to few observed data being available. In [36], a sparse representation method for mixed signals was proposed by exploiting overcomplete dictionaries that were subject to the sparsity constraint to jointly represent the covariance and elliptic covariance matrices of the array output. However, for 2D situations, there are few research works for joint noncircular and circular signal direction finding. Additionally, much of the work in array processing has also been focused on optimization problems, such as genetic algorithms [37][38][39][40].
In this paper, inspired by the method in [35], two novel 2D direction-of-arrival estimation algorithms using two parallel ULAs with a two-stage direction-of-arrival matrix (TSDOAM) method and a two-stage rank reduction (TSRARE) method, separately, are proposed for mixed noncircular and circular signals estimation. The direction-of-arrival circularity difference rather than the direction-of-arrival difference between the noncircular and circular signals for the 2D directions-of-arrival's estimation is utilized in the two proposed methods. The maximum number of distinguished mixed signals of the two proposed methods are identified compared to the conventional methods, which show that the detected number of the signals is more than that of the array elements. Moreover, when both of the 2D directions-of-arrival are incoming from small angle separation-even when both of them are from the overlapping direction-the estimation accuracy of the two proposed methods is better than Yin's and Xia's method.
The rest of this paper is organized as follows. The array signal model is introduced in Section 2. The TSDOAM and TSRARE methods are described in detail in Section 3. The maximum number of detective mixed signals is analyzed in Section 4. Simulation results are presented to verify the performance of the two proposed methods in Section 5. Conclusions are drawn in Section 6.

Array Signal Model
As shown in Figure 1, suppose that there are K = K n + K c (assume the number of mixed signals is known to the receiver) uncorrelated far-field sources that are mixed K n noncircular sources s n,k (t) and K c circular sources s c,k (t) from direction (θ k , β k ), k = 1, 2, . . . , K, impinging on two parallel ULAs with each array having M elements. The distance between the two arrays is λ/2, denoted as d y , and the interelement spacing d x on each array is also λ/2, where λ is the wavelength of the incident waves. The additive noises of two ULAs are circular Gaussian noises, which are uncorrelated with the incoming signals. The observed data vectors X(t) = [x 1 (t), x 2 (t), . . . , x M (t)] T and Y(t) = [y 1 (t), y 2 (t), . . . , y M (t)] T from two parallel ULAs are given by: where A is called the steering vector with each column denoted a(θ k ) and a(θ k ) = [a 0 (θ k ), . . . , a M−1 (θ k )] T , whose element can be expressed as , and the element v(β k ) has the form of e j 2π λ d y cos β k . S(t) = [s 1 (t), s 2 (t), . . . , s K (t)] denotes the radiating signal vector. N x (t) = [n x,1 (t), . . . , n x,M (t)] T and N y (t) = [n y,1 (t), . . . , n y,M (t)] T represent the circular Gaussian noise vectors of the two arrays, respectively.

The TSDOAM Method
A novel method called the TSDOAM method with the two-stage direction-of-arrival matrix (DOAM) method and direction-of-arrival circularity difference, is proposed in this part. According to Equations (1) and (2), and based on the assumption that the noise and the signals are uncorrelated and that the mixed signals are also uncorrelated, the auto-covariance matrix R xx and cross-covariance matrix R yx can be written, respectively, as follows.
where A n and A c denote the steering matrices associated with noncircular and circular signals, separately. σ 2 is the variance of the circular Gaussian noises, and R n and R c have the form of R n = diag{ E[s n,1 s * n,1 ], ..., E[s n,K n s * n,K n ]} and R c = diag{ E[s c,1 s * c,1 ], ..., E[s c,K c s * c,K c ]} , respectively. In practice, non-circularity and circularity are important properties of a random variable; their concept comes directly from the geometrical interpretation of a complex random variable. The source would be called a circular source if its statistical characteristics have a rotational invariance characteristic; otherwise, it would be called a noncircular source. Based on this, we only consider the rotational invariance characteristic of the first-and second-order statistical properties of the sources. For a complex random source s k , we define E[s k ], E[s k s * k ], and E[s 2 k ] as the mean, the covariance, and the elliptic covariance of the source s k , respectively. For an arbitrary phase ϕ k as follows: If the source's first-and second-order statistical properties are rotationally invariant, the source is identified as circular; otherwise, it is determined to be noncircular. Therefore, the elliptic auto-covariance matrix R xx can be expressed as follows: Notice that for a complex circular random variable h, E[hh] = 0 [6]. Therefore, the circular component and the circular Gaussian noise component of Equation (8) both equal zero, and the elliptic auto-covariance matrix R xx can be rewritten as: Similarly, the elliptic cross-covariance matrix R yx is computed as follows: We then estimate the 2D directions-of-arrival of noncircular signals with Equations (9) and (10). First, let {η n,1 , ..., η n,K n } and {v n,1 , ..., v n,K n } be the eigenvalues and corresponding eigenvectors of R xx , respectively, namely: The pseudo-inverse of R + xx is: Due to R n being a diagonal matrix and A n a column full-rank matrix, we attain the following formula with Equation (9): Combining Equation (13) with (10), we obtain an alternative expression of R yx : Right-multiplying both sides of Equation (14) by R + xx A n , we get: Substituting Equations (11) and (12) into Equation (15), From Equations (1) and (2), we can see that the dimensions of the observed data vectors X(t) and Y(t) are both M × 1. Then, based on Equations (3) and (4), it is easy to know that the dimensions of the auto-covariance matrix R xx and cross-covariance matrix R yx are both M × M. From Equation (8), we can see that the dimension of the elliptic auto-covariance matrix R xx is also M × M. Therefore, R xx is a square matrix. Specifically for the first proposed method in this paper, the R + xx -which is the pseudo-inverse of R xx -is equivalent to the inverse of R xx . That is to say, the R xx and R + xx in the first proposed method are both square matrices. Moreover, the diagonal elements of R + xx are nonzero elements , and R + xx is full rank. Therefore, is an identity matrix, and Equation (16) can be simplified as: From Equation (17), 2D directions-of-arrival of noncircular signals-which are obtained by performing eigenvalue decomposition (EVD) of R yx R + xx , denoted as the direction-of-arrival matrix, lie in A n and Φ n , respectively.
where ξ n,k and u n,k are the eigenvalues and the corresponding eigenvectors of R yx R + xx , respectively. It can be verified that the spanned subspace from the steering matrix A n and the signal subspace U n = [u n,1 , . . . , u n,K n ] are the same.
Unlike Yin's method obtaining the 1D angle by the spectrum search with a certain region (which entailed high complexity), here, define h n,k = u n,k /u n,k (1); we can get: Integrating the expression a i (θ n,k ) and v(β n,k ) with Equations (18) and (19), the estimated 2D directions-of-arrival of noncircular signals are achieved as follows: In the next stage, the 2D direction-of-arrival of the circular signals can be obtained with the estimates θ n,k and β n,k above.
With Equation (9), R n can be estimated as: where [41]. Let the kth diagonal element of R n be R n (k, k); we get: here, we assume that the noncircular signals are BPSK-modulated signals. Therefore, E[s n,k s n,k ] = σ 2 n,k e jϕ n,k , where σ 2 n,k = E[s n,k s * n,k ] and ϕ n,k are the noncircular phases. It is easily deduced that σ 2 n,k = |R n (k, k)|. Due to R n being diagonal, it follows that: Then, let R xx,1 = A n R n A H n and R yx,1 = A n Φ n R n A H n ; we attain: Define: We have another direction-of-arrival matrix R yx,2 R + xx,2 that is related to circular signals. Just as the way of attaining the 2D directions-of-arrival of noncircular signals, the 2D directions-of-arrival θ c,k and β c,k of circular signals are obtained with the direction-of-arrival matrix R yx,2 R + xx,2 using the same direction-of-arrival matrix (DOAM) method. Until now, the TSDOAM method is summarized as follows.
Step 1: Calculate R xx , R yx , R xx , and R yx from Equations (3)-(10); Step 2: Execute the EVD of R xx to get its pseudo-inverse matrix R + xx ; Step 3: Perform the EVD of R yx R + xx with Equation (18); Step 4: Attain θ n,k and β n,k using Equations (20) and (21); Step 5: Construct the matrix R n with the estimate R xx and A n ; Step 6: Construct the direction-of-arrival matrix R yx,2 R + xx,2 with Equations (27) and (28); Step 7: Repeat Step 2 to Step 4 for the θ c,k and β c,k .

The TSRARE Method
In this section, in order to make full use of array elements of two ULAs, another novel method called the TSRARE method-which is based on the direction-of-arrival circularity difference and the rank reduction (RARE) method-is proposed in the two-stage estimation procedure.
By concatenating the observed data vectors X(t) and Y(t), we get: where C = [c(θ 1 , β 1 ), . . . , c(θ K , β K )] is termed as the extended steering matrix, and: As the radiating mixed signals are uncorrelated with each other, the conjugated covariance matrix R of Z(t) can be written as: and the elliptic covariance matrix R of Z(t) is as follows: where C n and C c denote the noncircular and circular extended steering matrix, respectively. Next, we perform singular-value decomposition (SVD) of R to estimate the 2D directions-of-arrival of noncircular signals as follows: where Λ n denotes the K n × K n diagonal matrix containing K n nonzero singular values in the diagonal. It is verified that all of the columns of C n are orthogonal to all of the columns of U n,2 ; that is, C H n (θ n,k , β n,k )U n,2 U H n,2 C n (θ n,k , β n,k ) = 0, k = 1, . . . , K n .
Associated with Equation (30), Equation (34) can be rewritten as: where r = [1 v T (β n,k )] T , and: T(θ n,k ) = a H (θ n,k )U n,1 a(θ n,k ) a H (θ n,k )U n,2 a(θ n,k ) a H (θ n,k )U n,3 a(θ n,k ) a H (θ n,k )U n,4 a(θ n,k ) where U n,1 = U n,21 U H n,21 , U n,2 = U n,21 U H n,22 , U n,3 = U n,22 U H n,21 , and U n,4 = U n,22 U H n,22 , U n,21 and U n,22 are obtained by dividing U n,2 into the same two-dimensional matrices-namely U n,2 = [U T n,21 U T n,22 ] T . In order to use the RARE method, we define a(θ n,k ) as: where b n,k = e −j 2π λ d x cos θ n,k , and the matrix T is a function of b n,k . Then, the 1D directions-of-arrival θ n,k can be obtained by finding the values of b n,k such that det[T(b n,k )] = 0. Additionally, the polynomial of b n,k has the following form: where m n,p = a T (1/b n,k )U n,p a(b n,k ), p = 1, 2, 3, 4. m n,p is the polynomial of b n,k whose lth coefficient is given by the sum of the elements of the lth diagonal of U n,p , where l = −M + [Φ n, f ] i,l−i+1 , Φ n,1 = µ n,1 µ T n,4 and Φ n,2 = µ n,2 µ T n,3 . Hence, Equation (38) can be expressed as: The roots of the polynomial det[T(b n,k )] can be computed by exploiting the computationally-efficient polynomial root multiple signal classification (MUSIC) algorithm, and the 1D noncircular angles θ n,k are obtained as: Substituting the estimated θ n,k into Equation (35), we then seek out the minima of the following function [42]: From Equation (41), we obtain β n,k that are given by the eigenvector corresponding to the smallest eigenvalue associated with T(θ n,k ) as: Next, we estimate the 2D directions-of-arrival of circular signals based on the estimates θ n,k and β n,k . With Equations (30) and (32), R n,1 -denoted as the estimated R n in the TSRARE method-can be expressed as: [41]. Denoting the kth diagonal element of R n,1 as R n,1 (k, k), we have: Similarly, the noncircular signals here are BPSK-modulated signals. It is easily deduced that σ 2 n,k = |R n,1 (k, k)|. Because R n is diagonal, it follows that R n,1 , which is the estimate of R n as follows: n,1 , . . . , σ 2 n,K n ] = diag[|R n,1 (1, 1)|, . . . , |R n,1 (K n , K n )|]. (45) Then, let R 1 = C n R n,1 C H n ; we attain: Perform the EVD of R − R 1 as follows: where U c,x and U c,z are called the signal and noise subspaces associated with the signal eigenvalue matrix Λ c,x and noise eigenvalue matrix Λ c,z , respectively. Similarly, we get the following relationship: and estimate the 2D directions-of-arrival (namely θ c,k and β c,k ) of circular signals in the same way as the 2D directions-of-arrival of noncircular signals. At this point, the 2D directions-of-arrival of noncircular and circular signals have been achieved by the TSRARE algorithm. The simple summary of the TSRARE algorithm is as follows.
Step 1: Calculate R and R with the observed data Z(t) with Equations (31) and (32); Step 2: Perform the SVD of R to get U n,2 using Equation (33); Step 3: Calculate the roots of det[T(b n,k )] associated with noncircular signals from Equation (34) to Equation (39); Step 4: Attain θ n,k and β n,k using Equations (40)

and (42);
Step 5: Construct R 1 with the estimate R n,1 and C n ; Step 6: Perform the EVD of R − R 1 to get U c,z using Equation (47); Step 7: Repeat Step 3 to Step 4 for the θ c,k and β c,k .

Remark 1:
In practice, it can be noted that only a finite number of observed data is available.
Thus, R xx , R yx , R xx , R yx , R, and R must be estimated byR xx = 1 Z(l)Z T (l). When less observed data are available, such that the estimated covariance matrices are no longer strictly diagonal matrices, this leads to the estimation performance reduction in the small 2D angle separation condition.
Remark 2: As we know, one of the awkward problems about 2D direction-of-arrival estimation is the pair situation, which causes severe estimation error without exact pair process. However, the two proposed methods can pair the 2D directions-of-arrival of mixed signals automatically. This is because, in the two-stage estimation procedures, the TSDOAM method pairs the 2D directions-of-arrival by performing one EVD whose eigenvalues and eigenvectors are a one-to-one correspondence relationship, while the TSRARE method pairs them by decoupling the 2D directions-of-arrival into two successive 1D processes.

Location Discussion and Analysis
The maximum number of sources to estimate is analyzed in this section by the two proposed methods. Since the two proposed methods estimate the mixed signals separately, and the number of sources can be resolved is related to the dimensions of the elliptic auto-covariance matrix R xx (see Equation (8)) and the conjugated covariance matrix R (see Equation (31)), it follows that M > max {K n , K c } with the TSDOAM method and 2M > max {K n , K c } with the TSRARE method must be satisfied, respectively, so as to resolve all noncircular and circular signals, while Yin's [15] and Xia's [16] method estimate them simultaneously, so M > (K n + K c ) and 2M > (K n + K c ) should be satisfied, respectively. In other words, the maximum detectable signals of the TSDOAM method that are detected by two ULAs with 2M elements are (M − 1) noncircular signals plus (M − 1) circular signals, namely (2M − 2), and the TSRARE method can detect (2M − 2) noncircular signals plus (2M − 2) circular signals, namely (4M − 4); while Yin's and Xia's method can distinguish (M − 1) and (2M − 2) mixed signals, respectively. Therefore, the TSDOAM and TSRARE methods can detect twice the mixed signals as Yin's and Xia's method, respectively. In addition, the TSRARE method can identify twice the signals as the TSDOAM method.

Simulation Results
In this section, some simulation results are presented to show the performance of the TSDOAM and TSRARE methods, compared with the existing methods, which include Yin's method and Xia's method. Assume two ULAs with each array consisting of omnidirectional sensors spaced by a half wavelength of the mixed noncircular and circular signals; the distance between the two ULAs is spaced by a half wavelength as well. Additionally, the noncircular signals employ BPSK-modulated sources, while the circular signals are QPSK-modulated sources for the simulation.

2D Direction-of-Arrival Estimation Performance
In this part, the maximum detectable 2D directions-of-arrival are investigated. Four BPSK signals with θ n,k and β n,k impinge from {70 • , 60 • , 85 • , 100 • } and {65 • , 80 • , 85 • , 65 • }, separately, and four QPSK signals are emitted from the same 2D directions-of-arrival. The number M of isotropic sensors of each array is three. Figure 2 plots the paired results of eight radiating signals from 50 Monte Carlo trials with the signal-to-noise ratio (SNR) set at 30 dB and snapshots L = 2000, which show that the 2D directions-of-arrival of eight (4M − 4 = 8) signals are paired correctly with the TSRARE method; even a common β is shared in both noncircular and circular signals. However, the TSDOAM method and the methods in [15,16] can detect up to four (2M − 2 = 4) signals, four (2M − 2 = 4) signals, and two (M − 1 = 2) signals, respectively, which fail to distinguish the mixed eight signals due to the limited array elements.

The Effect of SNR
In this subsection, we compare the 2D direction-of-arrival performance of the TSDOAM and TSRARE methods versus SNR with existing methods in [15,16]. Furthermore, the average root mean square error (ARMSE) is defined for precision evaluation as: where ζ k stands for θ k or β k , andζ qk is the parameter to be estimated forθ k orβ k , while M c denotes the number of Monte Carlo runs. Two BPSK signals together with two QPSK signals incoming from (75 • , 50 • ), (100 • , 65 • ) and (75 • , 50 • ), (100 • , 65 • ), respectively, impinge onto the ULAs with each array having five omnidirectional sensors. The snapshots of this test are set to 500, and the variable SNR of the four incident signals varies from −10 dB to 30 dB. The ARMSE of the four methods derived from the 2000 trials are given in Figure 3. It can be seen from Figure 3 that Yin's method and Xia's method fail to work; however, the two proposed methods perform well with increasing SNR. This is because the two proposed methods estimate the 2D directions-of-arrival of noncircular and circular signals separately based on the circularity difference between noncircular and circular signals rather than the direction-of-arrival difference utilized in Yin's and Xia's methods. Moreover, the TSRARE method has a lower ARMSE than the TSDOAM method, which results from the fact that the array aperture is fully utilized with the TSRARE method during the two-stage estimation for the BPSK and QPSK signals.

The Effect of Snapshots
The 2D direction-of-arrival performance of the TSDOAM and TSRARE methods versus snapshots with Yin's and Xia's methods is verified in this part. The simulation parameters are the same as Experiment 2, except the SNR is fixed at 5 dB, and the number of the collected snapshots varies from 10 to 490. At each snapshot setting, 2000 independent runs are executed for each method to obtain their 2D direction-of-arrival estimation in the statistical sense, whose results are given in Figure 4. Similar conclusions and reasons can be drawn from Figure 4 that as snapshots are increasing, the curves of the two proposed methods work well, while those of Yin's and Xia's methods remain unchanged. In addition, the TSRARE method has better estimation performance than the TSDOAM method in both 2D directions-of-arrival.

The Effect of Angle Separation
In this part, we testify to the 2D direction-of-arrival estimation performance of the two proposed methods coupled with Yin's and Xia's methods versus angular separation. The number of each array of the two ULAs is also five. Additionally, four signals consist of two BPSK signals and two QPSK signals with θ k and β k incoming from {70 • , 100 • , (100 + ∆) • , (70 + ∆) • } and {50 • , 85 • , (85 + ∆) • , (50 + ∆) • }, separately, and ∆ is varied from 0 to 12. In addition, the SNR is fixed at 15 dB, and the snapshots L = 300.
The ARMSE versus angular separation is shown in Figure 5 with 2000 Monte Carlo trials used. From Figure 5, it can be seen that the two proposed methods outperform Yin's and Xia's methods in both 2D small angle separations; however, with the angle separation increasing, the 2D direction-of-arrival estimation performance achieved by the two proposed methods is inferior to that by Yin's and Xia's methods. This is because in small angle separation, the 2D direction-of-arrival information of circular signals included in the elliptic covariance matrix can be equivalent to that of noncircular signals, due to the fact that the 2D direction-of-arrival of the circular signal is close to that of the noncircular signal; thus, the estimation precision of the noncircular signals can be improved. A similar reason is suitable for the improved 2D direction-of-arrival estimation accuracy of circular signals; while Yin's and Xia's methods estimate the 2D direction-of-arrival of noncircular and circular signals simultaneously, which are based on the direction-of-arrival difference of noncircular and circular signals that inevitably gives rise to the degradation of the performance in small 2D direction-of-arrival separation [35]. As for large 2D direction-of-arrival separation, the estimation performance of the two proposed methods behaved with lower accuracy than Yin's and Xia's methods. The reason is that in the condition of a few snapshots, the incident mixed signals' covariance matrix and elliptic covariance matrix are not strictly diagonal matrices, which results in making it difficult to separate the mixed signals from the two-stage estimation procedures. Furthermore, the TSRARE method performs better than the TSDOAM method all the way, with the full usage of array elements in the two-stage estimation procedures.

Conclusions
Two novel 2D direction-of-arrival estimation methods, named, respectively, as the TSDOAM and TSRARE methods for mixed noncircular and circular signals' estimation with two parallel ULAs are presented in this paper. The direction-of-arrival circularity difference rather than the direction-of-arrival difference between the noncircular and circular signals for the 2D directions-of-arrival' estimation is utilized in the two proposed methods. The explicit derivation of the two proposed methods is described, and the maximum number of incident signals of the two proposed methods is analyzed, which shows that the detected number of signals is more than that of the array elements compared to the conventional methods. Simulation results demonstrate the usefulness of the two proposed methods.