Fast Noncircular 2D-DOA Estimation for Rectangular Planar Array

A novel scheme is proposed for direction finding with uniform rectangular planar array. First, the characteristics of noncircular signals and Euler’s formula are exploited to construct a new real-valued rectangular array data. Then, the rotational invariance relations for real-valued signal space are depicted in a new way. Finally the real-valued propagator method is utilized to estimate the pairing two-dimensional direction of arrival (2D-DOA). The proposed algorithm provides better angle estimation performance and can discern more sources than the 2D propagator method. At the same time, it has very close angle estimation performance to the noncircular propagator method (NC-PM) with reduced computational complexity.


Introduction
Two-dimensional direction of arrival (2D-DOA) estimation has been widely used in mobile communication systems, sonar, navigation, radar, etc. [1][2][3][4][5], which is an important research branch in array signal processing. Many 2D-DOA estimation algorithms have sprung up in recent years in order to improve the performance of angle estimation, which include the two dimensional multiple signal classification(2D MUSIC) algorithm [6], the 2D Unitary estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm [7], the modified 2D-ESPRIT algorithm [8], the matrix pencil method [9], the maximum likelihood method [10,11], the parallel factor (PARAFAC) algorithm [12], and so on [13][14][15][16][17][18][19][20]. However, those 2D-DOA estimation algorithms are confronted with the problem of the high computational complexity generally and they are very difficult to apply in engineering practice. As is known to us, the propagator method (PM) algorithm uses linear operations to replace the eigenvalue decomposition of the covariance matrix [21], and it has a great advantage in resolving the amount of calculation. Therefore, the 2D-DOA estimation based on PM is becoming a hot spot of research. For example, Wu et al. have developed the 2D-DOA estimation algorithm via the rotational invariance property of propagator matrix [22]. In [23], an improved PM algorithm is proposed for 2D-DOA estimation, which not only reduces the computational complexity, but also avoids the aperture loss.
Unfortunately, all the algorithms mentioned above did not consider the characteristics of the impinging signals. In fact, many noncircular signals such as the amplitude modulated (AM), binary phase shift keying (BPSK), minimum shift keying (MSK), and Gaussian MSK (GMSK) signals are used in wireless communication or satellite systems. In recent years, some scholars use non-circular signal characteristics to improve the performance of direction estimation, which contain the noncircular MUSIC (NC-MUSIC) algorithm [24], the NC-ESPRIT algorithm [25], and the noncircular parallel where A x = [a x (v 1 ), a x (v 2 ), · · · , a x (v K )] and a x (v k ) = [1, e −jπv k , · · · , e −j(M−1)πv k ] T , v k = cos φ k sin θ k , θ k is the elevation angle and φ k is the azimuth angle. Φ = diag(e −jπu 1 , e −jπu 2 , · · · , e −jπu K ) and u k = sin φ k sin θ k . S(t) = [s 1 (t), s 2 (t), · · · , s K (t)] T is the noncircular signal vector. In addition, the vector of strictly second-order noncircular signals can be expressed as [28]: S(t) = ΛS o (t), S o (t) ∈ C K×1 , S o (t) = S o * (t), and Λ = diag e jϕ 1 , e jϕ 2 , · · · e jϕ K , (e jϕ p = e jϕ q , f or p = q). n i (t) is the additive white Gaussian noise vector of the ith subarray. Therefore, the whole array output is where A = A y • A x is the MN × K steering vector matrix, • represents the Khatri-Rao product, and A y = [a y (u 1 ), a y (u 2 ), · · · , a y (u K )], a y (u k ) = [1, e −jπu k , · · · , e −j(N−1)πu k ] T , and

Euler Transformation
The real part and imaginary part of (t) x can be obtained by utilizing the real-valued property of noncircular signals and Euler's formula as follows: Then, we define a new virtual array data as follows:

Euler Transformation
The real part and imaginary part of x(t) can be obtained by utilizing the real-valued property of noncircular signals and Euler's formula as follows: x where n R (t) = [n(t) + n * (t)]/2, n I (t) = [n(t) − n * (t)]/(−2j), Then, we define a new virtual array data as follows: where Define two matrices as follows: , and we can get the following relationship: where is a real-valued matrix whose diagonal elements contain the needed angle information: Similarly, define two (M − 1) × M Toeplitz matrices as follows: Then, we construct two matrices J 3 and J 4 as follows: We also get the following relationship: where D = diag tan( πv 1 2 ), tan( πv 2 2 ), · · · , tan( πv K 2 ) is a real-valued matrix whose diagonal elements also contain the desired angle information.

