# Crosscorrelation and DOA Estimation for L-Shaped Array via Decoupled Atomic Norm Minimization

Authors: Yu Zhang; Yinan Sun; Gong Zhang; Xinhai Wang; Yu Tao
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

A novel two-phase method for two-dimensional (2D) direction-of-arrival (DOA) estimation with L-shaped array based on decoupled atomic norm minimization (DANM) is proposed in this paper. In the first phase, given the sample crosscorrelation matrix, the gridless DANM technique considering the noise and finite snapshots effects is employed to exploit the structure and sparse properties of the crosscorrelation matrix. The resulting DANM-based algorithm not only enables the crosscorrelation matrix reconstruction (CCMR) but also reconstructs the covariance matrix of the L-shaped array. Hence, sequentially, in the second phase, the conventional 2D DOA estimators for the L-shaped array can be adopted for the angle estimation. With appropriate 2D DOA estimators, the resulting proposed algorithms can not only achieve better performance but also detect more source number, compared with conventional crosscorrelation-based DOA estimators. Moreover, the proposed method, termed CCMR-DANM, not only has blind characteristic that it does not require the prior information of source numbers but also is more efficient than the existing CCMR-based counterparts. Numerical simulations demonstrate the effectiveness and outperformance of the proposed method.


---

## Body

## 1. Introduction

The problem of two-dimensional (2D) direction-of-arrival (DOA) estimation plays an important role in array signal processing and has attracted much interest in the area of wireless communications, radar and sonar [1–7]. For 2D DOA estimation, many array structures, such as rectangular arrays, circular arrays, and L-shaped arrays, have been developed. Among these arrays, since the L-shaped array can achieve better estimation performance than others, it has attracted a lot of attentions and many corresponding algorithms for 2D DOA estimation have been proposed in last several decades [8–15]. Moreover, these algorithms can be divided into three categories. The first is to separately estimate the angles corresponding to each uniform linear subarray based on the covariance matrix of each subarray with conventional 1D DOA estimators, such as MUSIC [16] and ESPRIT [17]. However, an extra pairing operation is needed in these algorithms [18]. The second is to jointly estimate the two angles based on the covariance matrix of the L-shaped array. They can detect more source numbers than the first ones and do not need an extra pairing [8, 9]. The last is based on the crosscorrelation of the L-shaped array, which is naturally contaminated by less noise, compared with the covariance matrix. As a result, the corresponding algorithms can achieve a better performance in low signal-to-noise ratio (SNR) [10–15]. It is worth noting that in practical applications, all these three kinds of methods need the prior information of source numbers and are employed with sample matrices (no matter the sample covariance matrices or the sample crosscorrelation matrices). Note that the sample matrix is calculated from finite collected snapshots of the observed signals which is contaminated with additive noise. Hence, the sample matrix cannot equip the ideal structure of the ideal matrix, which leads to performance degradation especially in low SNR and small number of snapshots [19].To overcome these disadvantages, the basic idea is to first reconstruct the ideal matrix from the sample matrix before employing DOA estimation methods. Inspired by this thought, several crosscorrelation matrix reconstruction- (CCMR-) based methods are proposed. Specifically, crosscorrelation atomic norm minimization (CC-ANM) utilising the 2D ANM technique [20] for CCMR is proposed in [21]. And [22] grafts the covariance fitting criterion [23] for CCMR and proposes a crosscorrelation gridless sparse iterative covariance-based estimation (CC-GLS) method. However, a high-dimensional two-level Toeplitz matrix is needed to be constructed in both methods, which leads to high computational complexity [18]. In contrast, we propose an efficient decoupled ANM- (DANM-) based method for CCMR in [24], while the method does not consider the finite snapshot effect, which results in performance degradation in practical applications.In this paper, given the sample crosscorrelation matrix, the DANM technique [18] considering the noise and finite snapshots effects is employed to exploit the structure and sparse properties of the crosscorrelation matrix. Moreover, the resulting DANM-based algorithm not only enables the crosscorrelation matrix reconstruction but also reconstructs the covariance matrix of the L-shaped array. Hence, the conventional 2D DOA estimators for the L-shaped array can be sequentially adopted for the angle estimation. The proposed two-phase method is more computationally efficient than the aforementioned two CCMR-based methods and can achieve better estimation performance compared with traditional crosscorrelation-based methods. Numerical simulations demonstrate the effectiveness and outperformance of the proposed method.The rest of this paper is organized as follows. Section2 presents the signal model and problem formulation for 2D DOA estimation with L-shaped array. Section 3 proposes an efficient DANM-based two-phase method for crosscorrelation and DOA estimation. Section 4 discusses the computational complexity and extends the proposed method to the sparse L-shaped array cases. Section 5 presents simulation results, followed by conclusions in Section 6.Throughout this paper,a, a and A denote a scalar, a vector, and a matrix, respectively. TrA denotes the trace of A. Ta denotes the Hermitian Toeplitz matrix with the first column being a. diagA is a column vector formed from the elements of the main diagonal of A and diaga generates a diagonal matrix with the diagonal elements constructed from a. vec· stacks all the columns of a matrix into a vector. Ia is an a-size identity matrix, and IΩ is a selection matrix with index set Ω. 0a and 1a are the a×1 zeros and one vectors, respectively. ⊗ is the Kronecker product. E· denotes expectation, and Var· denotes variance. We use ·T, ·∗, and ·H to denote the transpose, the conjugate, and the conjugate transpose operation, respectively.


