Gridless sparsity-based DOA estimation for sparse linear array

: A novel two-stage gridless sparsity-based method for sparse linear array direction-of-arrival (DOA) estimation is presented. First, the covariance matrix gridless sparse representation method based on atomic norm minimisation is proposed for corresponding structured Toeplitz matrix construction and source number detection. Then, the conventional MUSIC algorithm can be employed for DOA estimation. Compared with conventional subspace-based algorithms, the proposed two-stage method can be carried out without knowing source number and directly detect more source signals than sensors. Numerical simulations demonstrate the effectiveness and outperformance of the proposed method.


Introduction
Direction-of-arrival (DOA) estimation as an important role in array signal processing has attracted much interest in the last several decades. Conventional direction finding algorithms as MUSIC [1], ESPRIT [2], and ML [3] have been widely applied. However, the disadvantages cannot be ignored, e.g. requirement of source number, sensitiveness to noise. Moreover, consider applying to sparse linear array (SLA) case, they cannot detect more source signals than sensors. Although, the spatial smoothing (SS) technique [4] can be employed to extend the aperture, it can only handle the full augmentable sparse array case. Recent years, inspired by compressive sensing and sparse representation, several sparsity-based methods are presented for SLA DOA estimation [5][6][7]. However, the continuous angle space is reduced to a set of discrete grids in these methods to formulate the sparse representation, where the model error as basis mismatch is unavoidable. As a result, their performance degrade in practical applications [8]. More recently, to overcome this inherent drawback, several gridless sparsity-based methods spring up for SLA DOA estimation. Based on low-rank matrix reconstruction (LRMR), [9,10] proposes the covariance matrix reconstruction approach (CMRA) using nuclear (trace) norm minimisation for DOA estimation. Using atomic norm minimisation, the atomic norm soft thresholding for multiple measurement vectors (AST for MMV) method is presented in [11], while it is worked in the data domain. What's more, by utilising atomic norm minimisation in the correlation domain, we propose the denoising covariance matrixbased atomic norm minimisation (DCMANM) method for DOA estimation [12]. However, as source signals are considered to be correlated, the atomic norm for two-dimensional harmonic is applied which results in computation expensive.
In this paper, only uncorrelated source signals are considered, then the one-dimensional atomic norm is using for sparse representation formulation in the correlation domain. The resulting computational efficient covariance matrix gridless sparse representation (CMGSR) method is proposed for corresponding structured Toeplitz matrix construction and source number detection. Whereafter, the DOAs are determined from the corresponding Toeplitz matrix with the help of conventional MUSIC method. The whole proposed two-stage method as CMGSR + MUSIC can be carried out without knowing source number and directly detect more source signals than sensors.

Model formulation
Consider an SLA that consists of M omnidirectional sensors with locations [0, d 1 , …, d M − 1 ], where d m − 1 is the distance between the reference sensor and the mth sensor. In this paper, we assume d m − 1 is an integral multiple of half a wavelength, and d M − 1 = N(λ/2), where λ denotes the wavelength. Consequently, the SLA can be regard as a subset of a virtual N -elements ULA with half a wavelength spaced and the sensor index set of the SLA is defined as Ω ⊂ 1, …, N with Ω = M.
Assume K narrowband signals impinge onto the SLA from directions of θ = [θ 1 , …, θ K ]. Under the aforementioned assumptions, the output of the SLA is formulated as where A Ω θ = a Ω θ 1 , …, a Ω θ K is the manifold matrix of the SLA with a Ω θ i = 1, e j2π d 1 /λ sinθ i , …, e j2π(d M − 1 /λ)sinθ i T being the corresponding steering vector of ith source signal. A = a θ 1 , …, a θ K denotes the manifold matrix of the virtual ULA with a θ i = 1, e j2π 1/2 sinθ i , …, e j2π(N /2)sinθ i T . Γ Ω is a selection matrix with only the Ω i th position at the ith row being 1 and 0 otherwise. In this paper, source signals are uncorrelated and uncorrelated with the Gaussian White noises. Hence, the covariance matrix of the observed signals can be expressed as where p = p 1 , …, p K T denotes the source power parameter, σ is the noise variance parameter.

CMGSR
By vectorisation operation, one has J. Eng where ⊕ denotes the Khatri-Rao product. What's more, Hence, r Ω can be reformulated as where the second equality applies the fact that p ≥ 0 which means every element in p is no less than zero. W = Γ Ω ⊗ Γ Ω G, y = Ap and y N − 1 is the last N − 1 rows of y.
Defined a set of atoms as With a f = 1, e j2π f , …, e j2π(N − 1) f T . It is straightforward to find that y is a linear combination of elements from the atomic set A by setting f = (1/2)sin θ . To characterise the sparsity of y, it is following from [13] that the atomic norm of y is defined as where ϕ k = 0 is a prior information. To proceed, we give an equivalent formulation of ϕ k = 0 and we have the following lemma. Toeplitz matrix whose first column is equal to u.T u ≥ 0 indicates that it is a semi-definite positive matrix. Moreover, the representation is unique.
The proof is omitted here due to the space limitation and will be detailed in a journal version.
By minimising the atomic norm with eliminating the noise item (σvec(I M )) and prior information constraint, it seems to be easy to reconstruct y from r Ω . However, in practical applications, R Ω can be only estimated from the finite L snapshots and the estimation error is inevitable between R Ω and the sample covariance matrix R Ω , and so is between r Ω and r^Ω. Consequently, to reconstruct y from r^Ω, a sparse reconstruction formulation, in the form of minimising the atomic norm with estimation error constraint, is given as (see (9)) where J = blkdiag E 1 , …, E M denotes the denoising matrix with E m = e 1 , …, e m − 1 , e m + 1 , …, e M T .e m is the standard orthogonal basis of one at the mth position. β is an userdefined estimation error thresholds. In what follows, we give a computation of β.
Note that where Λ = R Ω − R Ω is the estimation error matrix and Λ ii = diag diag Λ denotes the diagonal matrix with the main diagonal elements of Λ. Denoted ε ij by the elements of Λ at the i, j th position. It follows from Proposition 1 in [5] that for sufficient large snapshots and i > j, ε ij is approximately circular complex Gaussian distributed with zero mean, and the variance is What's more, since R Ω and R Ω are Hermitian matrices, ε ji = ε ij * . Hence, we can easily extend Proposition 1 to the following corollary.
Corollary 1: For sufficient large snapshots, ε ij is approximately circular complex Gaussian distributed with zero mean, and the variance is This corollary indicates that Consequently, the threshold can be given as where μ is a user-specific weighting factor. What's more, the estimate of Var ε 2, 1 can be easily obtained through the process proposed in [5]. We refer interested readers to related paper for more details. With the determined β and rewriting ∥ y ∥ A into its equivalent semi-definite programming (SDP) characterisation [14], we have the following SDP formulation min T, y, t where T is a Hermitian Toeplitz matrix. Equation (15) can be solved efficiently using off-the-shelf SDP solver as CVX [15]. As the atomic set A is a continuous dictionary with infinite atoms from the view of sparse representation, to distinguish from conventional sparsity-based methods, we designate this proposed method as CMGSR method.

DOA estimation
With the CMGSR method carried out, not only y but also T are obtained. Moreover, it is following from [14] that T can be represented as T = ∑ i = 1 K p k a f k a H f k . Since CMGSR is carried out without the knowledge of source number, here with the obtained T, the number of source signals can be easily determined as the rank of T (the number of eigenvalues which are larger than a predefined threshold), and the conventional direction finding method as MUSIC can be employed together with T to retrieve the DOAs.
In conclusion, the whole proposed method for SLA DOA estimation has two stage as CMGSR + MUSIC. The CMGSR method is carried out for corresponding Toeplitz matrix construction and source number detection. Whereafter, the MUSIC method is employ with the corresponding Toeplitz matrix for DOAs retrieval. It is worth noting that the proposed CMGSR algorithm is flexible since other conventional direction finding methods like ESPRIT, ML, can be straightforwardly incorporated into it. In this paper, we focus on the MUSIC algorithm for illustration.

Related discussion
As the DOAs are retrieved from the N × N Toeplitz matrix, the maximum source number or equally the maximum rank of this Toeplitz matrix for unique decomposition is N − 1 based on the Vandermonde Decomposition Lemma [16]. In other words, the proposed method is able to retrieve more number of source signals than that of sensors in this SLA case.
According to Vandenberghe and Boyd [17], the worst computational complexity of CMGSR is O N 4.5 . As the MUSIC method has a complexity of O N 3 , the whole complexity of proposed method is O N 4.5 + N 3 . This conclusion shows that computational complexity of the proposed method stays much lower than that of the DCMANM, which is O N 8.5 . Although the proposed method has much higher computational complexity than the conventional subspace-based methods, it can be carried out without knowing source number and can directly detects more source signals than sensors.

Numerical simulations
In this section, we take several numerical experiments to demonstrate the effectiveness and outperformance of the proposed method. We use root mean square error (RMSE) to evaluate the DOA estimation accuracy. The proposed CMGSR method is implemented with μ = 1.
First, suppose that two uncorrelated source signals with different power 1, 2 impinge onto a four-element SLA from −10.1°, 10.1°. The sensor index set Ω = 1, 2, 5, 7 , i.e. it is a minimum-redundancy array (MRA), a special full augmentable SLA. The number of the collected snapshots is set to 200. SS-MUSIC [1,4], CMRA and the proposed method are carried out to estimate the DOAs, the simulation results are obtained through a total of 300 trails. The CRB [18] is also taken into account for comparison. The RMSEs of DOA estimation versus SNR are plotted in Fig. 1. As expected, the proposed method has the best performance though the overall SNR.
What's more, Fig. 2 depicts the RMSEs of DOA estimation versus the number of collected snapshots. The SNR is fixed at 0 dB, the number of the collected snapshots varies from 50 to 250 and 300 trials are simulated. As shown in Fig. 2, the proposed method exceeds all the other methods in DOA estimation precision and can obtain good DOA estimates even with as few as 50 snapshots.
Finally, consider a SLA with sensor index set Ω = 1, 2, 6, 8 , we attempt to employ the proposed method to estimate the DOAs and powers of seven source signals impinged onto this SLA from directions as −55°, − 30°, − 15°, 0°, 15°, 30°, 55° at SNR = 10 dB. As this SLA is not full augmentable, the conventional SS-MUSIC cannot detect those signals. However, as is shown in Fig. 3, the proposed method detects them well. The results are derived from 300 trials. Black circles are the positions of the true DOAs and corresponding powers, and green dots denote the estimated ones.

Conclusion
In this paper, a gridless sparsity-based method named CMGSR + MUSIC is presented for SLA DOA estimation. First, the CMGSR algorithm is proposed for corresponding Toeplitz matrix construction and source number detection. Then, the classical MUSIC algorithm can be employed for DOA estimation. The proposed method can be carried out without knowing source number and can directly detect more source number than sensors. Numerical simulations demonstrate the effectiveness and outperformance of the proposed method. Worthy of noting is that the computational complexity of proposed method equipped with CVX grows fast with the increase of antenna aperture. To obtain a computational efficient implementation of the proposed method, the alternating direction method of multipliers (ADMM) [19] needs to be studied for solving the proposed convex problem, which is posed as a future work.