DOA estimation based on sum–difference coarray with virtual array interpolation concept

The sum–difference coarray is the union of difference coarray and the sum coarray, which is capable to obtain a higher number of degrees of freedom (DOF) than the difference coarray. However, this method fails to use all information provided by the coprime array because of the existence of holes. In this paper, we introduce the virtual array interpolation into the sum–difference coarray domain. After interpolating the virtual array, we estimate the DOA by reconstructing the covariance matrix to resolve an atomic norm minimization problem in a gridless way. The proposed method is gridless and can effectively utilize the DOF of a larger virtual array. Numerical simulation results verify the effectiveness and the superior performance of the proposed algorithm.

DOF.The coprime array is able to identify up to O(MN ) sources with M + N − 1 sen- sors, where M and N are coprime.And the nested array can resolve up to O(N 2 ) sources with N physical sensors.Compared with the coprime array, the nested array has more significant mutual coupling effects.Thus, plenty of modified coprime array constructions have been proceeded [13][14][15][16][17].
Since the DCA derived from coprime array is discontiguous, the locations where there are no virtual elements in the coarray are called "holes".For the traditional subspacebased DOA estimation methods, only the longest consecutive part of the derived virtual array can be used [18].However, the information of the discontiguous virtual sensors is discarded.In order to make the coarray information fully utilized, the sparse recovery methods like l 1 norm minimization and LASSO algorithm [19] are proposed.But these methods discretize the parameter space in a dense grid which would produce the basis mismatch problem.In this case, the off-grid targets cannot accurately fall on the predefined spatial grid.The basis mismatch problem cannot be completely resolved until the gridless methods without setting grid are proposed.Currently, plenty of DOA estimation methods have been developed based on coarray interpolation by interpolating additional sensors into the positions of holes in the virtual array [20][21][22][23].In [20], a DOA estimation algorithm using coarray interpolation is proposed through nuclear norm minimization (NNM) for coprime arrays.And a gridless DOA estimation is put forward via a low-rank Toeplitz covariance matrix reconstruction approach [22].Moreover, a virtual array interpolation algorithm based on the atomic norm minimization (ANM) and Toeplitz matrix reconstruction is proposed in [23].
Recently, a novel method that generates the sum-difference coarray (SDCA) via vectorized conjugate augmented MUSIC (VCAM) algorithm has been proposed in [24].As the sum coarray (SCA) is able to fill a majority of the holes in DCA, this method can obtain larger consecutive coarray aperture, leading to much larger number of DOF.Although the virtual aperture of SDCA is more than twice the physical aperture, there still exist some holes.The nonconsecutive virtual sensors will be discarded.
In this paper, to fully utilize all information derived from SDCA, we introduce the concept of interpolation to generate a virtual ULA.The holes in the SDCA are filled up by interpolating the virtual arrays.Specifically, we interpolate the missing elements in the discontiguous virtual array into a virtual ULA at first.So that all information derived from virtual array can be fully used.Secondly, we derive an atomic norm of the output signals after virtual array interpolation.Accordingly, an ANM problem is formed to reconstruct the Toeplitz covariance matrix of the interpolated virtual array.And the ANM problem is resolved by semi-definite programming (SDP).At last, the reconstructed covariance matrix can be used to estimate DOA by employing the subspace method.Simulation results prove the superiority of the proposed method.
The rest of this paper is organized as follows.Section 2 introduces the signal model of coprime array and the generation of SDCA.Section 3 presents the proposed DOA estimation method with SDCA interpolation-based ANM algorithm.Simulation results are provided to demonstrate the effectiveness of the proposed method in Sect. 4. Finally, Sect. 5 gives the conclusion.
Notations: Lower-case (upper-case) bold characters are used to represent vectors (matrices).C and R is the set of complex number and the real number.The superscripts (•) * , (•) T , (•) H denote the complex conjugation, transpose, and conjugate transpose opera- tion, respectively.vec(•) stands for the vectorization operator.rank(•) and Tr(•) , respec- tively, imply the rank and trace of matrix.⊙ represents the Khatri-Rao product, ⊗ refers Kronecker product, and • refers the Hadamard product.|S| implies the cardinality of a set S.