2D-DOA Estimation
According to Equation (5), the estimation of covariance matrix R of y(t) is denoted by collecting L snapshots: From Equation (8), R can be denoted by In the noiseless case, R x1 p r = R x2 , an estimation matrix p r can be obtained by [21]:p We construct a new matrix p s = I K P H r , where I K is the identity matrix. In the noiseless case, the relationship between p s and B can be obtained by a unique non-singular matrix T as Substituting Equation (10) into Equation (6), we can get If we define P 1 = (J 1 p s ) † J 2 p s , we then have Equation (12) shows that the diagonal elements of the matrix G can be obtained by performing the eigenvalue decomposition of P 1 , and T is the corresponding eigenvector.
Then, we can get the estimation ofû k : wherep k is the kth diagonal element of the matrix G.
Similarly, Substituting Equation (10) into Equation (7), we can also get If we define P 2 = (J 4 p s ) † J 3 p s , we then have Then, we get the estimation ofv k :v wherer k is the kth diagonal element of the matrix D.
We note thatû k andv k share the same eigenvector T, so the pairing is automatically formed. Thus, 2D-DOA can be obtained byθ We have now achieved the essence of the proposed algorithm. The major algorithmic steps are as follows: (1) Construct the matrix y(t) from Equation (5), and compute the covariance matrix R of y(t) through Equation (8).
(2) Estimation of the propagator p r from Equation (9), and then construct the matrix p s .
(3) Construct the matrix J 1 p s and J 2 p s and perform the eigenvalue decomposition of Similarly, construct the matrix J 3 p s and J 4 p s and perform the eigenvalue decomposition of (17) and (18).

Remark 1.
In [23], the conventional PM algorithm divides the steering matrix A into two matrices A 1 ∈ C K×K and A 2 ∈ C (MN−K)×K , and A 2 is the linear transformation of A 1 , i.e., is the propagator operator. According to Equation (1), x(t) = AS(t) + n(t), and the covariance matrix In our paper, according to Equation (5), y(t) = x R (t) x I (t) = BS o (t) + n r (t), and we compute the covariance of y(t) ∈ R 2MN×1 to estimate the propagator. Apparently, the available array aperture of the proposed algorithm can be thought of as twice that of the conventional 2D-PM [23], so it has better angle performance than 2D-PM.

Remark 2.
In [23], define p c = I K P H m , and then p c A 1 = A, which means that the columns in p c span the same signal subspace as the column vectors in A. Divide p c into p c 1 ∈ C M(N−1)×K and p c 2 ∈ C M(N−1)×K , p c 1 , p c 2 are the first M(N − 1) rows and the last M(N − 1) rows of p c . Then, get the relationship, P + c1 P c2 = A 1 ΦA 1 , where Φ = diag(e −jπu 1 , e −jπu 2 , · · · , e −jπu K ). Perform the eigenvalue decomposition of P + c1 P c2 to obtain the diagonal elements of the matrix Φ. Similarly, reconstruct p c to p c , p c 1 , p c 2 being the first N(M − 1) rows and the last N(M − 1) rows of p c , and perform the eigenvalue decomposition of P c1 + P c2 to obtain the diagonal can discern more sources than that of the conventional 2D-PM [23].

Remark 3.
In the NC-PM algorithm [27], according to Equation (2), the extended array output data is denoted with ones on its anti-diagnoal and zeros elsewhere, and X * ∈ C MN×L stands for the complex conjugation of X, Y ∈ C 2MN×L . Compute the covariance of Y to estimate the propagator p nc . Similarly, the invariance equations for p nc are constructed to estimate the 2D-DOA. As is known to us, each computation amount of the complex multiplication is four times that of the real-valued one. In our algorithm, we use Euler transformation to convert complex arithmetic of noncircular to real arithmetic. For example, according to Equation (5), y(t) = BS o (t) + n r (t), y(t) ∈ R 2MN×1 , and the computation amounts of covariance of y(t) with snapshots L are much lower than that of Y [27]. Due to real-valued processing, our algorithm can save about 75% computational load compared with the NC-PM algorithm [27].

