Gridless Underdetermined Direction of Arrival Estimation in Sparse Circular Array Using Inverse Beamspace Transformation

Underdetermined DOA estimation, which means estimating more sources than sensors, is a challenging problem in the array signal processing community. This paper proposes a novel algorithm that extends the underdetermined DOA estimation in a Sparse Circular Array (SCA). We formulate this problem as a matrix completion problem. Meanwhile, we propose an inverse beamspace transformation combined with the Gridless SPICE (GLS) algorithm to complete the covariance matrix sampled by SCA. The DOAs are then obtained by solving a polynomial equation with using the Root-MUSIC algorithm. The proposed algorithm is named GSCA. Monte-Carlo simulations are performed to evaluate the GSCA algorithm, the spatial spectrum plots and RMSE curves demonstrated that the GSCA algorithm can give reasonable results of underdetermined DOA estimation in SCA. Meanwhile, the performance of the algorithm under various configurations of SCA is also evaluated. Numerical results indicated that the GSCA algorithm can provide access to solve the DOA estimation problem in Uniform Circular Array (UCA) when random sensor failures occur.


Introduction
Estimating the directions of arrival (DOAs) of signals with a sensor array is one of the critical research topics in the field of array signal processing, and plays an indispensable role in applications such as underwater object detection with sonar [1], radar target localization [2], wireless communications [3] and radio astronomy systems [4], and so forth.
In particular, the performance of the DOA estimation system is significantly constrained by the number of sensors. For example, classical high-resolution DOA algorithms including MUSIC [5], ESPRIT [6], and maximum likelihood (ML) type methods [7] can distinguish up to N − 1 sources with N sensors. However, increasing the number of sensors will also increase the system complexity and cost. Furthermore, occasional sensor failures in array [8,9] severely degrade the DOA estimation performance. In this complex scenario, the number of sources may be greater than or equal to the number of sensors, resulting in a problem of underdetermined DOA estimation.
A sparse array combined with a sparse recovery algorithm offers a novel perspective on solving this intractable underdetermined DOA estimation problem [10,11]. Notably, array configurations play an important role in the DOA estimation system. Various kinds of sparse arrays have been studied intensively under the framework of the sparse recovery algorithm, such as the coprime array [12], nested array [13] and super nested array [14], and so forth. However, sensors are located along a line in the above-mentioned array configurations, which discards the scenarios of sensors being placed in a two-dimensional plane.
Among various planar arrays, the Uniform Circular Array (UCA) has attracted much attention due to its advantages of covering a 360 • azimuthal field of view and being easy to conform to cylindrical structures. Therefore, the UCA is widely adopted in systems such as Massive MIMO [15], Ground-Based Radar [16], and the tracking of Unmanned Aerial Vehicles (UAVs) [17], and so forth. Yet the UCA degenerates into a Sparse Circular Array (SCA) when sensor failures occur. The underdetermined DOA estimation in SCA is still an open problem. Following the idea of the nested sparse linear array (NSLA), a DOA estimation algorithm in the nested sparse circular array (NSCA) has been proposed in [18]; it is worth noting that a L 1 norm based sparse recovery method is adopted in [18][19][20]. Apparently, in this regime where a pre-defined dense grid dictionary of array steering vectors evaluated at the concerned DOA angle range is necessary, a grid mismatch problem is caused. Moreover, hyperparameters introduced in the above algorithms significantly impact the performance, and the selection of those hyperparameters is not mentioned in those papers [18,21,22].
Obviously, a gridless hyperparameter-free method for underdetermined DOA estimation in SCA is urgently needed. The idea of gridless sparse recovery, firstly introduced in [23], has received great attention from the spectral estimation community. Various prominent algorithms have been developed such as atomic norm minimization (ANM) [24,25], enhanced matrix completion (EMaC) [26], and the covariance fitting type method named gridless SPICE (GLS) [27], and so forth. The GLS has been adopted intensively among the above due to its outstanding ability to cope with multiple measurement vectors (MMV) and a hyperparameter-free property.
Nevertheless, GLS is based on finding a Toeplitz or Hankel structured matrix to fit the sample covariance matrix and interpolates the missing samples simultaneously, which cannot be satisfied in SCA. The non-Vandermonde-structured steering vector of SCA creates a big obstacle to the application of the GLS algorithm. Fortunately, a method named beam space transformation (BT) which transforms the steering vector of the UCA into a virtual Vandermonde-structured steering vector has been proposed in [28] and utilized in [20,29]. Inspired by this, we propose a gridless hyperparameter-free algorithm based on the inverse beamspace (IBT) transformation of the Toeplitz matrix which offers a way to adopt the GLS algorithm in the SCA . The sample covariance matrix of SCA in element space is completed to a Toeplitz matrix in beam space, which also provides convenience for the application of the efficient Root-MUSIC algorithm [2]. Computer simulations are performed which demonstrate the ability of the proposed algorithm to handle the tricky underdetermined DOA estimation problem in SCA. We summarize the differences and connections between our work and other related works in Table 1. The main contributions of this work are summarized as follows: 1.
We propose an inverse beam space transformation (IBT) of the Toeplitz matrix in SCA scenario. The missing elements in the sample covariance of SCA are completed; 2.
A gridless hyperparameter-free algorithm is proposed to cope with the underdetermined DOA estimation problem in SCA. The efficient outstanding Root-MUSIC method based on the completed covariance matrix can be adopted; 3.
Numerical simulations are performed under various scenarios to evaluate the proposed GSCA (Gridless DOA Estimation in Sparse Circular Array) algorithm.
Boldface lowercase letters such as a, b denote vectors, and boldface uppercase letters such as A, B denote matrices, and [A] i,j denotes the (i, j)-th component of matrix A. I N is the N × N identity matrix. ∠z means taking the argument of the complex number z.
The remainder of this paper is organized as follows. Section 2.1 describes the signal model of SCA, and the GLS algorithm is introduced in Section 2.2. The inverse beamspace transformation is introduced in Section 2.3. The proposed algorithm is introduced in Section 3. The simulation results and related discussions are included in Section 4. Finally, Section 5 concludes this paper.