## 2. Problem Formulation

ConsiderK far-field narrowband source signals siti=1K impinging on an L-shaped array from distinct directions at element angles ϕii=1K and azimuth angles θii=1K, as shown in Figure 1. The L-shaped array consists of two ULAs of M omnidirectional sensors which are uniformly spaced with a spacing of d along the x axis and y axis, respectively. The observed signals of the L-shaped array can be expressed as [24]
(1)xt=Axst+nxt,yt=Ayst+nyt,t=1,2,⋯,L,where t indexes the snapshot; L denotes the number of collected snapshots; and st, nxt, and nyt denote the vector of source signals and the vector of additive noise corresponding to the x subarray and y subarray at the snapshot t, respectively. Ax=axθ1,ϕ1,⋯,axθK,ϕK and Ay=ayθ1,ϕ1,⋯,ayθK,ϕK are the array manifold matrices of the x subarray and the y subarray, whose ith columns are the steering vectors of the ith source which satisfy
(2)axθi,ϕi=1,e−j2π1/2cosθisinϕi,⋯,e−j2π1/2M−1cosθisinϕiT,ayθi,ϕi=1,e−j2π1/2sinθisinϕi,⋯,e−j2π1/2M−1sinθisinϕiT.Figure 1
L-shaped array configuration.Herein, the spaced distanced is assumed to be equal to half of the wavelength λ. Moreover, let fi,1=1/2cosθisinϕi, fi,2=1/2sinθisinϕi and Ω=f1,1,f1,2,⋯,fK,1,fK,2 denote the frequencies that correspond to the direction on x and y axes and the set of corresponding frequencies, respectively. Note that fi,1,fi,2⟷θi,ϕi is the one-to-one mapping. Once the estimation of fi,1,fi,2 is obtained, the corresponding θi,ϕi can be retrieved as
(3)ϕ^i=arcsin2f∧i,12+f∧i,22,θ^i=arccosf^i,2f^i,1.Hence, in the following paper, we consider estimation offi,1,fi,2 instead of θi,ϕi for notational simplicity. Moreover, in this paper, the source signals st are assumed uncorrelated with each other and the noise nxt and nyt are i.i.d. additive white Gaussian random processes satisfying N0,σ2IM and are statistically independent of st. Therefore, under the above assumptions, we have the crosscorrelation matrix between xt and yt can be expressed as
(4)Rxy=ExtyHt=AxRsAy+N1=Rc+N1=∑i=1Kriaxfi,1ayHfi,2+N1,where N1 is the matrix with only the first element of main diagonal being σ2 and zero otherwise. Rs=diagr is the source correlation matrix with r=r1,⋯,rKT≥0, i.e., Esit2=ri. Rc=AxRsAy is the noise-free crosscorrelation matrix, while, in practical applications, Rxy can only be estimated from the finite L snapshots by R^xy=1/L∑t=1LxtyHt, which not only is contaminated by the noise but also contains errors caused by the finite snapshot effect. The goal of this paper is to recover the noise-free crosscorrelation matrix Rc and then the unknown 2D DOAs θi,ϕii from the sample crosscorrelation matrix R^xy.


## 3. Proposed Method

In this section, we propose an efficient DANM-based crosscorrelation and DOA estimation method for L-shaped array. To this end, the standard DANM technique is firstly presented. Then, the sparse representation of the noise-free crosscorrelation matrix is presented, which enables the DANM technique to exploit its sparse property. Further, simultaneously considering the noise and finite snapshot effects, and the structure property of the crosscorrelation matrix, we propose the original DANM-based formulation for CCMR which is intractable since the structure constraint. To make the formulation tractable, an effective relaxation is proposed. Moreover, an estimation error constraint leading to easy setting of the user-specific parameter is also proposed.


### 3.1. Prior Art: DANM Technique

We now review the standard DANM technique for harmonic retrieval. Define a matrixZ as
(5)Z=∑k=1Kαkaxfk,1ayHfk,2,where αk∈ℂ. Then, based on (5), a matrix-form atom set of infinite size is defined as [18, 15]
(6)A=axf1ayHf2∣f1,f2∈−12,12.Accordingly, the atomic norm ofZ over the atom set A is defined as
(7)∥Z∥A=inf∑kαk∣∑kαkaxfk,1ayHfk,2,axfk,1ayHfk,2∈A;αk∈ℂ,which seeks the sparsest (under l1-norm measure) decomposition of Z over A. Consider the matrix Z is contaminated by the noise matrix N and the matrix at hand is Z^=Z+N. Then, according to the DANM theory [18, 25], the atomic decomposition yields the true structure in (5), through the following DANM formulation:
(8)Z~=argminZZAs.t.Z^−ZF2≤η,where η is a user-specified parameter for error tolerance. Moreover, (8) is equivalent to the following semidefinite positive (SDP) formulation
(9)z~1,z~2,Z~=argminz1,z2,Z12MTrTz1+TrTz2s.t.Z^−ZF2≤ηTz1ZZHTz2≽0.


### 3.2. Standard DANM-Based Formulation for CCMR

Note thatRc=∑i=1Kriaxfi,1ayHfi,2. Apparently, Rc has a sparse linear atomic representation over the matrix-form atom set in (6). And we introduce a new matrix-form atomic norm
(10)RcA+=inf∑iri∣∑iriaxfi,1ayHfi,2,axfi,1ayHfi,2∈A;ri≥0,∀i.Note that this norm differs from the original atomic norm of the standard DANM in (7), because of the extra constraint r≥0.GivenRc, it is possible to retrieve the components ri,fi,1,fi,2 of its sparest representation by calculating its atomic norm. Considering the obtained at hand is the sample crosscorrelation containing the estimation error, it boils down to
(11)R~c=argminRcRcA+s.t.FR^xy,Rc≤β.where FR^xy,Rc is the function quantifying the estimation error and β indicates the error tolerance threshold. It is similar to the DANM formulation introduced in (8), but defined on the norm ·A+ instead of ·A. To show the intricacy of this difference, we rewrite (11) in the following equivalent form:
(12a)R~c,r~=minRc,rRcA=r1(12b)s.t.FR^xy,Rc≤β(12c)Rc=AxdiagrAyH(12d)r≥0.Note that (12c) is implicit in the objective function (11) but becomes an explicit constraint because of the new nonnegative constraint (12d). Without (12d), the SDP implementation of the DANM in (9) can be used to reformulate (12) into a convex problem. However, because of the extra constraint on r in (12d), this problem becomes intractable, because r is intertwined with the other variable Rc in the form of (12c).


### 3.3. Effective Relaxation

To solve (12) in a tractable manner, we seek to relaxr≥0 to an effective form with respect to Rc. To this end, we note that
(13)Rc:,1=∑i=1Kriaxfi,1,Rc1,:H=∑i=1Kriayfi,2,where Rc:,1 and Rc1,: denote the first column and the first row of Rc, respectively. Moreover, according to the property of Toeplitz matrices [26], if r≥0, we have
(14)TRc:,1≽0,TRc1,:H≽0.Adopting (14) to replace (12d) and reformulating (12a)–(12c) into the original decoupled SDP form in (9), we reach the following effective SDP relaxation for (11):
(15)u~1,u~2,R~c=argminu1,u2,Rc12MTrTu1+TrTu2s.t.FR^xy,Rc≤βTRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0.With an appropriate definition ofFR^xy,Rc and β, (18) can be solved successfully via off-the-shelf convex solvers, such as CVX [27].


### 3.4. Estimation Error Constraint

To define a niceFR^xy,Rc which leads to easy setting of β, we denote the estimation error matrix as
(16)E=R^xy−Rc.Denote byεp,q the p,qth element of the estimation error matrix E; then, one has the following proposition:Proposition 1.
For adequately largeL, and ∀p+q−2>0, εp,q is approximately circular complex Gaussian distributed with zero mean, and the variance is
(17)Varεp,q=1L∑i=1Kri+σ22−∑i=1Kri2.p+q−2>0,1≤p,q≤M.Proof.
Note thatR^xy is estimated from the L collected snapshots. Denote R^p,q as the p,qth element of R^xy. xpt, yqt, sit, nxpt, and nyqt are similarly defined. We have ∀p+q−2>0,
(18)R^p,q=1L∑t=1Lxptyq∗t=1L∑t=1L∑i=1Ksite−j2πp−1fi,1+nxpt∑j=1Ksj∗tej2πq−1fj,2+nyq∗t=∑k=1K1L∑t=1Lskt2e−j2πp−1fk,1ej2πq−1fk,2+1L∑t=1L∑i=1K∑j=1j≠iKsitsj∗te−j2πp−1fi,1ej2πq−1fj,2+1L∑t=1L∑i=1Ksitnyq∗te−j2πp−1fi,1+∑j=1Ksj∗tnxptej2πq−1fj,2+1L∑t=1Lnxptnyq∗t=R^p,q0+R^p,q1+R^p,q2+R^p,q3,where R^p,q0, R^p,q1, R^p,q2, and R^p,q3 denote the first, the second, the third, and the fourth summand. Note that ER^p,q0=∑i=1Krie−j2πp−1fi,1ej2πq−1fi,2=Rcp,q, where Rcp,q is the p,qth element of Rc. Hence, the p,qth estimation error εp,q can be appropriately estimated as εp,q=R^p,q1+R^p,q2+R^p,q3. Since L is sufficiently large, εp,q is approximately circular complex Gaussian distributed according to the central limit theorem [28]. Moreover, as the incident signals and the additive noise are mutually independent, we can easily obtain that the expectation of εp,q is
(19)Eεp,q=ER^p,q1+ER^p,q2+ER^p,q3=0,and the variance of εp,q is
(20)Varεp,q=ER∧p,q12+ER∧p,q22+ER∧p,q32,where
(21)ER∧p,q12=1L∑i=1K∑j=1j≠iKrirj=1L∑i=1Kri2−∑i=1Kri2,ER∧p,q22=1L∑i=1Kriσ2+∑j=1Krjσ2=2L∑i=1Kriσ2,ER∧p,q32=1Lσ4.Substituting (21) into (20), we can directly have (17).DenotevecE as the vectorized estimation error matrix, and note that the distribution property of ε1,1 is not defined in Proposition 1; then, one has
(22)EJvecE22=∑p=1M∑q=1,p+q−2>0MEεp,q2=M2−1Varε2,1,where Varε2,1 is obtained by (20) and J=0M2,IM2−1 is the selection matrix by which Eε1,12 is dropped from E. Hence, according to Proposition 1 and (22), we can define the function FR^xy,Rc as
(23)FR^xy,Rc=JvecR^xy−Rc22,and then, β can be easily set as
(24)β=μM2−1Varε2,1,where μ is a user-specific weighting factor permitting (24) to be held in a high probability. It is worth noting that the threshold β is influenced by various factors, such as the variances of signals and noise, limited snapshots, and array geometry. In contrast, μ just introduces a scale to the threshold. Thus, μ has much smaller dynamic than β and is easier to choose, which is verified in the Numerical Results. In the following, an approximation of Varε2,1 is given.Firstly,∑i=1Kri+σ2 can be well estimated by averaging the diagonal elements of the observed covariance matrix along the x and y axes. It is worth noting that since the first diagonal element is common among the covariance matrices, it should be calculated only once. Hence, ∑i=1Kri+σ2 can be estimated by
(25)12M−1diagR∧x+diagR∧yT1M−R^y1,1,where R^x and R^y denote the sample covariance matrix along the x and y axes estimated from the L collected snapshots, respectively. In practical implementation, the diagonal elements of R^x can be calculated as diagR^x=diag1/L∑t=1Ldiagxtdiag∗xt, as well as diagR^y. R^y1,1 is the 1,1th element of R^y. Next, according to (4), the p,qth (p+q−2>0) element of Rxy as Rp,q can be expressed as Rp,q=∑i=1Kriej2πp−1fi,1e−j2πq−1fi,2. Thus, the expectation of the squared modulus of Rp,q is
(26)ERp,q2=∑i=1Kri2+∑i=1K∑i′=1i′≠iKriri′×cos2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2.The termcos2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2 varies for different p,q,i,i′. Moreover, the phase mod2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2,2π is uniformly distributed in −π,π with respect to uniform distributions of fi,1 and fi,2 within −0.5,0.5. Thus, the cosine term has a mean of zero. Further, if we take the average of the module of R^p,q,p+q−2>0, the cosine term will be eliminated and only the first term is retained. Hence, we can approximate ∑i=1Kri2 with the average of the module of R^p,q,p+q−2>0 as
(27)∑i=1Kri2≈JvecR^xy22M2−1.Finally, substituting (25) and (27) into (17), the approximation of Varε2,1 is obtained, and then, β can be determined via (24). With the obtained β and the function FR^xy,Rc in (26), we have the proposed CCMR-based DANM (CCMR-DANM) method for L-shaped array as follows:
(28)u~1,u~2,R~c=argminu1,u2,Rc12MTrTu1+TrTu2s.t.JvecR^xy−Rc22≤βTRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0.Compared with the SDP formulation of standard DANM in (9), beside the error tolerance constraint, an extra constraint in (14) is included in the proposed SDP formulation. It is an effectiveness relaxation of the positive constraint r≥0, which is introduced by the prior information of the crosscorrelation matrix. Moreover, the importance of this extra constraint is verified in simulations.


