A Low-Complexity Method for Two-Dimensional Direction-of-Arrival Estimation Using an L-Shaped Array

In this paper, a new low-complexity method for two-dimensional (2D) direction-of-arrival (DOA) estimation is proposed. Based on a cross-correlation matrix formed from the L-shaped array, the proposed algorithm obtains the automatic pairing elevation and azimuth angles without eigendecomposition, which can avoid high computational cost. In addition, the cross-correlation matrix eliminates the effect of noise, which can achieve better DOA performance. Then, the theoretical error of the algorithm is analyzed and the Cramer–Rao bound (CRB) for the direction of arrival estimation is derived . Simulation results demonstrate that, at low signal-to-noise ratios (SNRs) and with a small number of snapshots, in contrast to Tayem’s algorithm and Kikuchi’s algorithm, the proposed algorithm achieves better DOA performance with lower complexity, while, for Gu’s algorithm, the proposed algorithm has slightly inferior DOA performance but with significantly lower complexity.


Introduction
Direction-of-arrival (DOA) estimation, which has found its potential applications in the fields of sonar, radar, wireless communication, etc, is an important research branch of array signal processing [1]. Two-dimensional (2D) [2,3] direction-of-arrival (DOA) estimation with different structured arrays, such as two-parallel arrays [4][5][6][7], L-shaped arrays [8][9][10][11][12][13][14][15], and uniform rectangular array [16,17] has drawn a remarkable amount of attention. In [18], it has been proven that the estimation performance of the L-shaped array prevails over many other structured arrays. Therefore, there has been growing interest in 2D DOA estimation utilizing the L-shaped arrays. Tayem and Kwon [12] presented computationally efficient 2D angle estimation with a propagator method using one or two L-shaped arrays. Unfortunately, this method cannot pair the 2D angles automatically and may cause a matching failure problem. Consequently, a pair-matching method using a cross-correlation matrix was proposed to remove the aforementioned problem in [13]. However, the correct estimation of the incoming "virtual angles" [12] was the fatal problem at a low signal-to-noise ratio (SNRs) and with a small number of snapshots, which seriously affected the estimation performance of 2D DOAs.
A method [14] based on joint singular value decomposition (JSVD) of two cross-correlation matrices (CCMs), which mitigated the additive noise effect, was put forward to estimate elevation and azimuth parameters without additively pairing procedures. However, the method required heavy calculation due to SVD operation and "beamforming-like" spectral search operation. A two-dimensional angle matching algorithm based on the estimated signal covariance matrix is proposed in [19]. When the signal source is coherent, it can be achieved by minimizing a cost function constructed by the two covariance matrices. This method is robust to the CCM-ESPRIT algorithm. Tayem [20] divided two uniform arrays on the L-matrix into two subarrays and computed the cross-covariance matrices on the two uniform arrays, respectively. Then, adding the two mutual covariance matrices with their transpose matrix, respectively, we can obtain two new cross-covariance matrices. By segmenting these two new matrices, we can get the two-dimensional angle estimation with linear operation. However, the method still requires a two-dimensional angle matching process. By using the conjugate symmetry of two uniform linear array patterns on the L-array, the effective aperture of arrays can be extended in [21], and, then, the automatic matching of the two-dimensional angle parameters based on the PM algorithm and ESPRIT algorithm can be obtained, which avoids the cumbersome peak searching process. Therefore, the method not only has good direction finding precision, but also has the advantage of low complexity. A novel cumulants-based approach [22] to 2D DOA estimation for coherent non-Gaussian sources with two parallel ULA (uniform linear arrays) is presented. It has a lesser amount of calculation, which avoids constructing several FOC (fourth order cumulants)-based sub-matrices to form two full rank spatially smoothed matrices. When two close coherent signals are present, it is more effective and efficient than the FOC-FSS (fourth-order cumulants-based forward spatial smoothing) method in 2D DOA estimation in both white noise and color Gaussian noise situations. Wu [23] proposed a novel 2D direct-of-arrival and mutual coupling coefficients estimation algorithm for uniform rectangular arrays. The algorithm can achieve a better performance than those auxiliary sensor-based ones. It first built a general mutual coupling model that is based on banded symmetric Toeplitz matrices and then used the rank-reduction method to solve the 2D DOA estimation problem. With the obtained DOA information, the mutual coupling coefficients can be estimated.
Chen [24] derived a series of 2D DOA estimators with a new data vector that combines the received array data and its conjugate counterparts for mixed circular and non-circular signals based on a 2D array structure consisting of two parallel ULAs. However, it can give a more accurate estimation when the number of sources is within the traditional limit of high resolution methods and still work effectively when the number of mixed signals is larger than that of the array elements. In addition, it avoids the complicated 2D spectrum peak search and therefore has a much lower computational complexity. A multiresolution approach [25] for the DOA estimation of multiple signals based on a support vector classifier has been presented. This method defines a probability map of the incidence of an electromagnetic signal and performs a synthetic zoom in the angular sector iteratively. Then, it is able to estimate the DOAs of a number of sources greater than the maximum allowed by conventional eigenvalue decomposition methods for a fixed planar array geometry, and provide good results dealing with both a single signal and multiple signals.
In this paper, based on CCMs, a new pair-matching algorithm is presented to achieve 2D angles with low complexity. Firstly, the elevation angles are estimated by a linear operation of the cross-correlation matrix formed from an L-shaped array, and then the corresponding azimuth angles are achieved by the interrelationship between the elevation and azimuth angle without an additional paired procedure. Moreover, the Cramer-Rao bound (CRB) for 2D DOAs of an L-shaped array is studied. The complexity advantage of the proposed algorithm is analyzed, which is significant as sensors and snapshots increased. Furthermore, the theoretical error of the proposed algorithm is derived.
The rest of this paper is organized as follows. Section 2 presents the array signal model. Description of the proposed algorithm is introduced in Section 3. Section 4 analyzes the complexity of the proposed algorithm. The theoretical error analysis of the proposed algorithm is derived in Section 5. The analysis of the CRB of the L-shaped array is given in Section 6. The experimental results are compared with several existing approaches in Section 7. Finally, Section 8 concludes this paper.