Signal Model for SCA
As shown in Figure 1, we consider an SCA that is composed of N p physical sensors selected from a N-element UCA with radius R. The n-th angle coordinate of the element located on the UCA is given by: Let Ω be the coordinate index set of N p integers selected from integers {1, 2, · · · , N}, and the angle coordinates generated by Ω are represented as follows: where [Ω] i is the i-th smallest number in the set Ω. Assume D far-field narrowband sources with azimuthal DOAs φ = [φ 1 , φ 2 , · · · , φ D ] T impinging on the SCA. The k-th observed snapshot is modeled as: where s(k) = [s 1 (k), s 2 (k), · · · , s D (k)] T is the source signal vector, and n(k) ∼ CN (0, σ 2 n I N ) is the additive white Gaussian noise vector. A p (φ) corresponds to the manifold matrix of SCA [18] which is formulated as: where a p (φ d ) is the steering vector of the SCA, and the i-th element is given by: where R = R/λ is the radius normalized by wavelength. Moreover, the K observed snapshot vectors can be packaged into a matrix as Furthermore, the sample covariance matrix of SCA in element space is calculated as:

Covariance Matrix Recovery with GLS
In this subsection, we briefly review the GLS algorithm which is the underlying framework of our algorithm. The GLS algorithm proposed in [27] is a gridless extension of the sparse iterative covariance-based estimation (SPICE) [31] method. The core idea of SPICE is to perform DOA estimation based on covariance fitting. The cost function of covariance fitting is given as: and where R = 1 K YY H is the observed sample covariance matrix, and Y ∈ C N×K is the matrix of snapshots. In the ULA regime where the array steering vector has a Vandermonde structure, the covariance matrix R can be re-parameterized as R = Toep{u}, which is given by: After a series of mathematical simplifications of (7) and (8), the semidefinite problem (SDP) [32] is casted as (K < N): Once we solve problem (10) or (11), the estimation of covariance matrix is obtained from R = Toep{u }. Meanwhile, the covariance-based DOA estimation algorithm is being adopted.
In the sparse linear array (SLA) scenario, the above SDPs are extended to the following SDPs. In the case of (K < N), the SDP is given by: when K ≥ N, the SDP is given as: where Γ Ω ∈ {0, 1} N p ×N is the selection matrix with its entries being 1 only at the [Γ Ω ] n p ,[Ω] np , Similarly, the estimated covariance matrix is obtained as R = Toep{u }, which can be regarded as a completed covariance matrix of the virtual ULA. Moreover, we are able to perform DOA estimation of up to N − 1 sources with the above covariance matrix. As we can see, the GLS algorithm offers a way to solve the underdetermined DOA estimation problem.
However, the Toeplitz structured covariance matrix is satisfied by ULA or virtual ULA, which is an essential precondition of the GLS algorithm. However, in the scenarios of UCA or SCA, the non-Vandermonde structured steering vector creates a big obstacle for the application of GLS algorithm. Inspired by the BT method, we extend the GLS algorithm into the SCA scenario by IBT. The IBT method is introduced in the following subsection.