Coprime array configuration
Consider a coprime array as illustrated in Fig. 1, which is composed of two sub-arrays.One sub-array has N sensors with the inter-element spacing of Md, and the other consists of M sensors with the inter-element spacing of Nd.Here, M and N are a pair of coprime integers (M < N ) .And the unit inter-element spacing d is set to be half wavelength /2 .Since the two sub-arrays share the first physical sensor, the coprime array contains L = M + N − 1 sensors.The sensor positions can be expressed as Assume the first sensor is the reference, i.e., p 1 = 0.
Suppose that Q far-field narrow-band sources from directions {θ 1 , θ 2 , . . ., θ Q } impinge on the coprime array.Then, we can represent the received signal of the tth snapshot as where A = [a(θ 1 ), a(θ 2 ), . . ., a(θ Q )] represents the array manifold matrix with a(θ q ) = 1, e j2π p 2 sin(θ q )/ , . . ., e j2π p L sin(θ q )/ T being the steering vector corresponding to the direction θ q .s(t) = [s 1 (t), s 2 (t), . . ., s Q (t)] T is the signal vector.Referring to [25], the qth deterministic source signal can be denoted as s q (t) = A q e jω q t , where A q and ω q represent the deterministic complex amplitude and the frequency offset, respectively.In addition, n(t) ∼ CN (0, σ 2 n I) is the zero-mean additive white Gaussian noise vector with covariance matrix σ 2 n I , where I is the identity matrix.

The concept of SDCA
In this subsection, we utilize the VCAM algorithm which was proposed in [24] to generate the SDCA.For the T s samples of the sensor data outputs x m (t) and x n (t) , the time average function can be calculated by (1) where τ = 0 denotes the time lag, and R s * q s q (τ ) = |A q | 2 e jω q τ .Compared with s q (t) = A q e jω q t , R s * q s q (τ ) can be seen as a signal with the power |A q | 4 coming from direc- tion θ q .Notice that the mth sensor is chosen to be the reference (i.e., m = 1 ) for the sake of analysis.Then, (3) can be simplified to where n = 1, 2, . . ., L .Putting all the R x * 1 x n into a single vector yields Then, we replace τ with −τ and take the con- jugate to get By combining v x (τ ) and [v x (−τ )] * , the conjugate augmented correlation vector v(τ ) is given by where Ā = [A T , A H ] T = [ā(θ 1 ), ā(θ 2 ), . . ., ā(θ Q )] with its qth column being ā(θ q ) = [a T (θ q ), a H (θ q )] T .Select a different set of values for τ , i.e., τ = τ s , 2τ s , . . ., N s τ s , where τ s represents pseudo-sampling period and N s represents pseudo-snapshots.Hence, the pseudo-data matrix of v(τ ) can be denoted as where T with the qth row being W q = [e jω q τ s , e jω q 2τ s , . . ., e jω q N s τ s ] .The covariance matrix of v(τ ) can be obtained with the N s pseudo-snapshots as where Besides, the qth column of Ã can be derived as According to (10), ã(θ q ) is the virtual steering vector of the virtual array.Obviously, a * (θ q ) ⊗ a(θ q ) = e j2π(−p i 1 +p i 2 ) sin(θ q )/ ( p i 1 , p i 2 ∈ P ), and a(θ q ) ⊗ a * (θ q ) = e j2π(p i 1 −p i 2 ) sin(θ q )/ ( p i 1 , p i 2 ∈ P ) are corresponding to the DCA, a * (θ q ) ⊗ a * (θ q ) = e −j2π(p i 1 +p i 2 ) sin(θ q )/ ( p i 1 , p i 2 ∈ P ) represents the negative SCA, a(θ q ) ⊗ a(θ q ) = e j2π(p i 1 +p i 2 ) sin(θ q )/ ( p i 1 , p i 2 ∈ P ) denotes the positive SCA.Conse- quently, the obtained virtual array can be treated as an SDCA, which is composed of the DCA, the positive SCA, and the negative SCA.
Definition 1 (SDCA).Suppose a coprime array specified by P , its SDCA is defined as where S SD represents the self-difference set and S CD represents the cross-difference set.S − SD and S − CD are the mirrored versions.Therefore, ignoring the unit inter-element spacing d, we can get the difference coarray positions as In the same way, S SS , S CS , S − SS and S − CS represent the self-sum set, the cross-sum set, and their mirrored versions.Hence, the sum coarray positions can be denoted as Accordingly, by averaging the elements of y in (10) corresponding to the same posi- tions in the virtual array S v , we can obtain the virtual received data of the derived virtual array S v as where Ãv is the steering matrix of the derived virtual array S v .From [24], we know that the SDCA contains some missing elements called holes.The number of DOF is limited when applying the MUSIC technique to estimate DOA.Because the discontinuous virtual sensors are discarded, the information received by the virtual array cannot be efficiently utilized.