### 3.5. DOA Retrieval

By solving (28), we note that besides the estimation of Rc, the estimation of Tu1 and Tu1 also are obtained, which can be expressed in the form of an augmented matrix as
(29)RA=Tu1RcRcHTu2.Next, we develop three kinds of angle estimation methods based onTu1 and Tu2, Rc, and RA, respectively.First, sinceTu1 and Tu2 contain the unknown angle information, the Vandermonde decomposition-based methods such as MUSIC and ESPRIT [17, 29, 30] can be employed for angle estimation in each dimension. Of course, a pairing operation needs to be done to finally obtain the angle pairs [18].Second, consideringRc as the noise-free crosscorrelation matrix, conventional 2D DOA estimation algorithms based on crosscorrelation [18–15] can be grafted with the proposed CCMR-DANM method for angle estimation. In this paper, we adopt the JSVD algorithm [10] as a representation for comparison in simulations.Last, according to the decoupled atomic norm theory [18, 25], we have
(30)RA=Tu1RcRcHTu2=AxAy1MRsAxAyH=HRs′HH,where H=AxT,AyTT and Rs′=1/MRs. Hence, the conventional covariance-based 2D DOA estimation algorithms, such as [8, 9], can be incorporated into the proposed CCMR-DANM method.Remark 2.
Note that all the matrices we obtained through (28) are low rank and their rank values are equal to the source numbers. In the three kinds of proposed angle estimation methods, either eigenvalue decomposition (EVD) or singular value decomposition (SVD) is taken on these low-rank matrices. Hence, we can determine the source numbers by counting the number of eigenvalues or singular values, which is larger than a predefined threshold, e.g., 0.05λmax , where λmax is the maximum value of the eigenvalues. In this sense, the proposed methods can be done in a blind mode without knowing the source numbers a priori.Remark 3.
It is worth noting that the three kinds of proposed angle estimation methods have different behaviors in terms of estimation accuracy and computational complexity. The first Vandermonde decomposition-based method is computationally efficient but requires an extra pairing operation. The thirdRA-based method is relatively computational expensive than the other two. The second crosscorrelation-based method is not only computationally efficient but also can achieve automatic pairing. Moreover, note that although the proposed method works on sample crosscorrelation matrix, the proposed method grafting with the RA-based method can detect 2M−1 sources, which is larger than that of the conventional crosscorrelation-based DOA estimation methods.In summary, the proposed CCMR-DANM method is widely applicable in practical implementations.


## 3.1. Prior Art: DANM Technique

