Interpolation-Based Direction-of-Arrival Estimation for Coprime Arrays via Covariance Matrix Fitting

Coprime arrays can highly increase degree-of-freedom (DOF) by exploiting the equivalent virtual signal. However, since the corresponding virtual array constructed by the coprime array is always a non-uniform linear array (non-ULA), most existing direction-of-arrival (DOA) estimation algorithms fail to utilize all received information and result in performance degradation. To address this issue, we propose a novel interpolation approach for coprime arrays to convert the virtual array into a ULA with which all received information can be efficiently utilized. In this paper, we consider a weighted covariance matrix fitting criterion to formulate a semi-definite programming (SDP) problem with respect to the interpolated virtual signal. After that, we can reconstruct a Hermitian Toeplitz covariance matrix corresponding to the interpolated ULA in a gridless manner, and the number of detectable targets is ulteriorly increased with the reconstructed covariance matrix. The proposed approach is hyperparameter-free so that the tedious process of selecting regularization parameters is avoided. Numerical experiments validate the superiority of the proposed interpolation-based DOA estimation algorithm in terms of DOF characteristic, resolution ability and estimation accuracy compared with several existing techniques.


I. INTRODUCTION
Direction-of-arrival estimation, as an important research branch of array signal processing, has been widely applied in sonar [1], radar [2] and wireless communication [3]. If the common array geometry like ULA or UCA is under consideration, the number of identifiable sources is impossible to exceed the array scale according to the Nyquist sampling theorem, i.e., an array consisting of M physical sensors can only detect up to M −1 targets with the subspace-based algorithms [4], [5]. To break through such a limitation, a systematical sparse array called the coprime array [6], [7] was proposed recently and has attracted noticeable attention due to its superior performance. Comprising two sparse ULAs which are deployed in accordance with coprimality, the coprime array forms a larger aperture compared with the common array configurations which consist of the same number of The associate editor coordinating the review of this manuscript and approving it for publication was Wei Liu . sensors. Thus the resolution ability of the coprime array is relatively improved. Besides, compared with other sparse array configurations like minimum redundancy array (MRA) [8] and minimum hole array (MHA) [9], the sensor positions of the coprime array can be certainly determined as long as the coprime integer pair is given. Furthermore, by exploiting the second-order information of the received signal to construct an augmented virtual array, the DOF of the coprime array can reach up to O(MN ) with only O(M + N ) physical sensors.
To make full use of the DOF offered by the coprime array, almost all existing DOA estimation approaches are performed on the derived virtual signal. However, since the coprime array is partially augmented, holes are generated in the difference co-array, i.e., the corresponding virtual array is discontinuous. Some DOA estimation algorithms, such as spatial smoothing MUSIC [10], co-array ESPRIT [11], and covariance matrix sparse reconstruction [12], are only performed on the maximum continuous segment within the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ virtual array, and hence the utilization of DOF is greatly reduced. Capable of using the total discontinuous elements, sparse signal reconstruction (SSR) algorithm [13] can automatically find the DOAs via solving an optimization problem but suffers from basis mismatch that targets do not fall on the pre-defined spatial grid exactly. To address such a problem, the idea of joint-sparsity is used in [14] to correct the offgrid DOAs with deviation variables. Besides, Tan et al. [15] design a sparsity recovery problem from the perspective of super-resolution, and all possible sources are covered within a continuous range. The methods proposed in [14] and [15] can also operate the derived signal corresponding to the nonuniform virtual array. Actually, the holes in the non-uniform virtual array can be regarded as potential sensors, and all information received by the coprime array can be successfully utilized once the holes are filled up. Consequently, numerous array interpolation approaches are proposed to convert the discontinuous virtual array into a ULA. Based on the principle of matrix completion, the virtual array is interpolated through nuclear norm minimization in [16]. The OptSpace [17], a kind of matrix completion algorithm, is also utilized in the virtual array domain to enhance DOF. Whereas the matrix completion-based approaches will generate a performance loss because the partial correlation observation corresponding to the non-uniform segment is retained. To address such a problem, the idea of gridless Toeplitz matrix reconstruction [18], [19] is utilized to exploit the potential DOF more recently and it is demonstrated that the reconstructed covariance matrix coincides with the interpolated ULA better. The performance degradation caused by the non-uniform array configuration is nonnegligible, especially in complex electromagnetic and multi-target scenarios. From the perspective of Toeplitz matrix reconstruction in this paper as well, we propose a novel virtual array interpolation-based DOA estimation algorithm for the coprime array. We consider a weighted covariance matrix fitting criterion to formulate an SDP optimization problem with respect to the interpolated virtual signal, and then a Hermitian Toeplitz covariance matrix corresponding to the interpolated ULA is reconstructed in a gridless manner. After that, the available DOF exceeds the one of the original virtual array, and numerous ULA-based DOA estimation algorithms can be readily applied. Unlike the approaches in [16]- [19], the proposed approach dispenses with the tedious process of selecting regularization parameters and it is proved to be equivalent to promoting low-rankness of the covariance matrix via atomic norm minimization implicitly. The contributions of this paper are summarized as follows: • We derive the atomic norm of the multiple virtual measurements in both continuous and discontinuous cases, and build up the relationship between them.
• We utilize a covariance matrix fitting criterion to reconstruct a Hermitian Toeplitz matrix corresponding to the interpolated virtual array in a gridless manner.
• We explore the implication of the proposed method to promote low-rankness and analyze its feasibility of virtual array interpolation.
The rest of this paper is organized as follows. In section II, the signal model of the coprime array is presented. We then detail the proposed virtual array interpolation-based DOA estimation algorithm in section III. The effectiveness of the proposed algorithm is demonstrated via simulation results in section IV and the conclusions are finally drawn in section V.
Notations: Throughout this paper, we use lower-case and upper-case boldface characters to represent vectors and matrices respectively. The superscripts (·) T , (·) * and (·) H denote the transpose, complex conjugation and conjugate transpose, respectively. The notation vec(·) represents the vectorization operator that sequentially stacks each column of a matrix. We use ⊗ to stands for the Kronecker product and | · | represents the cardinality of a set.