Inverse Beamspace Transformation (IBT) of SCA
The beamspace transformation method is presented in [28], which provides a general way to reformulate the DOA estimation problem with UCA into virtual ULA. Let M denote the highest order mode that can be excited on a circle of normalized radius R at a reasonable strength, which is given as: where · is the round-down operator. The m-th, m ≤ |M| phase mode is excited by the normalized beamforming vector in terms of The resulting UCA far-field beam pattern of mode m is where p = cN − m and q = cN + m [29]. In order to make the first item j m J m (2π R)e jmφ in (16) be the dominant one, the number of antennas N needs to meet the following condition: By using the property J −m (2π R) = (−1) m J m (2π R) of Bessel functions, and the residual terms are being omitted, the UCA beam pattern for mode m can be expressed as: For brevity, we define the following matrix in terms of: and where N B = 2M + 1 is the number of beam [20]. The Vandermonde structured array steering vector in beamspace is defined as: By introducing C J , (16) can be represented as: which is Obviously we have: To sum up, the relation between a(φ) and a v (φ) is given as: where T B is defined as Apparently, T B transforms the steering vectors in beamspace into element space, which is exactly the reverse of the original beamspace transformation. Therefore, the above process is named Inverse Beamspace Transformation (IBT).
As we can see, (25) offers great convenience for handling the non-Vandermonde structured steering vector of UCA which is the basic framework of our algorithm.

Proposed Algorithm GSCA
The proposed algorithm named GSCA (Gridless DOA Estimation in Sparse Circular Array) is summarized in Algorithm 1. In order to visualize the principle of the algorithm, we draw the main steps of the algorithm in Figure 2.  Next, the GSCA algorithm is introduced in detail. The GSCA algorithm is mainly based on the GLS algorithm combined with the aforementioned IBT. Notably, the elements of SCA are selected from a UCA, with utilizing the predefined selection matrix Γ Ω , the relation between the physical snapshot of SCA and the complete snapshot of UCA is established as: Thus, the covariance matrix of SCA is written as: where R = E{yy H } = E{(As + n)(As + n) H } with E{sn H } = E{s H n} = 0 being used; meanwhile, R s and R n are the signal covariance matrix and noise covariance matrix, respectively. By taking advantage of the IBT, matrix R x can be reparameterized as: Similarly, with an application of the GLS algorithm, the SDPs arising in the SCA scenario are shown below. In the case of (K < N B ), the SDP is given as: When K ≥ N B , the SDP is given as: Once problem (31) or (32) is solved, the completed covariance matrix of the UCA in element space is obtained as R = T B Toep{u }T H B . The middle part of R is a Toeplitz structured covariance which can be regarded as a beamspace transformed covariance matrix of the UCA. Next, we focus on the middle part which is marked as R B = Toep{u }. Apparently, a classical Root-MUSIC algorithm [33] can be performed thanks to the Toeplitz structure of R B .
The eigenvalue decomposition (EVD) of R B is given as: where the signal subspace U s , noise subspace, U n and corresponding eigenvalues Λ s , Λ n have the following forms: As we all know, the noise subspace U n is orthogonal to the signal subspace U s , and U s spans the same subspace as the steering matrix which is written as It is obvious to formulate the following equation: The null spectrum is formed as: for z = e jφ , and p(z) is defined as: Moreover, the null spectrum is also able to reformulate into a polynomial as follows: where h i is calculated by: The D roots inside the unit circle with the largest magnitude are chosen.Then the DOAs are obtained by:

Simulation Results
In this section, computer simulations are carried out to demonstrate the performance of the proposed DOA estimation algorithms. The root-mean-square error (RMSE) is adopted, which is defined as: where P = 200 is the number of Monte Carlo trials. Additionally, the spatial spectrum is depicted in the polar coordinate to visualize the performance of DOA estimation. To obtain the spatial spectrum, we replace step 5 of the GSCA algorithm in Algorithm 1 (Root-MUSIC [33]) with a spatial spectrum search (MUSIC [5]). The spatial search step is 0.1 • , and the search range is [0 • , 360 • ).