We now review the standard DANM technique for harmonic retrieval. Define a matrixZ as
(5)Z=∑k=1Kαkaxfk,1ayHfk,2,where αk∈ℂ. Then, based on (5), a matrix-form atom set of infinite size is defined as [18, 15]
(6)A=axf1ayHf2∣f1,f2∈−12,12.Accordingly, the atomic norm ofZ over the atom set A is defined as
(7)∥Z∥A=inf∑kαk∣∑kαkaxfk,1ayHfk,2,axfk,1ayHfk,2∈A;αk∈ℂ,which seeks the sparsest (under l1-norm measure) decomposition of Z over A. Consider the matrix Z is contaminated by the noise matrix N and the matrix at hand is Z^=Z+N. Then, according to the DANM theory [18, 25], the atomic decomposition yields the true structure in (5), through the following DANM formulation:
(8)Z~=argminZZAs.t.Z^−ZF2≤η,where η is a user-specified parameter for error tolerance. Moreover, (8) is equivalent to the following semidefinite positive (SDP) formulation
(9)z~1,z~2,Z~=argminz1,z2,Z12MTrTz1+TrTz2s.t.Z^−ZF2≤ηTz1ZZHTz2≽0.


## 3.2. Standard DANM-Based Formulation for CCMR

Note thatRc=∑i=1Kriaxfi,1ayHfi,2. Apparently, Rc has a sparse linear atomic representation over the matrix-form atom set in (6). And we introduce a new matrix-form atomic norm
(10)RcA+=inf∑iri∣∑iriaxfi,1ayHfi,2,axfi,1ayHfi,2∈A;ri≥0,∀i.Note that this norm differs from the original atomic norm of the standard DANM in (7), because of the extra constraint r≥0.GivenRc, it is possible to retrieve the components ri,fi,1,fi,2 of its sparest representation by calculating its atomic norm. Considering the obtained at hand is the sample crosscorrelation containing the estimation error, it boils down to
(11)R~c=argminRcRcA+s.t.FR^xy,Rc≤β.where FR^xy,Rc is the function quantifying the estimation error and β indicates the error tolerance threshold. It is similar to the DANM formulation introduced in (8), but defined on the norm ·A+ instead of ·A. To show the intricacy of this difference, we rewrite (11) in the following equivalent form:
(12a)R~c,r~=minRc,rRcA=r1(12b)s.t.FR^xy,Rc≤β(12c)Rc=AxdiagrAyH(12d)r≥0.Note that (12c) is implicit in the objective function (11) but becomes an explicit constraint because of the new nonnegative constraint (12d). Without (12d), the SDP implementation of the DANM in (9) can be used to reformulate (12) into a convex problem. However, because of the extra constraint on r in (12d), this problem becomes intractable, because r is intertwined with the other variable Rc in the form of (12c).


## 3.3. Effective Relaxation

To solve (12) in a tractable manner, we seek to relaxr≥0 to an effective form with respect to Rc. To this end, we note that
(13)Rc:,1=∑i=1Kriaxfi,1,Rc1,:H=∑i=1Kriayfi,2,where Rc:,1 and Rc1,: denote the first column and the first row of Rc, respectively. Moreover, according to the property of Toeplitz matrices [26], if r≥0, we have
(14)TRc:,1≽0,TRc1,:H≽0.Adopting (14) to replace (12d) and reformulating (12a)–(12c) into the original decoupled SDP form in (9), we reach the following effective SDP relaxation for (11):
(15)u~1,u~2,R~c=argminu1,u2,Rc12MTrTu1+TrTu2s.t.FR^xy,Rc≤βTRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0.With an appropriate definition ofFR^xy,Rc and β, (18) can be solved successfully via off-the-shelf convex solvers, such as CVX [27].


## 3.4. Estimation Error Constraint

To define a niceFR^xy,Rc which leads to easy setting of β, we denote the estimation error matrix as
(16)E=R^xy−Rc.Denote byεp,q the p,qth element of the estimation error matrix E; then, one has the following proposition:Proposition 1.
For adequately largeL, and ∀p+q−2>0, εp,q is approximately circular complex Gaussian distributed with zero mean, and the variance is
(17)Varεp,q=1L∑i=1Kri+σ22−∑i=1Kri2.p+q−2>0,1≤p,q≤M.Proof.
Note thatR^xy is estimated from the L collected snapshots. Denote R^p,q as the p,qth element of R^xy. xpt, yqt, sit, nxpt, and nyqt are similarly defined. We have ∀p+q−2>0,
(18)R^p,q=1L∑t=1Lxptyq∗t=1L∑t=1L∑i=1Ksite−j2πp−1fi,1+nxpt∑j=1Ksj∗tej2πq−1fj,2+nyq∗t=∑k=1K1L∑t=1Lskt2e−j2πp−1fk,1ej2πq−1fk,2+1L∑t=1L∑i=1K∑j=1j≠iKsitsj∗te−j2πp−1fi,1ej2πq−1fj,2+1L∑t=1L∑i=1Ksitnyq∗te−j2πp−1fi,1+∑j=1Ksj∗tnxptej2πq−1fj,2+1L∑t=1Lnxptnyq∗t=R^p,q0+R^p,q1+R^p,q2+R^p,q3,where R^p,q0, R^p,q1, R^p,q2, and R^p,q3 denote the first, the second, the third, and the fourth summand. Note that ER^p,q0=∑i=1Krie−j2πp−1fi,1ej2πq−1fi,2=Rcp,q, where Rcp,q is the p,qth element of Rc. Hence, the p,qth estimation error εp,q can be appropriately estimated as εp,q=R^p,q1+R^p,q2+R^p,q3. Since L is sufficiently large, εp,q is approximately circular complex Gaussian distributed according to the central limit theorem [28]. Moreover, as the incident signals and the additive noise are mutually independent, we can easily obtain that the expectation of εp,q is
(19)Eεp,q=ER^p,q1+ER^p,q2+ER^p,q3=0,and the variance of εp,q is
(20)Varεp,q=ER∧p,q12+ER∧p,q22+ER∧p,q32,where
(21)ER∧p,q12=1L∑i=1K∑j=1j≠iKrirj=1L∑i=1Kri2−∑i=1Kri2,ER∧p,q22=1L∑i=1Kriσ2+∑j=1Krjσ2=2L∑i=1Kriσ2,ER∧p,q32=1Lσ4.Substituting (21) into (20), we can directly have (17).DenotevecE as the vectorized estimation error matrix, and note that the distribution property of ε1,1 is not defined in Proposition 1; then, one has
(22)EJvecE22=∑p=1M∑q=1,p+q−2>0MEεp,q2=M2−1Varε2,1,where Varε2,1 is obtained by (20) and J=0M2,IM2−1 is the selection matrix by which Eε1,12 is dropped from E. Hence, according to Proposition 1 and (22), we can define the function FR^xy,Rc as
(23)FR^xy,Rc=JvecR^xy−Rc22,and then, β can be easily set as
(24)β=μM2−1Varε2,1,where μ is a user-specific weighting factor permitting (24) to be held in a high probability. It is worth noting that the threshold β is influenced by various factors, such as the variances of signals and noise, limited snapshots, and array geometry. In contrast, μ just introduces a scale to the threshold. Thus, μ has much smaller dynamic than β and is easier to choose, which is verified in the Numerical Results. In the following, an approximation of Varε2,1 is given.Firstly,∑i=1Kri+σ2 can be well estimated by averaging the diagonal elements of the observed covariance matrix along the x and y axes. It is worth noting that since the first diagonal element is common among the covariance matrices, it should be calculated only once. Hence, ∑i=1Kri+σ2 can be estimated by
(25)12M−1diagR∧x+diagR∧yT1M−R^y1,1,where R^x and R^y denote the sample covariance matrix along the x and y axes estimated from the L collected snapshots, respectively. In practical implementation, the diagonal elements of R^x can be calculated as diagR^x=diag1/L∑t=1Ldiagxtdiag∗xt, as well as diagR^y. R^y1,1 is the 1,1th element of R^y. Next, according to (4), the p,qth (p+q−2>0) element of Rxy as Rp,q can be expressed as Rp,q=∑i=1Kriej2πp−1fi,1e−j2πq−1fi,2. Thus, the expectation of the squared modulus of Rp,q is
(26)ERp,q2=∑i=1Kri2+∑i=1K∑i′=1i′≠iKriri′×cos2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2.The termcos2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2 varies for different p,q,i,i′. Moreover, the phase mod2πp−1fi,1−fi′1+2πq−1fi′,2−fi,2,2π is uniformly distributed in −π,π with respect to uniform distributions of fi,1 and fi,2 within −0.5,0.5. Thus, the cosine term has a mean of zero. Further, if we take the average of the module of R^p,q,p+q−2>0, the cosine term will be eliminated and only the first term is retained. Hence, we can approximate ∑i=1Kri2 with the average of the module of R^p,q,p+q−2>0 as
(27)∑i=1Kri2≈JvecR^xy22M2−1.Finally, substituting (25) and (27) into (17), the approximation of Varε2,1 is obtained, and then, β can be determined via (24). With the obtained β and the function FR^xy,Rc in (26), we have the proposed CCMR-based DANM (CCMR-DANM) method for L-shaped array as follows:
(28)u~1,u~2,R~c=argminu1,u2,Rc12MTrTu1+TrTu2s.t.JvecR^xy−Rc22≤βTRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0.Compared with the SDP formulation of standard DANM in (9), beside the error tolerance constraint, an extra constraint in (14) is included in the proposed SDP formulation. It is an effectiveness relaxation of the positive constraint r≥0, which is introduced by the prior information of the crosscorrelation matrix. Moreover, the importance of this extra constraint is verified in simulations.


## 3.5. DOA Retrieval

By solving (28), we note that besides the estimation of Rc, the estimation of Tu1 and Tu1 also are obtained, which can be expressed in the form of an augmented matrix as
(29)RA=Tu1RcRcHTu2.Next, we develop three kinds of angle estimation methods based onTu1 and Tu2, Rc, and RA, respectively.First, sinceTu1 and Tu2 contain the unknown angle information, the Vandermonde decomposition-based methods such as MUSIC and ESPRIT [17, 29, 30] can be employed for angle estimation in each dimension. Of course, a pairing operation needs to be done to finally obtain the angle pairs [18].Second, consideringRc as the noise-free crosscorrelation matrix, conventional 2D DOA estimation algorithms based on crosscorrelation [18–15] can be grafted with the proposed CCMR-DANM method for angle estimation. In this paper, we adopt the JSVD algorithm [10] as a representation for comparison in simulations.Last, according to the decoupled atomic norm theory [18, 25], we have
(30)RA=Tu1RcRcHTu2=AxAy1MRsAxAyH=HRs′HH,where H=AxT,AyTT and Rs′=1/MRs. Hence, the conventional covariance-based 2D DOA estimation algorithms, such as [8, 9], can be incorporated into the proposed CCMR-DANM method.Remark 2.
Note that all the matrices we obtained through (28) are low rank and their rank values are equal to the source numbers. In the three kinds of proposed angle estimation methods, either eigenvalue decomposition (EVD) or singular value decomposition (SVD) is taken on these low-rank matrices. Hence, we can determine the source numbers by counting the number of eigenvalues or singular values, which is larger than a predefined threshold, e.g., 0.05λmax , where λmax is the maximum value of the eigenvalues. In this sense, the proposed methods can be done in a blind mode without knowing the source numbers a priori.Remark 3.
It is worth noting that the three kinds of proposed angle estimation methods have different behaviors in terms of estimation accuracy and computational complexity. The first Vandermonde decomposition-based method is computationally efficient but requires an extra pairing operation. The thirdRA-based method is relatively computational expensive than the other two. The second crosscorrelation-based method is not only computationally efficient but also can achieve automatic pairing. Moreover, note that although the proposed method works on sample crosscorrelation matrix, the proposed method grafting with the RA-based method can detect 2M−1 sources, which is larger than that of the conventional crosscorrelation-based DOA estimation methods.In summary, the proposed CCMR-DANM method is widely applicable in practical implementations.


## 4. Discussions

Two related issues regarding the proposed solution in Section3 are discussed in this section. First, the computational complexity of the proposed method is analyzed and compared with that of existing CCMR-based methods. Second, we extend the proposed method to the sparse L-shaped array cases.


### 4.1. Computational Complexity

Note that the CC-ANM, CC-GLS, and the proposed method are all two-phase CCMR-based methods. They have common operations in sample crosscorrelation matrix construction and DOA estimation. The only difference is the technique used for CCMR. Hence, we only compare the computational complexities of the CCMR operation in these three methods. According to Vandenberghe and Boyd [31], the computational complexity for solving the SDP formulation of (28) is OM2+2M24M2.5. Moreover, the computational complexities of CCMR in the proposed CCMR-DANM, CC-ANM, and CC-GLS are listed in Table 1.Table 1
The computational complexities of CCMR in different algorithms.
AlgorithmComputational complexityCC-ANMO3M2−2M2M4.5CC-GLSO2M2−2M2M4.5ProposedOM2+2M24M2.5Further, Figure2 presents the computational complexities of CCMR in different algorithms versus the number of antennas at each subarray. Obviously, the proposed CCMR-DANM has the least computational cost compared with CC-ANM and CC-GLS when M>4. Further, the computational gaps between the proposed and the other two algorithms become large with the growth of the number of antennas. Hence, the proposed method offers us a more efficient choice for CCMR, especially with large M.Figure 2
The computational complexities of CCMR in different algorithms versusM.


### 4.2. Extension to Sparse L-Shaped Array Cases

In the sparse L-shaped array case, we consider there are two sparse linear array (SLA) arranged along thex and y axes. Moreover, each SLA is a subarray of a virtual ULA. In other words, the manifolds of the SLAs termed Ax′ and Ay′ can be expressed as [22]
(31)Ax′=IΩxAx,Ay′=IΩyAy,where IΩx is the selection matrix for the SLA in the axis with Ωx being the index set indicating the index of remaining antennas and IΩy is similarly defined. Hence, we have
(32)Rx′y′=IΩxRc+N1IΩyH,and sequentially, the error estimation matrix can be defined as
(33)E′=R^x′y′−IΩxRcIΩy.Moreover, in this paper, we consider that the two SLAs are two redundancy linear arrays [22], which means the virtual crosscorrelation as Rxy can be fully determined by Rx′y′. Further, each element in E′ has similar distribution as that in E. In other words, we can similarly defined the function FR^x′y′,IΩxRcIΩy and the tolerant bound β′ as those in (23) and (24), respectively. Then, we have the proposed CCMR-DANM for sparse L-shaped array as
(34)u~1,u~2,R~c=minu1,u2,Rc12MTrTu1+TrTu2s.t.J′vecR^x′y′−IΩy⊗IΩxvecRc22≤β′TRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0,where J′ is similarly defined as J.


## 4.1. Computational Complexity

Note that the CC-ANM, CC-GLS, and the proposed method are all two-phase CCMR-based methods. They have common operations in sample crosscorrelation matrix construction and DOA estimation. The only difference is the technique used for CCMR. Hence, we only compare the computational complexities of the CCMR operation in these three methods. According to Vandenberghe and Boyd [31], the computational complexity for solving the SDP formulation of (28) is OM2+2M24M2.5. Moreover, the computational complexities of CCMR in the proposed CCMR-DANM, CC-ANM, and CC-GLS are listed in Table 1.Table 1
The computational complexities of CCMR in different algorithms.
AlgorithmComputational complexityCC-ANMO3M2−2M2M4.5CC-GLSO2M2−2M2M4.5ProposedOM2+2M24M2.5Further, Figure2 presents the computational complexities of CCMR in different algorithms versus the number of antennas at each subarray. Obviously, the proposed CCMR-DANM has the least computational cost compared with CC-ANM and CC-GLS when M>4. Further, the computational gaps between the proposed and the other two algorithms become large with the growth of the number of antennas. Hence, the proposed method offers us a more efficient choice for CCMR, especially with large M.Figure 2
The computational complexities of CCMR in different algorithms versusM.


## 4.2. Extension to Sparse L-Shaped Array Cases