CRB
In this section, we give the Cramer-Rao Bounds (CRB) of noncircular signal for rectangular planar array. According to Equation (5), the received data is where B = [b 1 , b 2 , · · · , b K ] ∈ C 2MN×K , and n r (t) is the noise vector. The Fisher information matrix (FIM) in relation to φ = [φ 1 , φ 2 , · · · , φ K ] and θ = [θ 1 , θ 2 , · · · , θ K ] can be calculated as follows [29]: According to [29], we know that the (i, j) ith element of F 11 is given by Likely, we can give the (i, j) ith element of F 12 , F 21 , F 22 : where e i denotes the ith column of the unit matrix, , σ 2 is the covariance of the noise. According to Equations (21)-(24), we can obtain: where ⊕ represents Hadamard product. Then, the CRB can be denoted as: We present the curves of CRB versus different signal to noise ratios (SNRs) and snapshots L in Figures 2 and 3. The source number K is fixed at 3 M and N represents the numbers of sensors on the x-axis and the y-axis. In Figure 2, the snapshot L is fixed at 200. It is obvious that, with the improvement of SNR, the value of CRB decreases accordingly. In Figure 3, we set SNR at 20 dB, and the curve shows that the value of CRB decreases with increase of L, and simulation results and theory analysis are consistent.

Complexity Analysis
In this section, we analyse the computational complexity of the algorithm specifically. First, estimation of the covariance matrix R requires RMS. Therefore, the overall computational complexity of our algorithm is

Complexity Analysis
In this section, we analyse the computational complexity of the algorithm specifically. First, estimation of the covariance matrix R requires RMS. Therefore, the overall computational complexity of our algorithm is

Complexity Analysis
In this section, we analyse the computational complexity of the algorithm specifically. First, estimation of the covariance matrix R requires O(4M 2 N 2 L) real-valued multiplications (RMS). In addition, the estimation of the matrix p s takes O(2MNK 2 + 4M 2 N 2 K + K 3 ) RMS. Then, the estimation and eigenvalue decomposition of the matrix P 1 and P 2 totally require O(4M 2 NK(N − 1) + 3M(N − 1)K 2 + 8M(M − 1)N 2 K + 2K 3 ) RMS. Therefore, the overall computational complexity of our algorithm is O(4M 2 N 2 L + 5MNK 2 + 16M 2 N 2 K + 3K 3 − 4M 2 NK − 3MK 2 − 8MN 2 K) RMS. As we know that each computation amount of the complex multiplication is four times that of the real-valued one, we can show the Chen's noncircular propagator algorithm [27] needs O(16M 2 N 2 L + 8MNK 2 + 16M 2 N 2 K + 3K 3 + 16M(N − 1)K 2 ) RMS, J's noncircular ESPRIT [25] needs The complexity comparisons with different parameters are shown in Figures 4 and 5. In Figure 4, the numbers of sensor M and N on the x-axis and the y-axis are set at 8 and 6, respectively. The source number K is fixed at 3. In Figure 5, the parameters N and K are the same as Figure 4, and the snapshot L is set to 100. From Figures 4 and 5, we can observe that the proposed algorithm has much lower computational load than J's NC-ESPRIT algorithm and Chen's NC-PM algorithm.   We can summarize the merits of the proposed algorithm as follows: (1) The proposed algorithm has much lower computational load than the NC-PM and NC-ESPRIT algorithms because the proposed algorithm uses Euler transformation to convert complex arithmetic of noncircular PM to real arithmetic.    We can summarize the merits of the proposed algorithm as follows: (1) The proposed algorithm has much lower computational load than the NC-PM and NC-ESPRIT algorithms because the proposed algorithm uses Euler transformation to convert complex arithmetic of noncircular PM to real arithmetic.  We can summarize the merits of the proposed algorithm as follows: (1) The proposed algorithm has much lower computational load than the NC-PM and NC-ESPRIT algorithms because the proposed algorithm uses Euler transformation to convert complex arithmetic of noncircular PM to real arithmetic. (2) The proposed algorithm has better estimation performance than the 2D-PM algorithm because the array aperture is doubled according to Equation (5).