Methods
In this section, we introduce the concept of interpolation into the SDCA to further increase the DOF.Based on the concept of array interpolation [23], additional virtual sensors are interpolated into the "holes" of the virtual array to create an interpolated ULA.Moreover, (11) an ANM problem is formulated in order to reconstruct Toeplitz covariance matrix, which can take full advantage of all information to estimate the off-grid DOAs.

Virtual array interpolation for SDCA
For the sake of intuitive understanding, an example is demonstrated in Fig. 2, where M = 3 and N = 7.
Clearly, we can obtain the nonuniform virtual array derived from the SDCA as shown in Fig. 2a, where the missing elements {±31, ±34, ±35} in S v are the holes.In order to gener- ate a virtual ULA without discarding the discontinuous sensors, the additional virtual sensors are filled into the positions of the holes, which are called the interpolated sensors.The interpolated virtual array S I is depicted in Fig. 2b, where all information received by the vir- tual sensors in S v are included.It should be noticed that the interpolated sensors expressed by hollow circles in Fig. 2b are nonfunctional and the corresponding outputs are zero.Accordingly, the received signals from the interpolated sensors are regarded as zero.Thus, the output y I of the interpolated virtual array S I can be initialized as where �•� i represents the virtual sensor at the position id.S I \ S v denotes the elements in S I but not in S v .In order to successfully use the interpolated virtual array, the unknown virtual signals received by the interpolated sensors should be recovered.
Similar to (15), the ideal received signals of the interpolated virtual array S I can be expressed as where A I = [a I (θ 1 ), a I (θ 2 ), . . ., a I (θ Q )] ∈ C |S I |×Q represents the steering matrix of the interpolated virtual array.Obviously, y I behaves like the received signal from a single snapshot.In addition, the rank of its covariance matrix is one and subspace-based DOA estimation techniques cannot be applied to identify multiple sources.To overcome this problem, we divide this interpolated virtual array S I into J = (|S I | + 1)/2 overlapping sub-arrays, and each sub-array contains J elements.Accordingly, y I can be divided into J segments {y 1 I , y 2 I , . . ., y J I } as shown in Fig. 3, and the equivalent signals of each sub-array y j I can be denoted as (16) �y . ., J ) , and a j I (θ q ) = [e −jπ d J −j+1 sin(θ q ) , e −jπ d J −j+2 sin(θ q ) , . . ., e −jπ d 2J −j sin(θ q ) ] T represents the steering vector of the qth source corresponding to the jth virtual sub-array, where d i represents the ith virtual sensor position in S I .
Specifically, the sensor positions of the reference virtual sub-array are {0, d, . . ., (J − 1)d} as shown in Fig. 3.By setting j = 1 , we can obtain the steering vector corresponding to the qth source Accordingly, the received signals of the reference virtual sub-array can be expressed as Hence, we can construct the following atomic set in continuous domain by using the steering vector of the reference virtual sub-array as The atomic norm of the output z received from the reference virtual sub-array is defined as the minimum number of atoms in A that can express z , i.e., where conv(A) represents the convex hull of the atom set A and inf is the infimum of infinite set.