In the sparse L-shaped array case, we consider there are two sparse linear array (SLA) arranged along thex and y axes. Moreover, each SLA is a subarray of a virtual ULA. In other words, the manifolds of the SLAs termed Ax′ and Ay′ can be expressed as [22]
(31)Ax′=IΩxAx,Ay′=IΩyAy,where IΩx is the selection matrix for the SLA in the axis with Ωx being the index set indicating the index of remaining antennas and IΩy is similarly defined. Hence, we have
(32)Rx′y′=IΩxRc+N1IΩyH,and sequentially, the error estimation matrix can be defined as
(33)E′=R^x′y′−IΩxRcIΩy.Moreover, in this paper, we consider that the two SLAs are two redundancy linear arrays [22], which means the virtual crosscorrelation as Rxy can be fully determined by Rx′y′. Further, each element in E′ has similar distribution as that in E. In other words, we can similarly defined the function FR^x′y′,IΩxRcIΩy and the tolerant bound β′ as those in (23) and (24), respectively. Then, we have the proposed CCMR-DANM for sparse L-shaped array as
(34)u~1,u~2,R~c=minu1,u2,Rc12MTrTu1+TrTu2s.t.J′vecR^x′y′−IΩy⊗IΩxvecRc22≤β′TRc:,1≽0,TRc1,:H≽0Tu1RcRcHTu2≽0,where J′ is similarly defined as J.


## 5. Numerical Results

In this section, we present numerical examples to verify the effectiveness of the proposed method and evaluate the performance of the proposed method for L-shaped array DOA estimation. Unless specifically stated, in simulations, we consider three source signals, i.e.,K=3, with DOAs (15°, 10°), (25°, 20°), and (35°, 30°), impinge onto the L-shaped array with M=8. The number of collected snapshots L and Monte Carlo trials N is set to 200 and 500, respectively. We use root mean square error (RMSE) to evaluate DOA estimation precision of the proposed method with comparison to the MUSIC algorithm for L-shaped array (replaced with MUSIC hereafter), the JSVD [10], the CESA [12], the CC-GLS [22], and the CRB [9]. We omit the CC-ANM for comparison since it has similar estimation performance as CC-GLS and it is more computationally expensive [22]. The RMSE of the azimuth and elevation angle estimations is defined as
(35)RMSE=1K∑i=1K1N∑n=1Nϕ∧i,n−ϕi2+θ∧i,n−θi2,where ϕ^i,n and θ^i,n denote estimates of ϕi and θi from the nth trial, respectively.First, to guarantee the proposed CCMR-DANM method works in a best status, it is important to choose an appropriate user-specific weighting factorμ. The RMSEs of the proposed method with different μ versus SNR are presented in Figure 3. The curves in Figure 3 indicate that with different μ, the resulting proposed methods have different estimation performance. Moreover, the proposed method equipped with μ=1 has the best performance. Hence, in following simulations, we set μ=1 to implement the proposed CCMR-DANM method.Figure 3
RMSE of the proposed method with differentμ versus SNR.Next, let us evaluate the estimation performance of the proposed CCMR-DANM method compared with the aforementioned methods. Moreover, the proposed method without the constraint (14), termed proposed w/o 2, is also simulated for comparison. Figure 4 shows the RMSEs of these methods versus SNR. As shown in Figure 4, the proposed method is effective and outperforms both the conventional L-shaped DOA estimation methods and the existing gridless CCMR-based methods. Moreover, there is a nonnegligible gap between the proposed method and the proposed w/o 2 method. It indicates that the performance of the proposed method degrades greatly when the constraint in (14) is dropped, which means the constraint in (14) is indispensable for the proposed method.Figure 4
RMSE versus SNR withL=300.To further examine the performance, we present the RMSEs of aforementioned methods versus the number of snapshots withSNR=0 dB in Figure 5. The result likewise indicates that the proposed method has better estimation accuracy than the conventional methods. Combining the estimation performance comparisons with the analysis in Remark 2, we can conclude that the proposed method not only has the better estimation performance compared with the conventional L-shaped DOA estimation methods but also can be implemented without knowing the source numbers. Moreover, combining the estimation performance comparisons with the analysis of computational complexity in Section 4.1, we can conclude that the proposed method not only is more efficient but also has better estimation performance, compared with the existing gridless CCMR-based methods, such as CC-GLS and CC-ANM.Figure 5
RMSE versus the number of snapshots withSNR=0 dB.Last, let us consider the sparse L-shaped array cases. Assume both of the SLAs are sampled from the8 element ULA with the selection matrix IΩ, where Ω=1,2,4,7,8. And the same three source signals with DOAs (15°, 10°), (25°, 20°), and (35°, 30°) impinge onto the sparse L-shaped array. As shown in Figure 6, the 2D DOAs can be clearly observed via the proposed method. In other words, the proposed method is applicable to the sparse L-shaped array cases.Figure 6
Angle estimation via the proposed CCMR-DANM in a sparse L-shaped array case withL=200, SNR=5 dB.


## 6. Conclusion

In this paper, a new CCMR-based two-phase method for 2D DOA estimation with L-shaped array via DANM is proposed. In the first phase, the DANM technique simultaneously considering the noise and finite snapshot effects and structure property is employed to reconstruct the noise-free crosscorrelation matrix from the sample counterpart. Then, in the second phase, conventional 2D DOA estimation methods for L-shaped array can be adopted for 2D-DOA estimation. Numerical simulations demonstrate the effectiveness and outperformance of the proposed method.



---

*Source: 6659679-2021-03-30.xml*