II. COPRIME ARRAY SIGNAL MODEL
As Fig. 1 shows, the coprime array consists of two sparse ULAs which are located at {0, Md, 2Md, · · · (N −1)Md} and {0, Nd, 2Nd, · · · , (M − 1)Nd} respectively. Herein, M and N are a pair of coprime integers, and the unit inter-element spacing d is chosen to be a half-wavelength λ/2. Owing to the coprimality, only the first sensors of these two layers are co-located, and thus the coprime array comprises M + N − 1 sensors in total. With the assumption that L far-field uncorrelated narrowband signals from θ = [θ 1 , θ 2 , · · · , θ L ] T impinge on the coprime array, the received signal at the t-th moment can be modeled as where A(θ ) = [a(θ 1 ), a(θ 2 ), · · · , a(θ L )] ∈ C (M +N −1)×L represents the manifold matrix of the coprime array, s(t) = [s 1 (t), s 2 (t), · · · , s L (t)] T ∈ C L×1 denotes the signal waveform vector, and n(t) ∼ CN (0, σ 2 n I) denotes the additive complex-valued zero-mean Gaussian white noise. Suppose that the location of the i-th sensor is u i , and a(θ ) is the steering vector with the i-th element taken as e −j(2π/λ)u i sinθ .
To explore the second-order information of the received signal, we express the covariance matrix among L sources as where p l denotes the power of the l-th source. Limited by the number of snapshots, the expectation operator is unavailable in practice, and it is usually replaced by the statistical average derived from K available snapshots. After that, the equivalent virtual signal can be formulated by vectorizing the covariance matrix as .Â is the manifold matrix of the virtual array with sensor locations determined by the following difference co-array After removing the identical elements within S, an equivalent virtual array is constructed as the subset of S Such an augmented virtual array is shown in Fig. 2(b) where M = 3 and N = 5 are set for illustration. By selecting the corresponding elements from y, the equivalent virtual signal of S V is expressed as whereÂ v ∈ C |S V |×L represents the manifold matrix of the derived virtual array andî denotes the subset formed by the corresponding elements of i. Note that S V is discontinuous since the coprime array is partially augmented, i.e., there exist holes within S V . To make the utmost of all information received by the coprime array, we intend to convert S V into a continuous virtual ULA and denote it as S I . The purpose of virtual array interpolation is to recover the unknown signal datas yielded by the potential sensors which are illustrated as the hollow circles in Fig. 2(c).

III. THE PROPOSED VIRTUAL ARRAY INTERPOLATION APPROACH FOR DOA ESTIMATION
In this section, a virtual array interpolation approach for the coprime array is proposed to further increase DOF.
Herein, we first derive the atomic norm of the multiple virtual measurements in both continuous and discontinuous cases. We then consider a weighted covariance matrix fitting criterion to reconstruct a Hermitian Toeplitz covariance matrix corresponding to the interpolated ULA S I . Moreover, we explore the implication of the proposed approach to promote low-rankness and analyze its feasibility of virtual array interpolation. Finally, a ULA-based DOA estimation algorithm is applied to identifying much more sources with the reconstructed covariance matrix.

A. ATOMIC NORM EXPLOITATION
To avoid the basis mismatch problem, we introduce a mathematical tool known as the atomic norm to characterize the sparsity of the virtual statistics. To elaborate the property of the atomic norm, we first consider a noise-free and snapshotinfinite scenario where the potential sensors are accurately interpolated as well. Besides, we define J = (| S V | +1)/2 and I = (| S I | +1)/2. Note that the virtual signal behaves in a single-snapshot manner. Based on the idea of spatial smoothing, the ideal interpolated virtual signal denoted as y I can be divided into I segments shown as Fig. 3. After stacking the virtual signal of each subarray as Y = y 1 I , y 2 I , · · · , y I I , it is equivalent to obtaining the multiple measurements yielded by the reference subarray located from 0 to (I − 1)d. Owing to the conjugate symmetry of y I , the multiple measurements can be also formulated as T y 1 I where T (·) represents the operator to form a Hermitian Toeplitz matrix with the inner vector as its first column. Such a Toeplitz-structure matrix can be directly regarded as the covariance matrix corresponding to the reference subarray [20] for the reason that the virtual signal contains signal power information instead of waveform information. Furthermore, it is evident that the available DOF is up to I − 1. Under the ideal assumption and DOF constraint, it is readily derived that rank T y 1 I = L < I . According to the Vandermonde decomposition theorem [21], T y 1 I can be uniquely decomposed as VOLUME 8, 2020 where a I (θ ) denotes the steering vector of the interpolated reference array with the i-th element taken as e −j(i−1)πsinθ and P = diag(p). We formally define a atom set as Obviously from (7), Y is a linear combination of L atoms in the atomic set A. To promote sparsity, we consider the atomic 0 norm as the sparse metric, and it is defined as the minimum number of atoms in A that can synthesize Y where inf denotes the infimum. Similar to the 0 norm, the minimization of atomic 0 norm is also an NP-hard problem. Hence, we present the following atomic norm convex relaxation [22] defined as the gauge function of the convex hull of A (denoted as conv(A)) Still considering the noise-free and snapshot-infinite scenario, we find that the signal y v yielded by the discontinuous virtual array S V is the observed samples of y I on an index subset , i.e., y I = y v where · i denotes the i-th element of the inner vector. We define a selection matrix ∈ R J ×I which is constructed by removing the i-th row (i / ∈ ) of an I × I unity matrix. Then the multiple measurements of S V can be expressed as where A = A I . Similarly we define a atom set of the discontinuous multiple measurements as So that the atomic norm of Y can be readily derived as Lemma 1: proof: Through the definition given in (10) and (13), the following equalities evidently hold According to Lemma 1, it is reasonable to recover the unknown statistics of potential sensors even though only partial information is observed, i.e., the discontinuous virtual array S V can be theoretically interpolated into a ULA S I .

B. VIRTUAL ARRAY INTERPOLATION VIA COVARIANCE MATRIX FITTING
In practice, we can only obtain the derived virtual signal formulated as (6) and the multiple virtual measurements expressed as (11). Through spatial smoothing, the covariance matrix of reference discontinuous array can be calculated as where the coefficient 1/J is omitted because it makes no difference after normalization. Due to the impact of noise, finite snapshots and discontinuity of y v , the deviation is accumulated in R v when performing (16). According to the property of the array signal, the covariance matrix corresponding to S V should hold the structure as R = A PA + σ 2 n I = T y 1 I H + σ 2 n I. (17) The noise term σ 2 n I will not change the special structure of R so that we can also reparameterize R as T (z) H . After that, T (y 1 I ) can be reasonably modeled as T (z) − λ min (T (z))I where λ min (·) denotes the minimum eigenvalue of the inner matrix.
To reduce the deviation caused by R v , we consider the following covariance matrix fitting criterion [23]- [25] for the purpose of covariance matrix reconstruction where · F denotes the Frobenius norm. Unlike [24] and [25] where singular matrix should be considered to make the criterion universal, the reference covariance matrix is nonsingular and positive semi-definite (PSD) when it turns to the virtual array domain, i.e., R −1/2 v definitely exists even though only one snapshot is captured. After a simple calculation, the covariance matrix fitting criterion (18) is reformulated as where tr(·) denotes the trace operator. Therefore, under the constrains T y 1 I ≥ 0 and σ 2 n ≥ 0, the covariance matrix fitting problem can be formulated as the following convex optimization problem to minimize (18) min The covariance matrix fitting problem (20) is an SDP that can be efficiently solved by the CVX software [26]. After solving (20), we can obtain the interpolated virtual signal y 1 I as well as the power of the additive noise.

C. PERFORMANCE ANALYSIS AND DOA ESTIMATION
The noise-free covariance matrix theoretically behaves in a low-rank manner because the incident sources should be less than the sensors within S I . Different from the approaches proposed in [16]- [19] which promote low-rankness via nuclear norm or trace norm minimization explicitly, the proposed method enforces low-rankness via atomic norm minimization in an implicit way. By substituting (16) to the first term of (20), (20) is equivalent to After a simple transformation, (21) can be reformulated as (22) where z = z and T (z ) = T (z) . (22) can be further transformed as According to the theorem 2 mentioned in [18], (23) is equivalent to the following atomic norm minimization problem By Lemma 1, (24) can be ulteriorly expressed as Therefore, the proposed covariance matrix fitting problem formulated in (20) is capable of recovering the unknown signal of the potential virtual sensors, i.e., S V can be effectively interpolated as S I by (20). Moreover, the proposed method can promote the low-rankness of the covariance matrix via atomic norm minimization in an implicit manner. The solution to (20) can be utilized to construct a noise-free PSD Hermitian Toeplitz matrix T y 1 I which corresponds to the interpolated virtual ULA. Hence, numerous proven ULA-based DOA estimation algorithms, including MUSICbased [10], [20], ESPRIT-based [11], [27] and sparsity-based [13]- [15] approaches can be efficiently applied to the virtual array domain for unambiguous DOA estimation. For example, herein we formulate the MUSIC spatial spectrum as where n denotes the noise subspace of T y 1 I . Actually T y 1 I is approximating to a rank-full matrix in the noisy finite-snapshot cases, and thus the noise subspace n is generally acquired by selecting the eigenvectors corresponding to the I − L smallest eigenvalues of T y 1 I . The number of targets L is regarded as a priori under the virtual domain DOF constraint L < I or can be effectively estimated via AIC, BIC [28] and SORTE [29]. The final DOA estimation result can be empirically obtained by collecting L dominant peaks of f MUSIC (θ).
We summarize the proposed virtual array interpolationbased DOA estimation approach in Algorithm 1. It is worth mentioning that the proposed covariance matrix fitting problem described in (20) is hyperparameter-free, i.e., the regularization parameter specified by users via cross-validation is unnecessary. Moreover, by modifying the location-indexed subset , the proposed interpolation approach can be readily extended to other partially augmentable array configurations [30] including MRA, MHA and the modified nested array (MNA) [31], [32].

IV. SIMULATION RESULTS
In this section, we present several numerical experiments to demonstrate the advantages of the proposed interpolationbased DOA estimation algorithm for coprime arrays. We choose a pair of coprime integers M = 3 and N = 5 respectively to deploy the coprime array which consists of M + N − 1 = 7 physical sensors in total, and the sensors are located at {0, 3d, 5d, 6d, 9d, 10d, 12d}. By performing the operation (6), a discontinuous virtual array ranging from −12d to 12d with holes located at {−11d, −8d, 8d, 11d} is constructed. The proposed algorithm is compared with the sparse signal reconstruction (SSR) algorithm [13] with step size 0.1 • , the nuclear norm minimization (NNM) algorithm [16] and the atomic norm minimization (ANM) algorithm [18]. The regularization parameters for the SSR algorithm and the ANM algorithm are both set to be 0.25 empirically.

A. DOF CHARACTERISTIC
In the first numerical experiment, we test the DOF characteristic. The signal-to-noise ratio (SNR) and the number VOLUME 8, 2020   of snapshots K are set to be 0dB and 500 respectively. Herein we first consider 9 narrowband uncorrelated sources uniformly distributed in the range of [−50 • , 50 • ], and the simulation result is illustrated in Fig. 4 where the true DOAs are illustrated as the vertical red dashed lines. It is discovered from Fig. 4 that all 9 sources can be evidently identified by the four tested algorithms with only 7 physical sensors, and thus the effectiveness of coprime arrays to increase DOF is validated. Note that the available DOF of the discontinuous virtual array S V is up to 10, but the one of the interpolated ULA S I reaches up to 12. When the incident sources are increased to 12, the SSR algorithm misses several targets as illustrated in Fig. 5(a) because it is performed on S V . However, obviously all 12 sources can be correctly identified by the NNM algorithm, the ANM algorithm and the proposed algorithm as shown in Fig. 4(b), (c) and (d) respectively. It is demonstrated that the interpolation-based DOA estimation algorithms can effectively exploit the potential virtual sensors to further increase DOF.

B. RESOLUTION ABILITY
In the second experiment, we test the resolution ability by distinguishing two closer targets which are located at − θ and θ respectively. We first present a representative spatial spectrum comparison in Fig. 6 where SNR = 0dB, K = 500, and θ = 1 • . It can be obviously seen from the spectrum of SSR illustrated in Fig. 7(a) that a pseudo peak shows up at around 2.2 • and the two dominant peaks deviate from the true DOAs. The NNM algorithm distinguishes the two targets with a bad grace shown as two obtuse peaks in Fig. 6(b). In contrast, two much sharper and more correctly identified peaks are yielded by the ANM algorithm and the proposed algorithm. We further define that the two targets are correctly resolved as long as there are two peaks showing up in the spatial spectrum and the estimation error is smaller than θ . Herein  θ varies from 0.2 • to 2 • with a step size of 0.2 • , and we consider two scenarios where SNR are 0dB and 30dB respectively. 500 Monte-Carlo trials are run for each testing point to calculate the resolution probability. Comparing Fig. 7(a) and (b), we discover that the resolution ability is improved as the SNR increases no matter which algorithm is utilized. The resolution probability yielded by the SSR algorithm is the worst and hardly reaches 1 as θ increases. It is attributed to the pseudo peaks and basis mismatch commonly generated by the SSR algorithm [13]. Compared with the ANM algorithm and the proposed algorithm, the NNM algorithm requires a relatively larger angular interval to achieve the same probability. It is because the idea of matrix completion is utilized in NNM where the partial information of discontinuous virtual array S v is retained in the interpolated covariance matrix [16]. In contrast, the ANM algorithm and the proposed algorithm recover a covariance matrix corresponding to the virtual ULA S I from the perspective of matrix reconstruction where the effect of partial correlation observations is alleviated. From  Fig. 7, we find ANM and the proposed algorithm achieve the best resolution probability and their performances are almost identical.

C. ESTIMATION ACCURACY
In the last experiment, root mean square error (RMSE) is considered to evaluate the estimation accuracy of the four tested algorithms, and the RMSE criterion is defined as whereθ l,q denotes the estimated DOA of the l-th source in the q-th Monte Carlo trial. Although 7 physical sensors are deployed, we consider 9 equal-power uncorrelated sources uniformly distributed in the range of [−50 • , −50 • ], and for every testing point, we run Q = 500 iterations of Monte Carlo simulations. The average CPU time for running SSR, NNM, ANM and the proposed algorithm on Intel i7-10510U laptop are 1.26s, 1.05s, 0.86s and 1.68s respectively. According to the derivation in [18], the Cramér-Rao bound (CRB) is also plotted for reference. From the RMSE performance versus SNR illustrated in Fig. 8(a), we can discover that the RMSE of the four tested algorithms all decrease as the SNR increases. Limited by the pre-defined spatial grids, SSR yields a relatively worse performance. Moreover, the ANM algorithm and the proposed algorithm behave in similar estimation accuracy, and they outperform the NNM algorithm in any SNR cases. The performance trend of the proposed algorithm is always consistent with the CRB prediction. Meanwhile, the RMSE comparison versus the number of snapshots as shown in Fig. 8(b) also validates the superiority of the proposed algorithm over the other three algorithms. It is VOLUME 8, 2020 demonstrated that the covariance fitting criterion considered in this paper gives a better performance in the context of virtual array interpolation.

V. CONCLUSION
In order to take full advantage of all information received by the coprime array, we propose a novel interpolation-based DOA estimation algorithm from the perspective of Toeplitz matrix reconstruction in this paper. The holes within the discontinuous virtual array are treated as potential sensors, and then we consider a weighted covariance matrix fitting criterion to reconstruct an augmented covariance matrix corresponding to the interpolated ULA. Additionally, we investigate the atomic norm of multiple virtual measurements in both continuous and discontinuous cases, based on which the feasibility of the proposed algorithm to promote lowrankness and recover the unknown signal of the potential sensors is demonstrated. After virtual array interpolation, the available DOF is further increased and numerous ULAbased DOA estimation approaches can be readily applied. Different from most existing methods, the proposed approach is hyperparameter-free so that it dispenses with the tedious process of selecting regularization parameters. The superiority of the proposed algorithm is validated through simulation results compared with several existing techniques.