Unfolded Coprime Linear Array with Three Subarrays for Non-Gaussian Signals: Configuration Design and DOA Estimation

In this paper, we investigate the problem of sparse array design for the direction of the arrival (DOA) of non-Gaussian signals and exploit the unfolded coprime linear array with three subarrays (UCLATS) to obtain physical sensors location. With the motivation from the large consecutive degree of freedom (DOF), we optimize the process of obtaining physical sensors location from two steps. Specifically, the first is to model the process of obtaining the longest consecutive virtual sum co-array from a given number of physical array elements into a global postage-stamp problem (GPSP), whose solution can be employed to determine the locations of the longest possible consecutive sum co-array (2-SC) and initial physical array. The second step is to multiply the location of the virtual sum co-array by appropriate coprime coefficients to generate UCLATS and then multiply the initial physical array position by the same corresponding coefficients to obtain physical sensors location. Besides, an algorithm is proposed to obtain DOA estimates, which employs the discrete Fourier transform (DFT) method and partial spectrum searching multiple signal classification (PSS-MUSIC) algorithm to obtain initial estimates and fine estimates, respectively, termed as the DFT-MUSIC method. Compared with the traditional total spectrum searching MUSIC (TSS-MUSIC) algorithm, the DFT-MUSIC method performs the same asymptotical performance of DOA estimation with less than 10% complex multiplication times, which can be verified by numerical simulations under the same condition.


Introduction
Array signal processing (ASP) exploits the sensor array to receive spatial signals in order to obtain discrete observation data. Compared with a single directional sensor, the sensor array has many advantages such as stronger spatial gain, more flexible beam control and higher spatial angle resolution; thus, it has been widely used in beamforming, direction of arrival (DOA) estimation, source separation and microphone arrays for acoustic imaging over the decades [1][2][3][4][5]. Especially, DOA estimation plays an important role in array signal processing, which is widely applied in wireless communication system, radar system and navigation [1,6]. Conventional DOA estimation methods mainly focus on uniform linear array (ULA) [7,8] or uniform planar array (UPA) for their simple and symmetric structure, whose inter-element spacing is required to be no longer than half wavelength to avoid angle ambiguity, which results in limited DOA estimation accuracy.
We design the structure of UCLATS, whose cDOF of 2-DC is also provided. The problem of sparse array design with non-Gaussian signals is investigated from GPSP perspective and the structure of UCLATS.

2.
We divide the process of obtaining the location of physical sensors into two steps, which optimizes the array location step by step, resulting in lower design difficulty.

3.
We devise DFT-MUSIC algorithm for a better balance between computational complexity and DOA estimation performance, which can be utilized for the non-Gaussian signals with the proposed array geometry.
In Section 2, we introduce 2-DC, 2-SC, 2-DCSC, FODC and the properties of UCLATS. Section 3 describes the GPSP and the procedure of proposed sparse array geometry design principle. Section 4 presents the DFT-MUSIC method. Section 5 shows the performance analysis. Section 6 elaborates the simulation results, and Section 7 draws the conclusions.
Notations: Matrices and vectors are denoted by bold lowercase and uppercase characters, respectively. ⊗ and represent the Kronecker product and Khatri-Rao product, respectively. (·) T , (·) H and (·) * stand for the transpose, conjugate transpose and complex conjugation operation, respectively. vec(·) means the vectorization operation of a matrix. Cum(·) represents the cumulant operator and round(·) is the rounding operation.

Preliminaries
In this part, 2-DC, 2-SC, 2-DCSC, FODC, the properties of UCLATS and the received signal model are introduced.

2-DC, 2-SC, 2-DCSC, FODC
Definition 1. The 2-DC location D can be defined as [11] where D + and D − represent the positive and negative elements of D, respectively. D sel f and D cross are self-difference co-array set and cross-difference co-array set, respectively. S denotes the physical sensors location set.
where S sel f and S cross denote self-sum co-array set and cross-sum co-array set, respectively.
where S and D represent the physical sensors and 2-DC location set, respectively.

Definition 4.
The FODC location D 4dc can be defined as [26,28] where S and D represent the physical sensors and 2-DC location set, respectively. Therefore, it can be concluded that the 2-DCSC and FODC of physical arrays can be transformed to each other from (3) and (4).

The Properties of UCLATS
The configuration of UCLATS has been shown in Figure 1, which contains three sparse uniform subarrays overlapped at the origin. The position of UCLATS can be given by where M and N are coprime integers, which represent the physical sensors of subarray 1 S 1 and subarray 2 S 2 , respectively. d = λ/2 and λ denote wavelength. In addition, the total physical sensors number of UCLATS is T = M + N + N/2 − 1. (a) The DOF of UCLATS has been presented in Table 1

Sparse Array Design Principle
Motivated from the advantages of larger cDOF generated by the 2-DCSC of physical array, according to Equation (4) and [30], the process of obtaining 2-DCSC of physical sensors can be optimized step by step. Specifically, the first is to model the process of obtaining the consecutive longest 2-SC from given number as a GPSP, whose solution can be utilized to determine the locations of initial physical arrays and virtual 2-SC. The second is to multiply virtual 2-SC by appropriate coprime coefficients to form UCLATS, consequently, the physical array location can be obtained by multiplying the corresponding coprime coefficients to initial physical arrays location.

Global Postage-Stamp Problem (GPSP)
According to [31,32], we describe the GPSP as: for given positive integers h and k, a set with k non-negative integers can be given by Remark 1. The elements in S k should be summed h times to achieve the consecutive integers range with 0, 1, 2, · · · , n h (S k ).
Remark 2. n h (S k ) needs to be as large as possible.
The solutions to the GPSP have been given in [30][31][32]; we also present the results shown in Table 2.

Sparse Array Design Principle Based on UCLATS
Suppose that the physical sensors are composed of three sparse linear subarrays with M, N and N/2 + 1(N < M) elements, respectively. According to Section 3.1, the sets S N , S M and S N/2 can be obtained from the solution to GPSP. For S N , if all elements are negative, the consecutive integers set that can be obtained from N sensors are where S N represents a set with N integers, n 2 (S N ) is the largest possible integer mentioned in (12), which can be obtained by Table 2. For S M and S N/2 , if all elements are positive, the consecutive integers set which can be obtained from M and N/2 + 1 sensors are The 2-SC of S N , S M and S N/2 can be associated with three co-subarrays of UCLATS if we multiply Equations (13) and (14) by appropriate coprime coefficients. The co-subarray 1 of UCLATS associated with (13) can be denoted as where d 1 = (n 2 (S M ) + 1)d represents the element spacing of co-subarray 1 and co-subarray 3, d = λ/2 and λ is wavelength. The co-subarray 2 and co-subarray 3 location sets associated with (14) can be expressed as where d 2 = (n 2 (S N ) + 1)d represents the element spacing of co-subarray 2. Besides, n 2 (S M ) + 1 and n 2 (S N ) + 1 are coprime integers.
The co-array S sc = S sM ∪ S sN ∪ S s N/2 can be considered as UCLATS. Considering that the 2-SC of S N , S M and S N/2 have been multiplied by d 1 or d 2 , then the location of physical sensors should also be adjusted by multiplication operations. We give the location set of physical sensors as follows, Figure 2 shows the example of the process mentioned above, where M = 5,

DOA Estimation Method
Conventional MUSIC methods cost expensive computational complexity due to the TSS process. In this section, we propose a DFT-MUSIC method for the proposed array geometry. Specifically, the DFT [33] method for the consecutive 2-DCSC part of the proposed array geometry is utilized to obtain the initial DOA estimates. Subsequently, the PSS-MUSIC algorithm is exploited to obtain the fine DOA estimates.

Initial Estimation by DFT Algorithm
According to [33], the DFT is a non-parametric spectrum analysis method whose DOA estimation performance depends on the number of consecutive virtual elements. Based on Section 2.2, we have illustrated that the 2-DC of UCLATS contains consecutive co-arrays; besides, the relationship between the 2-DCSC of S and the 2-DC of S sc has been also provided. Thus, the consecutive virtual elements part of proposed array geometry with range [−Wd, Wd] can be employed by DFT method to obtain initial estimation, whose total number is T = 2W + 1, the value of W has been given in (28).
Define the steering vector a c (θ k ) ∈ C T×1 as Define the normalized DFT matrix F ∈ C T×T , whose element corresponding to coordinate (p, q) can be expressed as Consequently, the new normalized steering vectorã c (θ k ) constructed by F and a c (θ k ) can be written asã whose q-th element is The analysis of the value in Equation (21) has been given in Appendix B, which can be concluded from two aspects. Specifically, when T sin θ k /2 = q k is an integer, only the value of q = q k -th element ofã c (θ k ) is not zero, when q k is not an integer, most elements in (21) are small values close to zero, except the q = round{q k }-th element.
Therefore, the initial estimates of θ k can be obtained by finding the approximate position of the non-zero elements inã c (θ k ), it can also be used to obtain the number of sources. In applications, the source estimates can be obtained by performing DFT on the vectorized received signal z 1 , which is obtained by sorting and removing the redundancy operation of z. Further, the vector after DFT is y ini = Fz 1 , whose spectrum corresponds to K largest peaks are marked as q ini k (k = 1, 2, · · · , K); then, the initial DOA estimates can be calculated by

Fine Estimation by PSS-MUSIC Method
According to Appendix B, when q k = T sin θ k /2 is not an integer, more accurate DOA estimates cannot be obtained by q ini k ; thus, PSS-MUSIC is introduced to tackle this problem. Generally, MUSIC methods achieve improved DOA estimation performance, but TSS results in expensive computational complexity. The DFT-MUSIC utilizes PSS under known initial estimates, which reduces complexity effectively. This part first utilizes the DFT method described in Section 4.1 to obtain the initial estimates and then employs the PSS-MUSIC method to obtain fine estimation. The process of spatial smoothing with T consecutive co-arrays is depicted in Figure 3, where T = 2W + 1, and the value of W is shown in (28). According to [35], spatial smoothing methods require a full rank FOC matrix. By employing spatial smoothing process, the consecutive co-array steering vector a c (θ k ) shown in Equation (18) is divided into W + 1 sub-arrays, and each sub-array contains W + 1 elements, where the position of i-th(1 ≤ i ≤ W + 1) subarray can be denoted by whose vectorized FOC matrix z 1i is constructed within row W + 1 − i to 2W + 1 − i of z 1 . The covariance matrix can be constructed by Sum the covariance matrices of all W + 1 subarrays, whose mean value corresponding to full rank spatial smoothing matrix can be calculated by Based on the eigenvalue decomposition (EVD), R can be decomposed as where E s is formed by the eigenvectors corresponding to the maximum K eigenvalues, and E n is composed of the rest eigenvectors. D s and D n are diagonal matrices, the diagonal elements of D s are made up of the largest K eigenvalues, and the diagonal elements of D n are composed of other eigenvalues. Based on the orthogonal relationship between the noise space E n and the steering vector, the spectrum function of the fine estimation in the range [θ ini k − ∆, θ ini k + ∆], ∆ = 1 • can be represented by where a sub (θ) is the steering vector of subarray i(1 ≤ i ≤ W + 1) described in Figure 3, which can be set as a sub (θ) = [1, e −j2πd sin θ/λ , · · · , e −jπWd sin θ/λ ] . Thus, the fine estimatesθ k , (k = 1, · · · , K) can be obtained by performing PSS via (27).

Performance Analysis
In this section, we compare the performance of the proposed array geometry and DFT-MUSIC method with other arrays and algorithms from the viewpoints of DOF and computational complexity.

Achievable Consecutive DOF
It is noteworthy to analyze the cDOF of the proposed array. Considering that the S sM , S sN and S s N/2 have been associated with the co-subarrays of UCLATS in Figure 2, according to Lemma 1, the cDOF of the 2-DC of S sc can be obtained by  However, for the proposed array, the following relationship is established: where cDOF proposed denotes the cDOF of proposed array geometry, which can be calculated by the 2-DCSC or FODC of physical sensors.
As mentioned in Definition 1 and Definition 2, the 2-DC are composed of self-difference co-array and cross-difference co-array, and the 2-SC is composed of self-sum co-array and cross-sum co-array. Considering that the processes of obtaining three subarrays of proposed geometry are separate, the cross-sum co-array of proposed geometry is ignored. Consequently, the 2-SC of S s contains S sc completely. We also give the non-negative part of 2-DCSC of the proposed array geometry, as shown in the bottom part of Figure 4. It is obvious that Equation (29) is established. Besides, the consecutive range of cDOF UCLATS is also contained by cDOF proposed .
The cDOF comparisons of ACA [14], SAFE-CPA [28], TL-NA [11], FL-NA [26] and proposed array geometry are depicted in Figure 5. Considering that the SOC-based array cannot exploit the properties of 2-DCSC adequately, the proposed array geometry, SAFE-CPA and FL-NA achieve better cDOF performance than TL-NA and ACA. Besides, the proposed array geometry performs better than SAFE-CPA and FL-NA, as the total number of sensors increases. In the case of the same number of sensors, different combinations of subarrays may produce different cDOF; we also present the largest possible combinations of cDOF for each array structure, as shown in Tables 3-5.     Remark 3. The cDOF of the proposed array geometry listed in Tables 4 and 5 is the minimum cDOF, which can be calculated from cDOF UCLATS , as shown in (28). Considering the existence of the sum co-array between the sub-arrays, the actual cDOF will be a bit larger. e.g., (T = 10, actual cDOF = 269 > 251).

Computational Complexity
For non-Gaussian signals, we propose a DFT-MUSIC method for the proposed geometry described in Section 4, which reduces computational complexity of MUSIC algorithm by shrinking the range of spectrum searching. The complexity comparison of TSS-MUSIC algorithm and DFT-MUSIC are listed to verify the efficiency, which can be calculated by the times of complex multiplications. For DFT-MUSIC method, calculating the FOC matrix requires O J(M + N + N/2 ) 4 , obtaining initial estimates needs O T 2 , and PSS process costs O{nW(2(W − K) + 1)}, where n and T represent the PSS times and the number of cDOF, respectively. Besides, W = (T + 1)/2. For the TSS-MUSIC method, the computational complexity can be calculated by O J(M + N + N/2 − 1) 4 + W 3 + n 1 W(2(W − K) + 1) , where n 1 represents the TSS times.
The subarrays sensors of the proposed array are set to be M = 5, N = 4, snapshots J = 1000, sources number K = 2. The step of searching grid is ∆ = 0.01 • . Thus, the searching times of TSS and PSS are n = (2 • /∆)K = 400, n = (180 • /∆) = 18, 000. According to the discussion mentioned above, the spectrum searching range of DFT-MUSIC and TSS-MUSIC methods are (θ ini k − 1 • , θ ini k + 1 • ), k = 1, · · · K and (−90 • , 90 • ), respectively. The value of the complex multiplications times has been given in Table 6. Because the DFT-MUSIC method performs PSS rather than TSS, it narrows the searching range effectively. Hence, the computational complexity can be reduced effectively, which means that the DFT-MUSIC is computationally efficient.

Simulations Results
In this section, we exploit 500 Monte Carlo simulations via the DFT-MUSIC algorithm to validate the performance of the proposed array geometry. The root mean square error (RMSE) is given by where θ k represents the elevation of k-th target,θ k,i is estimates of θ k in the i-th (i = 1, . . . , 500) simulation. Besides, the searching interval of DFT-MUSIC algorithm is ∆ = 0.01 • . Suppose that there are K uncorrelated far-field non-Gaussian sources with θ = [5 • , 40 • ] incident on different arrays listed in Table 4

RMSE Performance Comparison versus Snapshots
To validate the effectiveness of the proposed array and DFT-MUSIC method in parameter estimation, Figure 6 shows the RMSE performance comparison versus snapshots with proposed array geometry by utilizing the DFT-MUSIC method, where the total number of physical sensors is T = 10. As depicted in Figure 6, the increase in the number of snapshots means a more accurate FOC matrix of received signal model, which results in better DOA estimation performance.

RMSE Performance Comparison versus Sensors
The RMSE performance comparison of the proposed array geometry versus sensors via DFT-MUSIC method is given in Figure 7. As the number of sensors increases, the diversity gains of the received signal increase. Consequently, a better DOA estimation performance can be obtained. Besides, when the total number of sensors is T = 8, the location of the proposed array is {−13, −6.5, 0, 2.5, 6.5, 7.5, 12.5, 15}, and T = 9 corresponds to {−17, −8.5, 0, 2.5, 7.5, 8.5, 12.5, 17.5, 20}.

RMSE Performance Comparison of Different Arrays with DFT-MUSIC Method
This part depicts the RMSE performance comparison of the proposed array, ACA [14], SAFE-CPA [28], TL-NA [11], FL-NA [26] with DFT-MUSIC algorithm, where the total number of sensors is T = 10. Besides, the snapshots of Figure 8 is J = 1000,the SNR of Figure 9 is SNR = −10dB.  As shown in Figures 8 and 9, the DOA estimation performance of all arrays has better parameter estimation performance with the increase in snapshots and SNR. In addition, considering that ACA and TL-NA are all designed based on SOC, which can be applied to vectorized FOC methods by adding 2-DC operations. Consequently, the virtual array elements locations obtained from the first 2-DC operation are not fully utilized, resulting in much redundancy. The SAFE-CPA, FL-NA and the proposed array designed based on FOC have better parameter estimation performance than the array designed based on FOC because of larger cDOF, and the proposed array achieves the maximum cDOF, which results in the best DOA estimation performance.

RMSE Performance Comparison of Different Algorithms
To verify the superiority of the DFT-MUSIC method, Figure 10 shows the RMSE performance comparison of DFT-MUSIC, TSS-MUSIC, initial estimation of DFT, SS-PM and SS-ESPRIT algorithm, where the interval of searching is ∆ = 0.01 • and the number of snapshots is J = 1000. Besides, the searching ranges of DFT-MUSIC and TSS-MUSIC are (θ ini k − 1 • , θ ini k + 1 • ), k = 1, · · · K and (−90 • , 90 • ), respectively, and θ ini k represents the initial DOA estimates.
It can be concluded from Figure 10 that the DFT-MUSIC algorithm owns the same asymptotical DOA estimation performance as the TSS-MUSIC algorithm, whereas Section 5.2. points out that the computational complexity of DFT-MUSIC is much lower than the TSS-MUSIC algorithm, and the DFT-MUSIC method achieves a better balance between complexity and DOA estimation accuracy; hence, it has better application prospects.

Conclusions
In this paper, a sparse array design method for non-Gaussian signals was proposed, which was based on GPSP and the properties of UCLATS. Further, a sparse array design example was given to illustrate the procedure. The process of obtaining physical sensors location could be divided into two steps; specifically, the first step was to model the problem of obtaining the longest possible consecutive virtual 2-SC from a given number of physical array elements into a GPSP and then to utilize the solution to GPSP to obtain the position of the virtual 2-SC and the position of the initial physical array. The second step was to multiply the positions of virtual 2-SC by appropriate coefficients to generate UCLATS. Consequently, the position of physical sensors could be obtained by multiplying the corresponding coefficient to initial physical array. Besides, the DFT-MUSIC algorithm was proposed, which achieved the same parameter estimation asymptotic performance with less than 10% computational complexity as compared with the TSS-MUSIC algorithm under the same simulation conditions.

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

Appendix A
(a) Considering the coprime relationship between M and N, the structure of UCLATS can be rewritten as Without loss of generality, suppose that M > N. According to Definition 1, the difference co-array D UCLATS includes the self-difference co-array set D sUCLATS and crossdifference co-array set D cUCLATS , specifically, Case 1: N is odd, For the cross-difference co-array set D cUCLATS , which can be simplified as It can be also verified in Figure 4.

Appendix B
For the convenience of presentation, Equation (21) can be rewritten as, Considering the e −j T−1 T π(q− T 2 sin θ k ) = 1, the discussions mainly focus on the value of sin[π(q− T 2 sin θ k )] sin[ π T (q− T 2 sin θ k )] . We analyze the value of ã sort (θ k ) q from two aspects, which are elaborated in the Case I and Case II below. Besides, the value range of parameters mentioned above can be given as, Case I. If T 2 sin θ k = q k ∈ Z, considering q ∈ Z,then sin π q − T 2 sin θ k sin π T q − T 2 sin θ k = sin(kπ) where k = q − T 2 sin θ k ∈ − T 2 + 1, 3T 2 , and obviously, Actually, the value of k T can be determined by (A17), where q = q k = T 2 sin θ k . When q = q k , q k , q ∈ Z, sin( k T π) = 0, sin(kπ) sin( k T π) Case II. If T 2 sin θ k = q k / ∈ Z, considering q ∈ Z, then sin π q − T 2 sin θ k sin π T q − T 2 sin θ k = sin(kπ) Obviously, k, k T / ∈ Z,i.e., sin(kπ) sin( k T π) = 0. Based on the discussion of Case I, if q = round{q k }, the value of |q − q k | ≈ 0(T 1). Then, the value of (A23) is close to √ T. When q = round{q k }, q ∈ Z, q k / ∈ Z, considering the relationship sin(kπ) sin( k T π) ≤ 1, e −j T−1 2 ( 2π T q−π sin θ k ) = 1, T 1 (A24) The ã sort (θ k ) q will be a small value.