Simulation Results
In this section, we use Monte Carlo simulations to verify the performance of the algorithm. In the simulation, the rectangular planar array is configured with N subarrays, each subarray contains M sensors, L is the snapshots of the sources, and K is the number of the sources. We assume that there are K = 3 non-coherent sources, which are BPSK modulated in where The root mean squared error (RMSE) is used for performance assessment, which is defined as

Simulation Results
In this section, we use Monte Carlo simulations to verify the performance of the algorithm. In the simulation, the rectangular planar array is configured with N subarrays, each subarray contains M sensors, L is the snapshots of the sources, and K is the number of the sources. We assume that there are K = 3 non-coherent sources, which are BPSK modulated in Figures 4-8 The root mean squared error (RMSE) is used for performance assessment, which is defined as    Figure 5a,b, we know that the proposed algorithm has better RMSE performance than Li's algorithm [23]. Furthermore, it has close RMSE performance to Chen's algorithm [27]. However, we should know that our algorithm has much lower computational amount than J's NC-ESPRIT algorithm and Chen's NC-PM algorithm owing to the real-valued processing, which means that it is more suitable for a practical application system. Figure 8 presents RMSE performance comparisons at different snapshots L. Where M = 8, N = 6, SNR is varied from 0 dB to 20 dB. We can observe that the RMSE performance is improved with the increase of snapshot L. When L increases, we get more samples to estimate the propagator matrix more accurately, and so the angle estimation performance is enhanced.    Figure 5a,b, we know that the proposed algorithm has better RMSE performance than Li's algorithm [23]. Furthermore, it has close RMSE performance to Chen's algorithm [27]. However, we should know that our algorithm has much lower computational amount than J's NC-ESPRIT algorithm and Chen's NC-PM algorithm owing to the real-valued processing, which means that it is more suitable for a practical application system.      Where M = 8, N = 6, SNR is varied from 0 dB to 20 dB. We can observe that the RMSE performance is improved with the increase of snapshot L. When L increases, we get more samples to estimate the propagator matrix more accurately, and so the angle estimation performance is enhanced. Figures 9 and 10 present RMSE versus different values of M or N, respectively. The snapshot L is fixed at 200. In addition, it is indicated that RMSE performance is improved when M or N increases. Multiple sensors enhance the aperture of the array as well as diversity gain. Therefore, it can improve the angle estimation performance.
Figures 9 and 10 present RMSE versus different values of M or N, respectively. The snapshot L is fixed at 200. In addition, it is indicated that RMSE performance is improved when M or N increases. Multiple sensors enhance the aperture of the array as well as diversity gain. Therefore, it can improve the angle estimation performance.
The estimation performance for two closely spaced sources is also investigated. Figure 11 depicts the scatter plot of 2D-DOA estimation results for two closely spaced sources. Where M = 8, N = 10, SNR = 10 dB, the snapshot L is 200. It is shown that our algorithm works well for the closely spaced sources.   In addition, it is indicated that RMSE performance is improved when M or N increases. Multiple sensors enhance the aperture of the array as well as diversity gain. Therefore, it can improve the angle estimation performance.
The estimation performance for two closely spaced sources is also investigated. Figure 11 depicts the scatter plot of 2D-DOA estimation results for two closely spaced sources. Where M = 8, N = 10, SNR = 10 dB, the snapshot L is 200. It is shown that our algorithm works well for the closely spaced sources.   The estimation performance for two closely spaced sources is also investigated. Figure 11 depicts the scatter plot of 2D-DOA estimation results for two closely spaced sources. Where M = 8, N = 10, SNR = 10 dB, the snapshot L is 200. It is shown that our algorithm works well for the closely spaced sources.

Conclusions
We have presented a novel direction finding algorithm for uniform rectangular planar array. The characteristics of noncircular signal and Euler's transformation are exploited to get the real-valued rectangular array data in a new way. The proposed algorithm can reduce the computational amount since it does not refer to plural operation and the eigenvalues' decomposition of the covariance matrix. The theory analysis and simulation results verify that our algorithm is more suitable for real-time processing system in engineering.

Conclusions
We have presented a novel direction finding algorithm for uniform rectangular planar array. The characteristics of noncircular signal and Euler's transformation are exploited to get the real-valued rectangular array data in a new way. The proposed algorithm can reduce the computational amount since it does not refer to plural operation and the eigenvalues' decomposition of the covariance matrix. The theory analysis and simulation results verify that our algorithm is more suitable for real-time processing system in engineering.