Fig. 3 The virtual signals of J sub-arrays
where R Y I ∈ C J ×J is a full rank covariance matrix.Since the elements located at the holes in S I are set to be zeros, each element in R Y I suffers a deviation during the summa- tion process in (23).Therefore, R Y I can't be directly calculated.Encouragingly, according to [26], a new Hermitian Toeplitz matrix R is associated with R Y I , and they are related as R 2 = J R Y I .It is possible to construct the Hermitian Toeplitz covariance matrix R from the second-order statistics y I by Notice that the interpolated virtual array S I is symmetric about zero point.The complex- valued signals of the symmetrical pair in y I are conjugate.Thus, the matrix R is directly derived from the signal statistics y I , which includes all information received by the inter- polated virtual array.However, since there are several zero elements filled in y I which denote the holes, the diagonals in R corresponding positions are zeros.Thus, a binary vector b ∈ R J is defined to describe the presence of virtual sensors indexed in S I , where the element of 0 in b stands for the sensor positions to be interpolated and 1 otherwise.Therefore, according to (22), the Toeplitz covariance matrix reconstruction in ( 24) can be represented as the following ANM problem where T (z) denotes a Hermitian PSD Toeplitz matrix with the optimization variable z as the first column.Therefore, T (z) can be seen as the covariance matrix of signals corresponding to the reference virtual sub-array.Furthermore, notice that it contains all information of the interpolated virtual array [23].Here, B = bb T ∈ R J ×J is a binary matrix distinguishing the zero (unknown) elements and the nonzero (known) elements in R .This makes the nonzero elements in R consist of the elements in the reconstructed covariance matrix T (z) .And �•� F represents the Frobenius norm, δ is a threshold to restrict the discrepancies between the nonzero elements in R and the corresponding elements in T (z) .Moreover, T (z) 0 guarantees the Hermitian positive semi-definite Toeplitz structure of the optimum solution.Then, the optimization problem (25) can be written as where µ is a regularization parameter.(23) Since T (z) is a PSD Hermitian Toeplitz matrix, if k = rank(T (z)) < J , according to Vandermonde decomposition, T (z) can be decomposed as [27] Then, we can obtain the trace of T (z) as where Tr(•) represents the trace operator.Observing the trace of T (z) in (28) along with the definition of the atomic norm of z in (22), we can get Hence, the ANM problem in ( 26) can be equivalently transformed as, where η = 1 J µ .Note that the optimization problem (30) is convex and easy to be solved by the CVX software [28].After the optimal solution ẑ is obtained, the covariance matrix T (ẑ) of the interpolated ULA can be successfully reconstructed.At the same time, the virtual signals derived from the interpolated sensors which are initialized as zero can also be recovered.Since the signal covariance matrix T (ẑ) which corresponds to the interpolated virtual ULA is reconstructed, the subspace methods such as the MUSIC based [6,18] and the ESPRIT based [7] can be employed for DOA estimation.Here, we calculate the MUSIC spatial spectrum by the following formula where U N is the noise subspace of T (ẑ) , which is acquired by selecting the eigenvectors corresponding to the J − Q smallest eigenvalues of T (ẑ) .Finally, the DOA estimation can be acquired by searching for the locations corresponding to the Q largest peaks of the spectrum in f MUSIC (θ) .Accordingly, the number of targets that can be identified are up to J − 1.
The advantages of the proposed algorithm are summarized as follows.Firstly, the discontiguous virtual array is interpolated into a virtual ULA, and all of the information in the virtual array can be effectively utilized.Secondly, the interpolated virtual array of SDCA provides a larger coarray aperture and higher DOF.Furthermore, the atomic norm minimization problem can be performed to reconstruct the covariance matrix of the interpolated virtual array in a gridless manner, which avoids the basis mismatch problem.( 27) (30) ,

Results and discussion
In this section, several simulation experiments are illustrated to demonstrate the superiority of the proposed method in estimating the DOAs.We select the pair of coprime integers M = 3 and N = 7 to constitute the coprime array, consisting of M + N − 1 = 9 physical sensors with the positions {0, 3d, 6d, 7d, 9d, 12d, 14d, 15d, 18d} .The proposed method is compared with the dif- ference coarray interpolation-based method (DCA-I) algorithm [23] and SDCA algorithm [24].The regularization parameter η is set to be 0.25 in the simulations for the proposed method and DCA-I method.
In the first experiment, we compare the estimated spatial spectra for each algorithm.Here, we consider 15 sources uniformly distributed in the range of [−60 • , 60 • ] .The signal-to-noise ratio (SNR) is set to be 0 dB.Furthermore, the snapshots T s and the pseudo-snapshots N s are set as T s = N s = 100 .Figure 4 displays the DOA esti- mations of these methods.From Fig. 4a, we can see that the DCA-I algorithm can't detect the 15 sources correctly.Because the available DOF of the interpolated difference coarray is only 18.There exists a large DOA estimation bias under the low SNR and snapshots.In Fig. 4b, the SDCA algorithm produces several spurious peaks.It is because the discontinuous virtual sensors are discarded and the DOF is only 30.The DOA estimation ability is very poor in a low SNR and snapshots condition.In contrast, it is clear that only the proposed method can effectively detect all sources under the low SNR and snapshots as shown in Fig. 4c.Since the proposed method is capable of exploiting all information contained in the virtual array and the available DOF is further increased to 36.
In the second experiment, we compare the DOA estimation resolution by considering two closely spaced targets from θ 1 = −1 • and θ 2 = 1 • .We set SNR = 0 dB and T s = N s = 100 .In Fig. 5, we can see that the proposed method has much sharper peaks and estimates much more correctly than the others.It is because the proposed method provides a larger array aperture without discarding the discontinuous virtual sensors.It is evident that the proposed method is able to achieve higher resolution.
In the last experiment, we use the root mean square error (RMSE) to evaluate the DOA estimation performance of the three algorithms.The RMSE is defined as follows: - where K denotes the number of Monte Carlo trials and θq,k is the estimated DOA result of the qth source in the kth Monte Carlo trial.Accordingly, we consider Q = 7 and Q = 14 uncorrelated narrowband sources uniformly distributed in the range of [−60 • , 60 • ] , respectively, and we run K = 300 Monte Carlo simulations.The DCA-I algorithm can't resolve 14 sources, so its RMSE is not considered.In Fig. 6, we demonstrate the results of the RMSE versus SNR.In this simulation, we consider N s = T s = 300 and SNR var- ies from −5 to 20 dB.When Q = 14 , the proposed method obtains smaller RMSE than the SDCA method with the SNR increasing.When Q = 7 , the DCA-I method performs the worst of all, because the DCA-I method only uses the difference coarray.And the proposed method performs better than others, since the proposed method exploits all the virtual sensors of the sum coarray and difference sensors leading to higher DOF and larger array aperture than the others.In addition, Fig. 7 displays the RMSE with (32 the variation of the number of snapshots, the SNR is 0 dB and N s = T s varies from 100 to 500.Clearly, we can see that as the snapshots increase, the RMSE of the proposed algorithm decreases gradually.At the same time, the proposed method performs better than the other two algorithms.Since the proposed method makes full use all information derived from virtual sensors, which makes its performance better than others.Generally, the proposed method significantly outperforms others under a wide range of SNR or snapshots.

Conclusions
In order to effectively utilize all information derived from the virtual array, the interpolation concept is introduced into the SDCA.We propose an interpolation-based DOA estimation by reconstructing the Toeplitz covariance matrix via atomic norm minimization.The holes in the virtual array are interpolated to convert the nonuniform virtual sum-difference coarray to a ULA.Then, ANM problem is formulated to solve the Toeplitz covariance matrix reconstruction in a gridless manner.The reconstructed covariance matrix can be utilized to estimate DOA.After the SDCA interpolation, the achievable number of DOF is further increased and more sources can be identified.The simulation results validate that the proposed method provides a superior performance compared with other methods.

Fig. 2
Fig. 2 An example of the virtual coarray, where M = 3 and N = 7 .a S v , the SDCA of the coprime array.b S I , the interpolated virtual array of the SDCA

Fig. 4
Fig. 4 Spatial spectra of 15 sources, SNR = 0 dB .a The DCA-I method.b The SDCA method.c The proposed method

Fig. 5 Fig. 6
Fig. 5 Spatial spectra of two closer targets, SNR = 0 dB .a The DCA-I method.b The SDCA method.c The proposed method

Fig. 7
Fig. 7 RMSE performance versus number of snapshots