Selection of N, N p and Ω
In this subsection, the simulation results are presented to illustrate the selection of N, N p , and Ω. We explored the effect of various N p on the performance of DOA estimation under a selected N. The SNR and K are chosen as 15 dB and 1024, and the number of Monte Carlo trials P is 200. The number of sources D and the number of physical elements N p are equal in order to satisfy the underdetermined scenario. Obviously, multiple label sets Ω will be generated under each pair of (N, N p ). In order to exclude the influence of the particularity of Ω on the results, label set Ω is randomly generated in each trial. We set a threshold of RMSE (43) to evaluate the simulation results. The simulation results are shown in parts a-c of Table 2, respectively. (The minimum N p that succeeds under each N is bolded.)  From the above results, we can roughly draw the following empirical conclusions to select N p , which is given by: where · is the round-up operator. To sum up, we choose the minimum N p that succeeds when N = 7 or 9. Apparently, when sensor failure occurs, the locations of the faulty sensors are random. Considering the circular symmetry, there are three different array configurations when N = 7 and N p = 5; and the number of different array configurations is seven when N = 9 and N p = 6. We have plotted these arrays in Figures 3 and 4, respectively. We perform the following simulations to evaluate the performance of the GSCA algorithm under different array configurations when the SNR and the number of snapshots K are fixed at a moderate value (SNR = 15 dB, K = 1024). The number of sources D is an integer variable selected from N p , N to satisfy the underdetermined scenario. The simulation results are shown in Figure 5a and 5b, respectively. Obviously, when the simulation parameters {SNR, K, D, N, N p } are the same, the RMSEs of different array configurations are almost at the same level, which means the robustness of our method to various array configurations with the same N and N p . Therefore, in the following simulations, we selected one type of array from each group to study the effects of SNR and K.  In this simulation, we study the performance of DOA estimation versus SNR and K. The number of sensors in UCA: N is set to 7, and the number of physical sensors N p is set to 5 with label Ω = {1, 2, 3, 5, 7} (Type 2 in Figure 3). A schematic of this SCA is shown in Figure 3b. The normalized radius R is set to 0.65, in this case, N B is 7 based on (14). Notably, the source number D is set to 5 and 6 in order to perform underdetermined DOA estimation, and DOAs are set to be equidistantly distributed in [0 • , 360 • ) [34].
The RMSE versus SNR and number of snapshots K are plotted in Figure 6, respectively. As shown in Figure 6a, the RMSE is gradually dropping as the SNR increases; yet the RMSE is slightly dropping as K increases. In addition, the RMSE increases as the number of sources D increases when the SNR and K are fixed (Figure 6b). Furthermore, Figure 7a,b depicts the normalized spatial spectrums under D = 5 and D = 6, respectively. The SNR is set to be 10 dB, and the number of snapshots K is 1024. Apparently, as the number of sources D increases, the number of outlier peaks increases when SNR and the number of snapshots K are fixed.  In this simulation, the number of sensors in UCA, N, is set to 9, and the number of physical sensors N p is set to 6. The index vector of the physical sensor is Ω = {1, 2, 3, 5, 7, 9} (Type 5 in Figure 3). The normalized radius R is set to 0.7, and N B is calculated as 9 with (14). The array structure is shown in Figure 4e. In contrast with the simulation of N = 7, N p = 5, the RMSE curves and the normalized spatial spectrum are drawn in Figures 8 and 9, respectively. Moreover, the RMSE curve of D = 6 in Figure 8a converges to 1.2 • while the RMSE curve of D = 6 in Figure 6a converges to 2.8 • . Furthermore, comparing Figures 7b and 9a, it is evident that the number of outlier peaks decreases when D is fixed as 6. Those benefits come with the enlarged array aperture R and the increased number of physical sensors N p .

Complexity Analysis
The major computational complexity of the proposed GSCA algorithm corresponds to the step (31) or (32). A well-known off-the-shelf SDP solver SDPT3 [35] is employed to solve our algorithm. The SDPT3 is an interior-point-based method that has the computational complexity of O(n 2 1 n 2.5 2 ) [36], where n 1 denotes the number of variables and n 2 × n 2 is the dimension of the positive semi-definite matrix in the SDP. In our cases, n 1 = N B + N 2 p and n 2 = 2N p , such that the computational complexity of solving SDP is O((N B + N 2 p ) 2 (2N p ) 2.5 ). In addition, the EVD (34) step also contributes a large part of the computational complexity, which is O(N 3 B ). The complexity of polynomial rooting steps (37)-(41) is O(DN B ). Thus the major computational complexity of the proposed GSCA algorithm is . We evaluate the algorithm under CPU I7-10510U at 2.30 GHz and 12 GB RAM. The average CPU running time of 200 Monte-Carlo trials is given in parts a and b in Table 3, respectively.

Conclusions
In this paper, we have proposed the GSCA algorithm to perform underdetermined DOA estimation in SCA. The GSCA algorithm takes advantage of the inverse beamspace transformation (IBT), together with the GLS algorithm; in this way, the covariance matrix of SCA is completed to a Toeplitz matrix in beamspace; meanwhile, the Root-MUSIC is adopted and DOAs are obtained. We have performed computer simulations, and results demonstrate that the proposed algorithm is able to produce reasonable results of underdetermined DOA estimation in SCA. Furthermore, the GSCA algorithm still works well in various array configurations, which means the tricky DOA estimation problem in UCA with random sensor failures can be handled. In the future, we will work on improving the performance of the algorithm and strive to extend our algorithm to the two-dimensional DOA estimation problem.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: