Improved Coarray Interpolation Algorithms with Additional Orthogonal Constraint for Cyclostationary Signals

Many modulated signals exhibit a cyclostationarity property, which can be exploited in direction-of-arrival (DOA) estimation to effectively eliminate interference and noise. In this paper, our aim is to integrate the cyclostationarity with the spatial domain and enable the algorithm to estimate more sources than sensors. However, DOA estimation with a sparse array is performed in the coarray domain and the holes within the coarray limit the usage of the complete coarray information. In order to use the complete coarray information to increase the degrees-of-freedom (DOFs), sparsity-aware-based methods and the difference coarray interpolation methods have been proposed. In this paper, the coarray interpolation technique is further explored with cyclostationary signals. Besides the difference coarray model and its corresponding Toeplitz completion formulation, we build up a sum coarray model and formulate a Hankel completion problem. In order to further improve the performance of the structured matrix completion, we define the spatial spectrum sampling operations and the derivative (conjugate) correlation subspaces, which can be exploited to construct orthogonal constraints for the autocorrelation vectors in the coarray interpolation problem. Prior knowledge of the source interval can also be incorporated into the problem. Simulation results demonstrate that the additional constraints contribute to a remarkable performance improvement.


Introduction
Direction-of-arrival (DOA) estimation has been a popular research field in array processing for several decades. This problem aims at retrieving the directional information of sources from the array of received signals, and plays an important role in a variety of practical scenarios [1][2][3]. Conventional DOA estimation methods such as multiple signal classification (MUSIC) [4] and estimation of signal parameters via rotational invariance technique (ESPRIT) [5] basically rely on the spatial properties of the signals impinging on the antenna array, while ignoring the temporal and frequency properties. However, many modulated signals used in communication, radar, and sonar systems exhibit a cyclostationarity property. These cyclostationary signals are not periodic with respect to time, but their statistical characteristics vary periodically with time. This property can be exploited to effectively eliminate interference and background noise [6][7][8][9]. The cyclostationarity-based DOA estimation methods were started by Gardner in [6,7], and the proposed cyclic MUSIC and cyclic ESPRIT algorithms used a cyclic correlation matrix (CCM) in place of the zero-lag covariance matrix. An extended cyclic MUSIC was proposed in [10]; besides the CCM, an additional conjugate cyclic correlation matrix (CCCM) was also utilized to obtain an increase in degree-of-freedom from O(N) to O(2N − 1). All the mentioned methods are proposed for uniform linear arrays (ULA), and fail to estimate the DOAs when the number of sources is larger than the DOFs of the ULA. Therefore, the main task in this paper is to integrate the cyclostationarity property with the spatial domain-which means that the second-order statistics incorporated with the cyclic frequency and the time lags are utilized, instead of the conventional covariance that only contains spatial information-and to improve the performance of DOA estimation under the circumstance that the sources are much more numerous than the sensors.
In order to estimate the DOAs of as many signals as possible under a given number of sensors, many sparse linear array structures have been proposed. The traditional one is the minimum redundancy array (MRA), which can generate consecutive virtual apertures by utilizing the coarray information [11,12]. However, the sensor locations are determined by exhaustive searching and no general expressions are available for the locations of the physical sensors in an MRA. This leads to an intensive investigation of the sparse linear arrays, such as the (super) nested arrays [13,14] and the (generalized) coprime arrays [15,16], which have a more concise and flexible geometry for sparse array configuration. These sparse arrays, consisting of two uniform linear arrays with different interelement spacing, open a new approach to array processing. They can resolve O(N 2 ) sources from merely N sensors. One distinctive advantage of the coprime array is that a negligible mutual coupling effect is introduced because the interelement spacing of the coprime array is several times longer than the half-wavelength, which makes it more attractive than the nested array. Intrinsically, these sparse arrays all intend to generate consecutive virtual apertures by exploiting the difference coarray. The arrays with a filled coarray are called fully augmentable arrays and when there are holes in their coarray, as in the case of coprime arrays, they are called partially augmentable arrays [17]. Three kinds of algorithms have been developed for these two types of augmentable arrays.
The first kind used compressive sensing (CS)-inspired l 1 norm minimization techniques [18][19][20] by assuming a sparse representation of the incoming DOAs on a prespecified discrete spatial grid. However, the CS-based methods suffer from basis mismatch effects since the true DOAs are unlikely to lie on the prespecified grid, no matter how fine it is chosen to be [21]. The second techniques are the spatial smoothing (SS) algorithms, which can construct a positive definite matrix in the difference coarray domain that contains noise subspace information. Once the spatially smoothed matrix is obtained as in [15,22] or by the direct coarray augmented method [23], the MUSIC method [22] and the ESPRIT method [24] algorithms can be used to perform the DOA estimation. An additional version of the smoothed matrix is proposed in [25], which is optimized through low-rank denoising to avoid the estimation of source numbers. The main limitation of these SS algorithms is that they can only utilize the consecutive virtual lags of the coarray; when there are holes as in the partially augmentable arrays, the lags out of the consecutive part of the coarray are ignored. In order to utilize the full coarray information, the Hermitian Toeplitz matrix completion-based array interpolation techniques were proposed [17,26]. Benefitting from the development of matrix completion theory [27,28], the simple nuclear norm minimization problem can be formulated to recover a low-rank matrix with small set of entries. Thus, the holes in the coarray can be filled via matrix completion. A few papers have taken into consideration the perturbation effect caused by the finite number of snapshots-an uncertainty is modelled into the matrix completion problem when assigning the existing values to their corresponding locations [29,30]. However, the error bound is set blindly and an improper setting may deteriorate the performance of the matrix completion. Moreover, [29] solved two sequential convex optimization problems, which was time-consuming and can be replaced by a modern structured matrix completion algorithm [31], while [30] utilized the redundant sample covariance to optimize the target low-rank matrix and did not make full use of the difference coarray information, leaving more entries in the matrix to be filled.
In this paper, we explore the coarray interpolation techniques with cyclostationary signals. Firstly, we use the cyclic correlation matrix and the conjugate cyclic correlation matrix as the second-order statistics and build up the difference coarray and the sum coarray models. Then, the spatial spectrum sampling operations and the derivative (conjugate) correlation subspaces are defined, and we demonstrate that the vectorization of the CCM and the CCCM lie in the correlation subspace and the conjugate correlation subspace, respectively. With the relationship between the vectorization of the sampling correlation matrices (the CCM and the CCCM) and their corresponding autocorrelation vectors in the coarray domains, we construct orthogonal constraints for the autocorrelation vectors. Moreover, prior knowledge of the source interval can be incorporated into the orthogonal constraints. Finally, we formulate a Toeplitz completion problem for the difference coarray model and a Hankel completion problem for the sum coarray model, and integrate the orthogonal constraints into the structured matrix completion problems. Numerical results demonstrate the superior performance of the proposed algorithms.
The main contribution of this paper can be summarized as follows: • We investigate the array interpolation techniques with the cyclostationary signals, building up the difference coarray model and the sum coarray model by utilizing the structure of the CCM and CCCM in a coprime array, and then formulate the array interpolation as a Toeplitz completion problem and a Hankel completion problem.

•
We define the spatial spectrum sampling operations, and prove that they can extract the entire power spectrum associated with the CCM and the CCCM. After vectorization of the spatial sampling operations, we show that the vectorization of the sampling correlation matrices lies in the subspaces, which are called the (conjugate) correlation subspaces.

•
We build orthogonal constraints for the autocorrelation vectors in the coarray domains according to the (conjugate) correlation subspaces. Prior knowledge of the source interval can be incorporated into the orthogonal constraints, which can help to improve the performance of the structured matrix completion remarkably.
The rest of this paper is organized as follows. In Section 2, we introduce the cyclostationary signal model and build up the difference coarray and sum coarray models based on the coprime array. In Section 3, the spatial spectrum operations and the derivative (conjugate) correlation subspaces are defined; then, the orthogonal constraints-which can incorporate the prior knowledge of the source interval-are constructed and integrated into the structured matrix completion problems. Section 4 presents the simulation results for comparison, and conclusions are drawn in Section 5.
Notation: Throughout this paper, scalars, vectors, matrices, and sets are denoted by lowercase letters, lowercase letters in boldface, uppercase letters in boldface, and letters in blackboard, respectively. The superscripts * , T, and H denote the complex conjugate, the transpose, and the complex conjugate transpose. The Moore-Penrose pseudoinverse is denoted with superscript †. vec(·) and E(·) represent the vectorization and expectation operations. The symbol ⊗ denotes the Kronecker product. |A| denotes the cardinality of a set A.
[A] i,j indicates the (i, j)th entry of A. The triangle bracket x S n represents the value corresponding to the support n ∈ S.

Cyclostationary Signal Model and the Coprime Coarray Model
In this section, we firstly introduce the cyclostationary signal model received by the coprime linear array, and then construct the coarray models based on the structure of the CCM and the CCCM.

Cyclostationary Signal Model
For vector s(t) formed by D cyclostationary signals, the CCM R α ss (τ) and the CCCM R α ss * (τ) are defined as where · t denotes infinite-time averaging with respect to t. Then, s(t) k (k = 1, 2, . . . D) is said to be cyclostationary if [R α ss (τ)] k,k or [R α ss * (τ)] k,k does not equal zero at frequency α for some lag parameter τ. The value of the cyclic frequency α is related to the carrier frequency and the baud rate of the signal-it is usually integer multiples of the carrier frequency or the baud rate. Some cyclostationary signals may contain both nonzero cyclic correlation and nonzero conjugate cyclic correlation.
Consider D narrowband sources, exhibiting second-order cyclostationarity with a common cycle α, impinging on the antenna array, which contains P physical sensors with the pth sensor located at z p d, where z p is an integer and d = λ/2 (λ being the carrier wavelength of the sources). The signal received on the pth sensor is where s k (t) is the kth zero-mean desired signal from direction θ k and n p (t) is the zero-mean additive Gaussian white noise at the pth sensor. There are two assumptions that must hold: (a) the impinging sources are mutually cyclically uncorrelated, which means R α ss (τ) and R α ss * (τ) are diagonal matrices; and (b) the sources are uncorrelated with the noise. Under assumptions (a) and (b), the cyclic correlation function and the conjugate cyclic correlation function between the pth and the qth sensor are given by R α where α(θ k ) = [a 1 (θ k ), a 2 (θ k ), . . . a P (θ k )] is the steering vector, whose pth element is defined as a p (θ k ) = e j2π( f c /c)z p d sin(θ k ) . Therefore, the P × P CCM and the CCCM can be constructed with the elements for all indices p and q in Equations (4) and (5) as follows: In practice, only a finite number of snapshots are available. Let x(t) denote the snapshots; the sampling correlation matrices, namely, the CCM and the CCCM, are obtained by

Coprime Coarray Model
In order to generate the augmented array, we firstly define the difference coarray and the sum coarray as the following: Definition 1. Let S = z p , 1 ≤ p ≤ P denote the set of sensor positions of a sparse array (normalized with respect to d); then, the difference coarray C d and the sum coarray C s are defined respectively as Then, a steering vector associated with the coarray can be formed using , and exhibit a linear relationship with the vectorization of the CCM and the CCCM given by J d and J s are the transform matrices, which are defined as in [32]: Let r Cd(s) be the signal received at a virtual sensor array with sensor positions given by C d(s) . If a consecutive virtual aperture can be formed by using the difference coarray or the sum coarray definition, we can resort to spatial smoothing techniques [15,23] to recover their ranks, due to the fact that r Cd(s) is a rank one vector. Otherwise, a partially augmentable array will be formed, and, in order to utilize the full coarray information, we intend to fill the holes within the augmented uniform linear arrays (ULA) V d(s) , defined as follows: To give an intuitive impression of these configurations of the coprime array, we give an example by assuming a coprime array with physical sensor positions at S = {0, 3, 5, 6, 9, 10, 12, 15, 20, 25}. The difference coarray and sum coarray are shown in Figure 1 with red circles and blue diamonds, respectively. Holes within the coarrays are marked with crosses.
Then, a steering vector associated with the coarray can be formed using Let d(s) r be the signal received at a virtual sensor array with sensor positions given by If a consecutive virtual aperture can be formed by using the difference coarray or the sum coarray definition, we can resort to spatial smoothing techniques [15,23] to recover their ranks, due to the fact that d(s) r is a rank one vector. Otherwise, a partially augmentable array will be formed, and, in order to utilize the full coarray information, we intend to fill the holes within the augmented uniform linear arrays (ULA) () ds , defined as follows: To give an intuitive impression of these configurations of the coprime array, we give an example by assuming a coprime array with physical sensor positions at = {0,3,5,6,9,10,12,15, 20, 25} . The difference coarray and sum coarray are shown in Figure 1 with red circles and blue diamonds, respectively. Holes within the coarrays are marked with crosses.

(Conjugate) Correlation Subspaces and the Proposed Coarray Interpolation Algorithms
The coarray interpolation techniques are proposed on the basis of matrix completion theory, and the Hermitian Toeplitz structure revealed in [23] is enforced on the matrix under estimation for regulation. In order to further improve the performance of the coarray interpolation, we explore the additional structure of the autocorrelation vectors, proposing definitions of spatial spectrum sampling operations and the derivative (conjugate) correlation subspaces, and then integrate the derivative orthogonal constraints into the coarray interpolation algorithms.

(Conjugate) Correlation Subspaces
Inspired by the compressive covariance sensing [33] that the array covariance matrix has a sparse representation in a dictionary constructed from the steering vectors, and the method for constructing orthogonal constraints by utilizing the spatial feature [34][35][36], we consider constructing a subspace in which the correlation matrix lies so that an orthogonal constraint can be constructed for the correlation matrix. Firstly, we give the definitions of the spatial spectrum sampling operations.

Definition 4.
For a linear array of P physical sensors, and D cyclostationary signals impinging on the antenna array from directions θ k , the spatial spectrum sampling operations for the CCM SS(R α xx (τ)) and the CCCM CSS(R α xx * (τ)) are defined respectively as follows: We observe that for an ULA with P physical sensors, if P is large enough, the P-points averaged inner product of the two steering vectors Proof. See Appendix A When the nonuniform linear array is a coprime array, the inner product α H (θ i )α(θ k )/P with θ i = θ k is much larger than the other values when θ i = θ k , and we can still approximate it with a Kronecker delta function. Thus, Equations (15) and (16) can be written as It can be seen from (18) and (19) that the integral parts will be nonzero only when i = k, and the final results of (18) and (19) will equal to R α xx (τ) and R α xx * (τ), respectively, which means that the spatial spectrum sampling operations can extract the entire power spectrum associated with the CCM and the CCCM. With the observations of (18) and (19), we then vectorize (15) and (16), and utilize the following property [37]: We thus obtain the following equations: Sensors 2018, 18, 219 It can be seen from (21) and (22) that the vectors (vec(R α xx (τ)), vec(R α xx * (τ))) are equal to themselves after being multiplied by the matrix (not an identical matrix) constructed by the integral with respect to θ i . We can draw a conclusion that the vectorization of the correlation matrices is in the column subspaces of these matrices. We refer to these matrices as the correlation subspace matrix M CS and the conjugate correlation subspace matrix M CCS with the following definitions: Definition 5. The correlation subspace matrix M CS and the conjugate correlation subspace matrix M CCS are The integral interval can be [−π/2, π/2] or can be set according to prior knowledge about the sources. The vectorization of the CCM and the CCCM is in the column subspaces of M CS and M CCS , respectively. Since M CS and M CCS are low-rank matrices, we can simplify the representations of the subspaces by using the eigenvectors associated with the positive eigenvalues, denoted by Q CS and Q CCS . These subspaces are dependent on the array configuration and prior knowledge about the sources. The dimensions of these subspaces are determined by the positive numbers of the eigenvalues of the (conjugate) correlation matrices. The eigenvalues of the (conjugate) correlation matrices are shown in Figure 2. We take the coprime array with the configuration illustrated in Figure 1 as an example.
It can be seen from (21) and (22) 22 PP  matrices with The  Figure 2. We take the coprime array with the configuration illustrated in Figure 1 as an example.  As shown in Figure 2, when we have no prior knowledge about the sources and the integral interval is set to [−π/2, π/2], the numbers of the positive eigenvalues of the correlation subspace matrix and the conjugate correlation subspace matrix are 43 and 35, respectively, which coincide with |C d | and |C s |. When a shorter source interval is available, we can use less eigenvectors to construct the subspaces, but the exact dimension can only be determined by numerical tests.
After the subspaces are specified, we can construct constraints for the sampling correlation matrices according to the fact that the vectorization of the (conjugate) correlation matrices is orthogonal to the complementary subspaces of Q CS and Q CCS , so we have where I − Q CS Q † CS and I − Q CCS Q † CCS are the projection matrices on the orthogonal subspaces. According to (25) and (26), the denoised correlation matrices can be obtained by projecting the sampling correlation matrices to their individual subspaces such that the estimation error caused by the finite length of the snapshots can be eliminated. Here, we emphasis the construction of constraints for the structured matrix completion problems rather than the operations for sampling correlation matrix denoising. Thus, we transform these constraints into the coarray domains according to the relationship between the autocorrelation vectors and the vectorization of the (conjugate) correlation matrices shown in (12) and (13), obtaining