Array Signal Model
As illustrated in Figure 1, K far-field narrowband plane wave signals s i (t), i = 1, . . . , K, impinge on the L-shaped array structured by two uniform orthogonal arrays in the x-z plane. Each array consists of N identical omni-directional sensors separated by λ/2 inter-element spacing d, namely, d = λ/2, where λ is the wavelength of the incident waves. The ith source has an elevation angle θ i and an azimuth angle ϕ i . The observed signal vectors at the sub-arrays along the x-axis and z-axis are written in matrix form as T are the N × 1 received signal vectors along the x-axis and z-axis, respectively.
T are the Gaussian white noise vectors along the x-axis and z-axis, respectively. In addition, are denoted as N × K array manifold matrices of the x-axis and z-axis, respectively. N × 1 array manifold vectors cos θ i along the x-axis and z-axis, respectively. We suppose that the source signals are non-Gaussian and uncorrelated to each other; the Gaussian noises with zero-mean and variance σ 2 are statistically independent to the signals.

The Proposed Algorithm
Firstly, a cross-correlation matrix R xz is obtained as follows: where . From Equation (3), it can be noted that the additive noise is removed by the cross-correlation operation. Let R xz1 and R xz2 be the first and last N − 1 columns of R xz , and we have where and A z2 (θ) denote the first and last N − 1 rows of A z (θ), respectively. When the incoming signal covariance matrix R s is diagonal matrix, Equation (5) can be rewritten as By combining Equations (4) and (6), a new 2N × (N − 1) matrix R is defined as follows: Then, we partition the matrix A xe (ϕ, θ) in Equation (7) as where Similarly, we partition R in Equation (7) sub-matrix R 2 , which has the following relationship In practice, the propagator matrix P is achieved by minimizing the following cost functions where · 2 F signifies Frobenius norm. The estimate of P is as follows: To maximize usage of array information, we introduce an extended propagator matrix P e as follows: where I K is the K × K identity matrix. In the noiseless case, right-multiplying by A xe1 (ϕ, θ) of Equation (13), we obtain Next, we partition P e into two N × K sub-matrices P e1 and P e2 , and Equation (14) can be rewritten as According to Equation (15), we get Then, we introduce a new matrix ψ that can be expressed as In Equation (18), by performing eigen-value decomposition (EVD) of ψ, eigenvaluesα i and corresponding eigenvectors A 1 that correspond to the diagonal elements of Λ H (θ), and the estimate of A xe1 (ϕ, θ) can be achieved, respectively. Here, we denote where Ω is a permutation matrix with Ω −1 = Ω.
Then, according to the expression of Λ H (θ), the elevation angles are as follows: In addition, using P e11 to denote the first (N − 1) rows of P e1 , P e12 to denote the last (N − 1) rows of P e1 , P e21 to denote the first (N − 1) rows of P e2 , and P e22 to denote the last (N − 1) rows of P e2 , respectively, we define With the assumption that A = A xe (ϕ, θ)Ω, we know that P e11 A 1 , P e12 A 1 , P e21 A 1 , P e22 A 1 are the first (N − 1) rows, the second to N-th row, the (N + 1)-th to (2N − 1)-th row, the last (N − 1) rows of A , respectively, so which contribute to whereΦ = ΩΦΩ -1 with Φ = diag(e −j(2π / λ)d cos ϕ 1 , . . . , e −j(2π / λ)d cos ϕ K ) . In addition, the azimuth angles lie in the diagonal elementsβ i of B + 1 B 2 as follows: At this point, 2D elevation and azimuth parameters have been automatically paired by EVD operation. The summary of the proposed algorithm is shown as follows: Step 1: Compute R xz and R from Equations (3) and (7).

Complexity Analysis
As for the complexity, we analyze on the basis of the matrix complex multiplication, which mainly involves in auto-correlation or cross-correlation matrix construction, EVD or SVD operation, pseudo-inverse operation, and "beamforming-like" spectral search. Define the search step of azimuth ϕ ∈ [0, 180 • ] with ∆ϕ = 0.01 • . The major computations of the proposed algorithm is , while Tayem's algorithm [12], Kikuchi's algorithm [13], and Gu's algorithm [14] , respectively, where L denotes the number of snapshots. Due to sample snapshots L >> N > K, therefore, the proposed algorithm has lower complexity than others. Figure 2a,b shows the complexity comparison between the proposed method and other methods. From both Figure 2a,b, we find that the proposed method has lower computational load than others as sensors and snapshots increase.

Theoretical Performance Analysis
The perturbation is caused by noise in the proposed method, and we analyze on the basis of the matrix perturbation theory [26,27].
Let X = X + ∆X, Z = Z + ∆Z, and the covariance matrix with perturbation be expressed as where ∆R xz is the perturbation of the covariance matrix.
Then, R xz1 = R xz1 + ∆R xz1 , R xz2 =R xz2 + ∆R xz2 where ∆R xz1 , ∆R xz2 is the first and last N − 1 columns of ∆R xz . From Equation (7), we can get where ∆R 1 , ∆R 2 are the first K rows and the last (2N − 1) rows of ∆R, respectively. From Equation (10), we get (P + ∆P) H (R 1 + ∆R 1 ) = R 2 + ∆R 2 , according to P H R 1 = R 2 and, neglecting the second-order term ∆P H ∆R 1 , we can get ∆P H The extended propagator matrix P e is as follows: and P e1 = P e1 + ∆P e1 , P e2 = P e2 + ∆P e2 , where ∆P e1 , ∆P e2 are the first and last N rows of ∆P. According to Equation (18), ψ = ψ + ∆ψ, ψ = P + e1 P e2 . Similar to Equation (29), we can get ∆ψ = P + e1 (∆P e2 − ∆P e1 ψ). By performing EVD of ψ with perturbation, the influence to eigenvalues α i can be expressed as where v i and u i stand for the left and right orthogonal eigenvectors associated with α i of ψ, respectively.
Let φ i = arg(α i ). Then, Equation (20) can be written as θ i = arccos(φ i λ/2πd), θ i = θ i + ∆θ i . The perturbations of elevation ∆θ i can be obtained according to the theorem of first-order approximation of multivariate function [28]. Specific content is as follows.
For z close to x, the first-order approximation of f near x can be represented as: where ∇ f (x) denotes the gradient of f and is a column vector. Thus, we can get where For Equations (21) and (22), the perturbations are and A 1 = A 1 + ∆A 1 , where ∆A 1 is the estimation error of A 1 . According to Equation (25), Φ = Φ + ∆Φ, and ∆Φ = B + 1 (∆B 2 − ∆B 1 Φ) can be obtained with a similar method to Equation (29). It can be easily obtained that the perturbations of diagonal elements β i are the diagonal elements of ∆Φ. Let ζ i = arg(β i ), and the perturbations of the azimuth can be expressed as follows, similar to elevation: Therefore, the root mean-squared error of two-dimensional direction of arrival estimations are ). (37)

Cramer-Rao Bound (CRB) Analysis
In the case of L-shaped array configuration, the Cramer-Rao bound (CRB) of 2D DOAs is considered here. Rewrite the received data from L-shaped array as The Fisher information matrix (FIM) F with respect to ϕ = [ϕ 1 , ϕ 2 , . . . , ϕ K ] and θ = [θ 1 , θ 2 , . . . , θ K ] can be written as Note that the (i, j)-th element of F 11 [29] is given by Similarly, we get the (i, j)-th element of F 12 , F 21 and F 22 , respectively, as follows: where Re(·) denotes the real part, e i denotes the i-th column of the identity matrix, trace(·) denotes the trace of a matrix and M ij denotes the (i, j)-th element of M,Ȧ ς m ,Ȧ ς (m = i, j, ς = ϕ, θ), R s and γ has the form ofȦ In Equation (47), Q has different expressions for different type of noises as below: for white noise, P, for unknown noise, where I N denotes the N × N identify matrix, and the (p, q)-th element of the unknown noise covariance matrix P is 0.8 |p−q| e j(p−q)π 2 . According to Equations (40) to (43), we obtain where ⊗ denotes the Hadamard matrix product. Then, the CRB matrix C can be expressed as and we can obtain the CRB of azimuth and elevation parameters as follows: where C i,i denotes the (i, i)-th element of C. Therefore, we define the CRB for the parameters of the i-th source as

Experimental Results
Simulation experiments are conducted in this part. In all experiments, the elements spacing of L-shaped array is λ/2.
In the first experiment, we examine the scattergram of 2D elevation and azimuth of the proposed algorithm compared with that of the Kikuchi algorithm in both white noise and unknown noise situations. The number of isotropic sensors N is 5. Two uncorrelated equal power signals with elevation θ i and azimuth ϕ i incoming separately from (55 • , 65 • ) and (75 • , 80 • ). In addition, their SNRs are set to 20 dB and the number of snapshots are fixed at 300. Five hundred independent trials are carried out. Figures 3 and 4 show that 2D DOA statistic performance of the proposed algorithm is better than the Kikuchi algorithm, especially in an unknown noise situation. In addition, pairing failures are emerging in Figures 3b and 4b. The reason is that the noise factor in the proposed algorithm has been removed, and the difference between "virtual angles" is small in the Kikuchi algorithm when pair-matching is required.    In the second experiment, the proposed algorithm in theoretical analysis and experimental studies, Tayem's algorithm, Kikuchi's algorithm, Gu's algorithm and CRB are compared in terms of root mean square error (RMSE) with respect to SNRs and snapshots in white noise situations. Define RMSE as  From Figures 5 and 6, it can be noted that the theoretical estimation performance of the proposed algorithm is better than the experimental at low SNR, and, with the increase of SNR and snapshots, they gradually overlap together. In addition, the proposed algorithm is better than Tayem's algorithm, Kikuchi's algorithm, but slightly inferior to Gu's algorithm at low SNR and with a small number of snapshots. As the SNR and snapshots increased, the estimation performance of the proposed algorithm is close to Gu's algorithm with lower computational cost, which avoids SVD of the cross-correlation matrix R and "beamforming-like" spectral search.
In the third experiment, the proposed algorithm in theoretical analysis and experimental studies, Tayem's algorithm, Kikuchi's algorithm, Gu's algorithm and CRB are compared in terms of RMSE with respect to SNRs and snapshots in an unknown noise situation. The parameters configured in this experiment are the same as the second experiment. Figures 7 and 8 show the 2D DOA statistic performance in an unknown noise situation. Apparently, as shown in Figures 7 and 8, similar conclusions can be drawn. From Figures 7 and 8, it can be noted that the trend of theoretical and experimental estimation performance of the proposed algorithm is the same as Figures 5 and 6. Then, we can get that the DOA estimation performance of Tayem's algorithm and Kikuchi's algorithm deteriorates seriously because Tayem's algorithm and Kikuchi's algorithm are sensitive to the type of noise. In addition, the estimation performance of the proposed algorithm is roughly the same as Gu's algorithm at low SNR and with a small number of snapshots. At high SNR and with a large number of snapshots, the estimation performance of the proposed algorithm is very close to Gu's algorithm with lower computational cost.

Conclusions
A novel low-complexity method for 2D angle parameter estimation is proposed in this paper. The explicit description of the proposed method is derived to achieve the automatic pairing 2D angle parameters. In addition, the theoretical performance analysis and CRB of 2D DOAs is given. Simulation results show the effectiveness of the proposed algorithm in contrast to other algorithms, especially at low SNR and with a small number of snapshots.