Coarray Interpolation Algorithms
For the coprime array, the increased DOFs are achieved by utilizing the coarray information. When there are holes in the coarray, the lags out of the consecutive part of the coarray cannot be involved in implementing the DOA estimation. In order to reduce the estimation error and increase the DOFs, coarray interpolation algorithms are proposed which can utilize all the coarray information. It is based on the results of the direct coarray augmentation [23] method, which can obtain the same noise subspace as the spatial smoothing method [15,22] by rearranging the autocorrelation vector to form a Hermitian Toeplitz matrix. As to the holes in the coarray C d(s) , they can be filled based on the property that the matrix to be recovered has low-rank terms, related to the signal components. It is the typical form of matrix completion, which can recover a low-rank matrix via solving a nuclear norm minimization problem.
In the context of coarray interpolation for cyclostationary signals, we can formulate a Toeplitz completion problem for the difference coarray model according to the fact that the CCM has a Toeplitz structure, as shown in (6). We do not require it to be a Hermitian matrix due to the finite snapshot effect with the cyclostationary signals. The autocorrelation vector on the uniform grid V d is denoted by r Vd , which can be considered the antidiagonal elements of the Toeplitz matrix (denoted by T (r Vd )). The sampling autocorrelation vector on the nonuniform grid C d is denoted by r Cd , which is obtained from J † d vec( R α xx (τ)) according to Equation (12). In order to integrate the orthogonal constraint (27) into the Toeplitz completion problem, we start with the autocorrelation vector denoising problem: Then, we relax the minimization term with r Vd Cd − r C d = ε, where ε is the error term. Equation (29) is equivalent to Finally, we incorporate (30) into the nuclear norm minimization problem, which results in where · * is the nuclear norm, and τ ε 2 can be considered the error-penalizing term with penalty parameter τ. The aim of (31) is to recover a low-rank Toeplitz matrix with the values r Vd Cd , which is close to r Cd and is orthogonal to a subspace constructed using Q CS and J d . The main drawback of other coarray interpolation algorithms is that the error bound is set blindly when assigning r Cd to r Vd Cd ; the error may vary a lot according to the environment, and improper setting may deteriorate the matrix completion performance. Our algorithm avoids this problem by estimating the error automatically, and the existence of the error also allows r Vd Cd to satisfy the orthogonal constraint, which is constructed by exploring the additional structure of r Vd Cd and enables our algorithm to obtain a more accurate estimate. The penalty parameter τ can tradeoff the influence of the low-rank term and the error term, which can be set empirically according to the environment: When the signal to noise ratio (SNR) is high and the number of snapshots is large enough, we can set a large τ to penalize the error term; otherwise, a small value is set for τ in pursuit of a solution with a smaller nuclear norm. However, τ is set in a short interval in the simulations as it is not very sensitive to the environment. Now we turn to the sum coarray model for interpolation, the CCCM of the ULA possesses a Hankel structure according the expression in (7), which means that each ascending skew-diagonal from left to right is constant. The autocorrelation vector on the uniform gird V s is denoted by r Vs , which can be considered as the diagonal elements of the Hankel matrix (denoted as H(r Vs )). The sampling autocorrelation vector on the nonuniform grid C s is obtained by r Cs = J † s vec( R α xx * (τ)) according to Equation (13). The sum coarray interpolation problem is similar to the difference coarray interpolation case, and can be formulated as The minimization functions and the constraints in problems (31) and (32) are convex, so we resort to the cvx toolbox [38] to solve the constrained minimization problems. After T (r Vd ) or H(r Vs ) is obtained, the singular value decomposition (SVD) is then implemented to obtain the noise subspace, and subspace-based methods can be used to estimate the DOAs. Remarks:

1.
Even though the recovered matrices T (r Vd ) and H(r Vs ) can be considered as the correlation matrices of a ULA with N = (|V| + 1)/2 sensors, the DOFs cannot achieve N − 1 due to the fact that the filled lags do not provide additional information on the sources. The actual freedom is governed by the nonuniform grid C d(s) . For the coprime array illustrated in Figure 1, the numbers of the nonuniform grids for the difference coarray and the sum coarray are |C d | = 43 and |C s | = 35, respectively. Thus, for a modulated signal that has both the cyclostationarity and the conjugate cyclostationarity property, the difference coarray model-based Toeplitz completion algorithm can resolve more sources than the sum coarray model-based Hankel completion algorithm.

2.
The dimensions of the (conjugate) correlation subspaces need to be set carefully. When there is no prior knowledge of the sources, we set the angle interval to [−π/2, π/2], the dimension of the correlation subspace is chosen as |C d |, and |C s | is set as the dimension of the conjugate correlation subspace. When we have prior knowledge of the sources, we can implement the eigenvalue decomposition (EVD) of the (conjugate) correlation subspace matrices, and the dimensions are set according to the numbers of the resultant positive eigenvalues.

3.
The proposed coarray interpolation algorithms are not only suitable for the coprime array, but are also suitable for other partial augmentable arrays as it satisfies the recovery condition of structured matrix completion revealed in [31].

4.
The correlation subspace matrix (23) and the conjugate correlation subspace matrix (24) cannot be calculated analytically due to the fact they are unintegrable with respect to θ i . For a symmetric source interval [−β/2, β/2], the following numerical approximations are used instead: where the number of discrete samples is set to N = 1001 in the simulation. Even though we use all possible spatial grids to construct the (conjugate) correlation subspaces, this doesn't affect the fact that our algorithms can resolve gridless sources, which is the advantage over the sparsity-aware algorithms.
The entire process of the proposed algorithms is summarized in Algorithm 1.

Algorithm 1.
The proposed improved coarray interpolation algorithms with cyclostationary signals.

Input
The received signal vector x(t), cyclic frequency α, sources prior Output The optimized low-rank (conjugate) correlation matrix

Step 3
Construct the (conjugate) correlation subspace Q CS or Q CCS Step 4 Optimize (31) for the difference coarray or (32) for the sum coarray

Simulation Results
In our simulations, the coprime array was a sparse array with ten sensors located at S = {0, 3, 5, 6, 9, 10, 12, 15, 20, 25}, which is illustrated in Figure 1. After the holes are filled, the partial augmentable array can be considered as a ULA with N = 26 elements for both the difference coarray model and the sum coarray model. The performance of the proposed coarray interpolation algorithms for the difference coarray model (CI-DC) and the sum coarray model (CI-SC) were compared with several modern algorithms exploiting the coprime array, including the spatial smoothing MUSIC (SS-MUSIC) [22] and the sparsity-aware algorithms for the difference coarray model (SA-DC) and the sum coarray model (SA-SA) [20]. The root mean squared error (RMSE) was used to measure the performance, and is defined as , whereθ k denotes the estimated DOA of the kth source and I = 200 is the number of independent repetitions. The sources were binary phase shift keying (BPSK)-modulated including the signals of interest (SOIs) with a 4 Mb/s bit rate and interference with a 3.2 Mb/s bit rate. The proposed algorithms and the sparsity-aware algorithms exploiting cyclic statistics were under the optimal parameters (τ = 0.125 µs, α = 4 MHz). Cyclic MUSIC [9] was utilized in our algorithms after the low-rank (conjugate) correlation matrices were obtained.
In the first experiment, we tested the interference elimination capacity and the increased DOFs of the proposed algorithms. In the first scenario, D = 17 equal-power signals including 13 SOIs and 4 interference sources impinged on the array. All the sources were uniformly distributed between −60 • and 60 • . The interferences arrived from {θ 2 , θ 4 , θ 14 , θ 16 }. The true DOAs, including the SOIs and the interference sources, are illustrated by the vertical dashed lines in Figure 3a-d. Two thousand snapshots were sampled at the frequency of 32 MHz. The signal-to-noise ratio and the interference-to-noise ratio were 0 dB. The source angle interval was set to [−π/2, π/2]. The eigenvalues of M CS and M CCS are presented in Figure 2. We did not take the sparsity-aware algorithm into account due to its similar performance with our algorithm. Figure 3a-c demonstrate the performance in terms of the spatial spectrum of the SS-MUSIC, CI-DC, and CI-SC algorithms in the first scenario. It should be mentioned that the consecutive range of the difference coarray was from −17 to 17, indicating that the maximal DOF for the SS-MUSIC is 17 under the array configuration of S. It can be seen from Figure 3a that the SS-MUSIC method does not have the signal-selective capacity because it forms peaks at the angles of the interference sources and it has large pointing errors at θ 1 and θ 2 . The proposed CI-DC in Figure 3b and the CI-SC in Figure 3c, as expected, can null out the four DOAs of the interference sources and have correctly determined all the true DOAs of SOIs. In the second scenario, we increased the number of SOIs from D = 13 to D = 19. The sources were uniformly distributed between −60 • and 60 • . The SNR was 10 dB and 2000 snapshots are sampled. We tested the increased DOF of the proposed CI-DC algorithm shown in Figure 3d. The reason why we choose CI-DC is that the difference-coarray-based algorithm has more nonuniform grids than the sum coarray model; thus, it can resolve more sources than CI-SC. Figure 3d shows that the CI-DC algorithm can correctly estimate all the 19 SOIs and the achieved DOF is beyond the limit of SS-MUSIC, which can only utilize the consecutive part of the difference coarray.
Two thousand snapshots were sampled at the frequency of 32 MHz. The signal-to-noise ratio and the interference-to-noise ratio were 0 dB. The source angle interval was set to [ / 2, / 2]  Figure 2. We did not take the sparsity-aware algorithm into account due to its similar performance with our algorithm. Figure 3a-c demonstrate the performance in terms of the spatial spectrum of the SS-MUSIC, CI-DC, and CI-SC algorithms in the first scenario. It should be mentioned that the consecutive range of the difference coarray was from −17 to 17, indicating that the maximal DOF for the SS-MUSIC is 17 under the array configuration of . It can be seen from Figure 3a that the SS-MUSIC method does not have the signal-selective capacity because it forms peaks at the angles of the interference sources and it has large pointing errors at 1  and 2  . The proposed CI-DC in Figure 3b and the CI-SC in Figure 3c, as expected, can null out the four DOAs of the interference sources and have correctly determined all the true DOAs of SOIs. In the second scenario, we increased the number of SOIs from 13 D  to 19 D  . The sources were uniformly distributed between 60  and 60 . The SNR was 10 dB and 2000 snapshots are sampled. We tested the increased DOF of the proposed CI-DC algorithm shown in Figure 3d. The reason why we choose CI-DC is that the difference-coarray-based algorithm has more nonuniform grids than the sum coarray model; thus, it can resolve more sources than CI-SC. Figure 3d shows that the CI-DC algorithm can correctly estimate all the 19 SOIs and the achieved DOF is beyond the limit of SS-MUSIC, which can only utilize the consecutive part of the difference coarray. In the second experiment, we compared the RMSE versus the number of snapshots and SNR with the difference coarray model between our CI-DC algorithm (including two different source interval conditions and without the orthogonal constraint condition) and the SA-DC algorithm. Both algorithms utilize the CCM as the second-order statistics and have interference elimination capability. The prespecified grids for SA-DC are from [ 90 ,90 ]  with a sampling interval of 0.1 .
The regularization parameter for SA-DC was empirically chosen to be 0.4, and the penalty parameter  in our algorithm was set to 24 when the number of snapshots is less than 1000 in Figure 4a or the SNR is less than 12 dB in Figure 4b; otherwise,  was set to 26. Five equal-power BPSK SOIs uniformly distributed between 9  and 9 impinged on the array. In Figure 4a, the In the second experiment, we compared the RMSE versus the number of snapshots and SNR with the difference coarray model between our CI-DC algorithm (including two different source interval conditions and without the orthogonal constraint condition) and the SA-DC algorithm. Both algorithms utilize the CCM as the second-order statistics and have interference elimination capability. The prespecified grids for SA-DC are from [−90 • , 90 • ] with a sampling interval of 0.1 • . The regularization parameter for SA-DC was empirically chosen to be 0.4, and the penalty parameter τ in our algorithm was set to 24 when the number of snapshots is less than 1000 in Figure 4a or the SNR is less than 12 dB in Figure 4b; otherwise, τ was set to 26. Five equal-power BPSK SOIs uniformly distributed between −9 • and 9 • impinged on the array. In Figure 4a, the SNR is fixed at 0 dB, and the number of snapshots is fixed at 400 in Figure 4b. The other settings were the same as those in the first experiment. As can be seen from Figure 4a,b, even though the SA-DC algorithm can utilize the whole difference coarray information, its performance is inferior to the proposed algorithm due to the basis mismatch caused by the prespecified grids. The proposed CI-DC algorithm with source interval [−10 • , 10 • ] outperforms the condition with source interval [−90 • , 90 • ] and the case without the orthogonal constraint, especially when the SNR is low and the snapshots are limited; this indicates that a smaller source interval helps to improve the performance remarkably. When the number of snapshots is larger than 1000 in Figure 4a and the SNR is larger than 16 dB in Figure 4b, the RMSE of the proposed algorithm with different constraint conditions is similar. SNR is fixed at 0 dB, and the number of snapshots is fixed at 400 in Figure 4b. The other settings were the same as those in the first experiment. As can be seen from Figures 4a,b, even though the SA-DC algorithm can utilize the whole difference coarray information, its performance is inferior to the proposed algorithm due to the basis mismatch caused by the prespecified grids. The proposed CI-DC algorithm with source interval [ 10 ,10 ]  outperforms the condition with source interval [ 90 ,90 ]  and the case without the orthogonal constraint, especially when the SNR is low and the snapshots are limited; this indicates that a smaller source interval helps to improve the performance remarkably. When the number of snapshots is larger than 1000 in Figure 4a and the SNR is larger than 16 dB in Figure 4b, the RMSE of the proposed algorithm with different constraint conditions is similar.
(a) (b) In the third experiment, we compared the RMSE versus the number of snapshots and the SNR with the sum coarray model between our algorithm with three different constraint conditions and the SA-SC algorithm. The sources and the parameters in both algorithms were the same as those in the second experiment. It can be seen from Figure 5 that the CI-SC with a small source interval performs the best, which indicates that the orthogonal constraint incorporated with a smaller source interval helps to regulate the estimate toward a more accurate solution. This observation coincides with the result in the second experiment. In the third experiment, we compared the RMSE versus the number of snapshots and the SNR with the sum coarray model between our algorithm with three different constraint conditions and the SA-SC algorithm. The sources and the parameters in both algorithms were the same as those in the second experiment. It can be seen from Figure 5 that the CI-SC with a small source interval performs the best, which indicates that the orthogonal constraint incorporated with a smaller source interval helps to regulate the estimate toward a more accurate solution. This observation coincides with the result in the second experiment. U are two submatrices of the same dimension. The spectrum searching scheme was utilized to find all the peaks. Twenty-one equal-power BPSK SOIs uniformly distributed between 70  and 70 impinged on the same coprime array as specified in Figure 1. The SNR was set to 20 dB and the number of snapshots was 2000. The other parameters were set to be the same as those in the first experiment. After the CI-DC and CI-SC algorithms were completed, the SVD of CE R was implemented to obtain the noise subspace and the spatial spectrum was calculated based on (36). The spatial spectrum of the extended cyclic MUSIC method is illustrated in Figure 6. It can be seen that all the 21 SOIs were correctly determined using the extended cyclic MUSIC method; this result validates the hypothesis that the extended cyclic MUSIC method, utilizing the solutions from the structured matrix In the last experiment, we intended to resolve more sources than a single coarray interpolation scheme by utilizing both the obtained T (r Vd ) and H(r Vs ). We used the extended cyclic MUSIC method [10] and simply constructed the extended cyclic correlation matrix as follows: Then, the SVD of R CE was implemented to obtain the noise subspace U n , and the spatial spectrum of the extended cyclic MUSIC method was given by where U n = [U n1 ; U n2 ] is the noise subspace, and U n1 and U n2 are two submatrices of the same dimension. The spectrum searching scheme was utilized to find all the peaks. Twenty-one equal-power BPSK SOIs uniformly distributed between −70 • and 70 • impinged on the same coprime array as specified in Figure 1. The SNR was set to 20 dB and the number of snapshots was 2000. The other parameters were set to be the same as those in the first experiment. After the CI-DC and CI-SC algorithms were completed, the SVD of R CE was implemented to obtain the noise subspace and the spatial spectrum was calculated based on (36). The spatial spectrum of the extended cyclic MUSIC method is illustrated in Figure 6. It can be seen that all the 21 SOIs were correctly determined using the extended cyclic MUSIC method; this result validates the hypothesis that the extended cyclic MUSIC method, utilizing the solutions from the structured matrix completion, can achieve more DOFs than the single coarray interpolation scheme. However, rather than simply constructing the extended cyclic matrix with the respective solutions from Toeplitz completion and Hankel completion, there is potential to resolve more sources by directly solving the extended cyclic correlation matrix completion problem with structure as in (35). completion, can achieve more DOFs than the single coarray interpolation scheme. However, rather than simply constructing the extended cyclic matrix with the respective solutions from Toeplitz completion and Hankel completion, there is potential to resolve more sources by directly solving the extended cyclic correlation matrix completion problem with structure as in (35).

Conclusions
In this paper, improved coarray interpolation algorithms were proposed for cyclostationary signals. We integrate the cyclostationarity with the spatial domain and performed DOA estimation under the circumstance that the sources are more numerous than the sensors. By exploring the structure of the CCM and the CCCM in a coprime array, we built up the difference coarray and sum coarray models and formulated the coarray interpolation as a Toeplitz completion problem and a Hankel completion problem. In order to further improve the performance of the coarray interpolation, we defined the spatial spectrum sampling operators and the derivative (conjugate) correlation subspaces, then constructed orthogonal constraints for the autocorrelation vectors in the structured matrix completion problems. Prior knowledge of the source intervals was also incorporated into the problem, which helped to improve the performance significantly. Numerical results validated the effectiveness of the proposed algorithms and demonstrated their superiority in terms of interference elimination capacity, increased DOFs, and estimation accuracy.
When P is large enough, f (x) will be approximated as Consequently, we can draw the conclusion that f (θ, θ 0 ) will approximate a Kronecker delta function when P is large enough. However, for a nonuniform linear array, we can still approximate f (θ i , θ k ) with a Kronecker delta function due to the fact that f (θ i , θ k ) with θ i = θ k is much larger than the values when θ i = θ k .