Central DOA Estimation Method for Exponential-Type Coherent Distributed Source Based on Fourth-Order Cumulant

(e performance of direction-of-arrival (DOA) estimation for sparse arrays applied to the distributed source is worse than that applied to the point source model. In this paper, we introduce the coprime array with a large array aperture into the DOA estimation algorithm of the exponential-type coherent distributed source. In particular, we focus on the fourth-order cumulant (FOC) of the received signal which can provide more useful information when the signal is non-Gaussian than when it is Gaussian. (e proposed algorithm extends the array aperture by combining the sparsity of array space domain with the fourthorder cumulant characteristics of signals, which improves the estimation accuracy and degree of freedom (DOF). Firstly, the signal-received model of the sparse array is established, and the fourth-order cumulant matrix of the received signal of the sparse array is calculated based on the characteristics of distributed sources, which extend the array aperture. (en, the virtual array is constructed by the sum aggregate of physical array elements, and the position set of its maximum continuous part array element is obtained. Finally, the center DOA estimation of the distributed source is realized by the subspace method. (e accuracy and DOF of the proposed algorithm are higher than those of the distributed signal parameter estimator (DSPE) algorithm and least-squares estimation signal parameters via rotational invariance techniques (LS-ESPRIT) algorithm when the array elements are the same. Complexity analysis and numerical simulations are provided to demonstrate the superiority of the proposed method.


Introduction
Direction-of-arrival (DOA) estimation is one of the important research fields in array signal processing [1][2][3]. At present, the angle parameter estimation technology based on subspace method has been relatively mature, and the typical representatives include multiple signal classification (MU-SIC) method [4] and estimation signal parameters via rotational invariance techniques (ESPRIT) method [5], but these algorithms are all based on the point source model. However, in actual communication, a large number of multipath phenomena are caused due to the reflection and scattering of a complex environment, which makes the signal source expanded at a certain angle in space. Under this condition, using the traditional algorithm for point source model, the estimation performance will be seriously deteriorated [6][7][8]. erefore, it is essential to study the model and algorithm of the distributed source.
Because of the spatial distribution characteristics of the distributed source, the difference between the mathematical models of point source and the distributed source is mainly reflected in the difference of the steering matrix which causes the point source algorithm to be no longer applicable to the distributed source model. e idea of the MUSIC method based on the point source model is applied to a distributed source model by using the principle of subspace orthogonality in literature [9], namely, the DSPE algorithm and LS-ESPRIT algorithm. However, to obtain the central DOA and angle spread, a two-dimensional spectral-peak search is required and the computation complexity is high. In order to reduce the complexity, the algorithm is proposed, where an approximate rotation invariant relation is constructed to solve the central DOA using a special array layout [10][11][12], which averts spectral-peak search. e paper [13] estimates the rotation invariant matrix by constructing a propagation operator to solve the central DOA without eigenvalue decomposition and with lower computational complexity. It is proved that, for any centrosymmetric array, the deterministic angular signal distribution function (DADF) vector has symmetric rotation invariance in literature [14]. is property is used to construct a polynomial, and we can obtain DOA estimation by finding the root of it. e paper [15] realizes the central DOA estimation of two-dimensional incoherent distributed sources with low complexity by using three parallel arrays and noncircular characteristics of signals. When it comes to large-scale MIMO, a beam space algorithm is applied to settle the DOA estimation of incoherently distributed sources [16]. e aforementioned articles [4][5][6][7][8][9][10][11][12][13][14][15][16] are all based on the research of uniform arrays. Sparse arrays, typically represented by coprime arrays [17] and nested arrays [18], have the advantage of large array aperture than that of the uniform array under the circumstance of the same number of array elements and are worthy of further application in DOA estimation algorithm of the distributed source model. Based on the nested matrix model, virtual array elements are constructed by vectorizing the covariance matrix in documents [19]. An efficient approach for estimating the DOA of coherent signals using coprime arrays is proposed in the literature [20]. However, vectorization converts the covariance matrix into a single snapshot receiving a model of a virtual array. e parameters are solved by spatial smoothing, which causes the aperture of the virtual array to be sacrificed. For Gaussian signals, the first-order and second-order statistics can completely describe their statistical characteristics. However, at present, most of the artificial communication signals, like binary phase-shift keying (BPSK) and M-ary amplitude shift keying (MASK), are non-Gaussian [21,22], and their higher-order cumulants contain information, from which, we can gain more useful information to further improve the estimation performance. In the point source model, a MUSIC-like method is proposed based on the fourth-order cumulant matrix (FCM) [23]. FOC is utilized to expand the array aperture, which effectively improves the estimated number of sources compared with the method based on the second-order statistics. A redundancy-elimination fourthorder cumulant DOA estimation algorithm is proposed, which removes the redundant information of the FOC matrix in the MUSIC-like algorithm and significantly reduces the computational complexity [24,25]. An algorithm is proposed for MA system identification using the third-and fourthorder cumulants [26]. e above fourth-order cumulant methods are all proposed for point source models, whose performance will deteriorate when applied to distributed source based on sparse arrays.
To solve the above problems, this paper combines the characteristics of exponential distribution with the fourthorder cumulant matrix to improve the performance of DOA estimation.
e main contributions of this paper are as follows: (1) e specific angular distributed form of the distributed source is not considered in the aforementioned algorithms; therefore, in this paper, the exponential-type coherent distributed source is considered, and the corresponding DOA estimation algorithm with high performance is designed by analyzing the characteristics of steering vector and combining with the large aperture characteristics of mutual matrix. At last, a novel parameter estimation algorithm of exponential-type coherent distributed source based on the fourth-order cumulant under sparse array is proposed (2) e proposed algorithm extends the array aperture by combining the sparsity of array space domain with the fourth-order cumulant characteristics of signals, which improves the estimation accuracy and degree of freedom. Firstly, the fourth-order cumulant matrix of the received signal is calculated based on the characteristics of distributed sources, which extends the array aperture. en, the virtual array is constructed by the sum aggregate of physical array elements, and the position set of its maximum continuous part array element is obtained. Finally, the center DOA estimation of the exponential-type coherent distributed source is realized by the subspace method. Compared with the contrast algorithm, the proposed algorithm improves the estimation accuracy and the number of estimable sources under the condition of the same number of array elements e rest of the paper is organized as follows. In Section 2, the exponential-type coherent distributed source model under sparse array is given. According to the signal model, a DOA estimation algorithm of distributed source based on the fourth-order cumulant is proposed in Section 3. In Section 4, performance simulation experiments are carried out and simulation results are given. Our conclusions and directions for future research are presented in Section 5. Table 1 shows the key notations and explanations in this paper.

Signal Model
e coprime array proposed in the literature [17] is considered in this paper. e coprime linear array is composed of two uniform linear subarrays, whose element numbers are M A and M B , respectively. M A and M B are mutual primes, and the spacing of subarrays is (M B λ/2) and (M A λ/2), respectively, where λ represents the wavelength of the carrier. e topological structure of the coprime linear array is shown in Figure 1.
It is assumed that there exist K target distributed sources transmitting plane waves from center directions θ � (θ 1 , θ 2 , . . . , θ K ) to the measuring array, which contains M � M A + M B − 1 array elements. P represents the aggregation of the locations of all physical array elements, where P � 0, d 1 , d 2 , . . . , d M− 1 and d m represents the distance between the m + 1 th array element and the first array element. en, the signal output vector can be expressed as where a(θ) � 1e − j2π(d/λ)sin θ · · · e − j2π(M− 1)(d/λ)sin θ T is M × 1-dimensional array steering vector, s k (θ, t; μ k ) is the angular signal density function of the kth source, μ k � (θ k , σ k ) is the angle parameter of the kth coherently distributed source, θ k and σ k , respectively, represent the center DOA and angle spread, and n(t) represents the additive white Gaussian noise, whose mean is zero and the variance is σ 2 n . For the coherently distributed source, the angular signal density function can be expressed as where s k (t) is the kth complex random source which reflects the time characteristics of s k (θ, t; μ k ), and g k (θ; μ k ) represents the corresponding deterministic angle signal distribution function which reflects the spatial distribution characteristics of s k (θ, t; μ k ). Furthermore, formula (1) can be rewritten as where b(μ k ) � π/2 − (π/2) a(θ)g k (θ; μ k )dθ defined as the generalized direction vector of the kth coherently distributed source. e angular distribution is an exponential distribution, namely, where ρ k is directly related to the angle extension parameter σ k of the distributed source (0 < ρ k ≤ 1). e closer ρ k is to 1, the smaller the angle spread is. In this case, the distributed source generalized direction vector of the exponential distribution can be expressed as erefore, the received signal matrix can be expressed as

Method Introduction.
Generally speaking, the FOC can be defined in two ways, cum( , respectively. For distributed sources, because of the particularity of the steering vector, the definition of the conventional fourth-order cumulant, namely, the former, is no longer applicable to the distributed source. Refer to Appendix A for a detailed analysis. Instead, the FOC of the received signal can be defined as the later, namely, where x m represents the received signal of the mth array element, , thus, C x can be calculated as follows: where the symbol ⊗ represents the Kronecker product and the symbol [ * ⊗ * ] mixed represents the mixed-Kronecker product which is defined as follows: � y 11 z 11 · · · y 1n z 11 y 11 z 12 · · · y 1n z 1q ⋮ ⋱ ⋮ ⋮ ⋱ ⋮ y 11 z p1 · · · y 1n z p1 y 11 z p2 · · · y 1n z pq y 21 z 11 · · · y 2n z 11 y 21 z 12 · · · y 2n z 1q Since it is assumed that n(t) is Gaussian noise and independent of the signal, therefore, its FOC is identical to zero, which can be substituted into the above equation; then, we can obtain where C s is the FOC matrix of the signal s and is K × K dimensions. When the signal is independent stationary non-Gaussian, the following formula is given: and independent sources have zero cross-cumulants only and we have where H can be considered as the steering matrix of the virtual array, whose dimension is extended to M 2 × K.
Considering the ((u − 1)M + v) th row and the kth column element of the h(θ k , ρ k ), we have that is, the delay of distributed source signal from the center arrival angle θ k at (d u + d v ) relative to the reference position.
It is equivalent to the existence of a virtual array element at We define the sum aggregate of physical array elements: which represents the location of all array elements in the extended virtual array. We take the coprime array in Section 2; for example, when M A � 2, M B � 3, the set of positions of physical sensors is P � 0, 2, 3, 4 { }. Virtual array elements are constructed by FOC used in this paper. e positions of sensor elements and physical sensors are shown in Figure 2.
We can observe that both fourth-order cumulant construction methods can extend the array aperture, and the expansion ability is not the same. Because the traditional construction method of the fourth-order cumulant can construct more virtual arrays, it is widely used in the point source model, but not in the distributed source model. Although the construction method of the fourth-order cumulant used in this paper is not as capable as that of the traditional method in extending the array aperture, it also obtains a larger DOF than the physical array. It is vital that the method used in this paper is applicable to the distributed source model. e position of the largest continuous part of the virtual array element is denoted by P′, which contains M elements. Assume that h i m (θ, ρ) is the i m th element of h(θ, ρ); then, the transformed virtual array steering vector is which is M-dimensional column vectors, and the corresponding array that received signal covariance can be expressed as follows: which is M × M-dimensional matrix. us, the subsequent computational complexity is reduced by taking the elements of C vula out of C x . Conducting the eigenvalue decomposition for C vula , we can obtain where Λ � diag(λ 1 , λ 2 , . . . , λ M }, λ 1 ≥ · · · ≥ λ K > λ K+1 � · · · λ M and E � e 1 , e 2 , . . . , e M . Moreover, the eigenvectors corresponding to the former K large eigenvalues constitute the subspace of the signal U s , while the eigenvectors corresponding to the latter M − K small eigenvalues constitute the subspace of the noise U n . erefore, according to the principle of the DSPE algorithm, the following formula can be constructed to search the two-dimensional spectral peak, and the center DOA and distribution parameter of the distributed source can be solved:

Summary of the Method's Steps.
rough the introduction of the above algorithm, the steps of this method are summarized as follows: (1) Use the array receiving the signal to estimate the fourthorder cumulant matrix C x according to equation (8) (2) Construct the virtual array according to formula (16), and obtain the position set of its maximum continuous part array element (3) Get rid of the redundancy of C x , and construct the matrix C vula according to equation (18) (4) Conduct the eigenvalue decomposition for C vula and obtain the noise subspace (5) Apply equation (20) to solve the center DOA and distribution parameters of distributed source

Complexity Analysis.
According to the introduction and summary of the above algorithm, the following two aspects make up the complexity of this algorithm: the solution of the fourth-order cumulant and the application of the subspace decomposition algorithm. e complexity of the fourth-order cumulant mainly includes formula (14). e complexity is O(TM 4 ), where T represents the number of snapshots. e subspace decomposition algorithm is an alternative method. e complexity of solving the DOA estimation by using the DSPE method includes the complexity of eigenvalue decomposition and that of a two-dimensional spectral-peak search.
e former is O(M 3 ) and the latter is where J θ and J ρ represent the number of grids in a twodimensional spectral-peak search. J θ � (180°/Δθ) + 1 and J ρ � (1/Δρ) + 1, where Δθ and Δρ represent the step length of researching. LS-ESPRTI algorithm has advantages over the DSPE algorithm in reducing complexity, whose complexity is Table 2 shows the computation complexity elaborately.
To get more intuition about the computational complexity, we experiment on computation complexity with different numbers of array elements and different numbers of the search grid. It can be seen from Figures 3 and 4 that the proposed method increases the computational complexity to some extent but improves the estimation accuracy and degree of freedom, so the complexity is within the accepted range.

Simulation Results
In this section, the simulation experiments of the proposed method are given. e root mean square error (RMSE) is regarded as the performance index of the proposed algorithm, and the definition of RMSE, which is the deviation between the estimated DOA and the true DOA, is given as follows: where θ i is the true DOA of the distributed source, θ ∧ i,j is the estimated DOA coordinate of the i th distributed source by the j th Monte Carlo experiment, and N is the number of Monte Carlo. e Cramer-Rao Lower Bound (CRLB) describes the best performance of the parameter estimation method in theory [27,28]. We derive the CRLB about parameter vector μ as where Π ⊥ V � I − V(V H V) − 1 V H . Refer to Appendix B for the detailed derivation of (22). e CRLB derived from the   In this part, we draw the parameter coordinates estimated by each Monte Carlo experiment into a scatter diagram and describe the effectiveness of the estimation by its aggregation degree. e subarrays M A � 2M B � 3. e DOA of the three distributed sources is evenly distributed from − 20°to 20°, and the parameter ρ � 0.7, ρ � 0.8, and ρ � 0.9, respectively. e step length of searching is Δθ � Δρ � 0.01. We conduct simulation experiments under the following three conditions: (1) SNR � 0 dB; (2) SNR � 10 dB; (3) SNR � 20 dB. For 200 independent trials, the simulation results are shown in Figures 5-7, corresponding to the three cases mentioned above, respectively. e black dots are the estimation results of all Monte Carlo experiments and the cross in red represents the true value. We observe that the simulation results in discrete distributions around true value; moreover, with the increase of SNR, the clustering degree of scattering becomes stronger; therefore, the proposed algorithm can give the correct estimations when different distributed sources have different parameters.

Simulation 2. Performance simulation under underdetermined condition.
Underdetermined condition refers to the case where the number of sources is larger than the number of array elements. e proposed algorithm improves the degree of freedom by constructing a virtual array, which can estimate more sources than the number of array elements. e subarrays M A � 2M B � 3. e DOA of the seven distributed sources is evenly distributed from − 50°to 50°, and the parameter ρ � 0.8. e number of snapshots is T �10000 and the step length of searching is Δθ � Δρ � 0.01. We conduct simulation experiments under the following two conditions: (1) SNR � 10 dB; (2) SNR � 20 dB. e normalized spectral peak is shown in Figure 8 when the SNR is 10 dB and in Figure 9 when it is 20 dB. In the figure, the red line is the spectral peak search curve, and the blue dashed line represents the true incoming wave direction. We observe that the position of the spectral peak is very close to the real value, so the proposed algorithm has the ability to estimate under undetermined conditions.

Simulation 3. RMSE of DOA estimates versus SNRs.
In order to research the performance of the proposed algorithm under different SNRs, we compare it with the DSPE algorithm and the LS-ESPRIT algorithm. e DOA of the three distributed sources is evenly distributed from − 20°t o 20°, and the parameter ρ � 0.8. e number of snapshots is set as T � 2000. e subarrays M A � 2M B � 3 in the proposed algorithm. e step length of searching is Δθ � Δρ � 0.01. To show the superiority of the proposed algorithm, the number of arrays is set as 4 in the DSPE algorithm and the distance between adjacent sensors is 0.5λ. Figure 10 shows the RMSE of the central DOA estimations with SNR ranging from − 4 dB to 20 dB based on the 500 Monte Carlo experiments. We observe that, under the same SNR, the estimation accuracy of the proposed method is     In this part, we experiment on the effect of different snapshots on the performance of the algorithm. e SNR is 20 dB, and the other parameters are identical to those in Simulation 3. Figure 11 shows the RMSE of the central DOA estimations with the number of snapshots ranging from 400 to 8000 based on 500 Monte Carlo experiments. We observe that, under the same number of snapshots, the estimation accuracy of the proposed method is higher in contrast with that of the DSPE algorithm and the LS-ESPRIT algorithm when the array elements are the same and closer to the CRLB.

Simulation 5. RMSE of DOA estimates versus different parameters ρ.
In this part, we experiment on the effect of different parameters ρ on the performance of the algorithm. e number of snapshots is fixed to 2000 and the SNR is set as 10 dB. e other parameters are the same as those in  Figure 12. We observe that, with the increase of the parameter ρ, the performance of RMSE decreases. When the parameter ρ � 1, the distributed source model is reduced to the point source model; thus, one conclusion can be drawn; that is, the point source model is a particular case of a distributed source model.  Figure 13. It is obvious that the smaller the step length Δθ, the better the performance, but at the same time, the greater the complexity. erefore, in practical applications, it is necessary to balance the weight of complexity and accuracy to obtain the best estimation effect.

Conclusion
In order to achieve high precision and high DOF estimation of multiple distributed sources, this paper introduces fourthorder cumulant into distributed sources and presents a central DOA estimation method for the exponential-type coherent distributed source based on the fourth-order cumulant. Firstly, an exponential-type coherent distributed source signal model under sparse array is established, and according to the signal model, we combine the sparsity of array space domain with the fourth-order cumulant characteristics of signals. In this way, we extend the array aperture and degree of freedom and improve the estimation accuracy. en, in order to get rid of redundancy, we get the maximum continuous part of the extended virtual array element. Finally, we can obtain the DOA through the subspace decomposition algorithm. Compared with the contrast algorithm based on the second-order statistics, the proposed algorithm improves the estimation accuracy and the number of estimable sources under the condition of the same number of array elements. Complexity analysis and numerical simulations are provided to demonstrate the superiority of the proposed algorithm.

A. Detailed Reason Why the Former FOC Is No Longer Applicable to Distributed Source
e conventional FOC of the received signal can be expressed as where x m represents the received signal of the m th array element. en, defining C x as an M 2 × M 2 -dimensional matrix whose element in ((m 1 − 1)M + m 3 ) th row and ((m 2 − 1)M + m 4 ) th column is cum(x m 1 , x m 2 , x * m 3 , x * m 4 ), thus, C x can be calculated as follows: where the symbol ⊗ represents the Kronecker product. Since n(t) is Gaussian noise and independent of the signal, therefore, its FOC is identical to zero, which can be substituted into the above equation; then, we can obtain where C s is the FOC matrix of the signal s and is K × K dimensions. When the signal is independent stationary non-Gaussian, the following formula is given: and independent sources have zero cross-cumulants only when i � j � k � l; we have cum(s i , s j , s * k , s * l ) ≠ 0. Hence, and we have C x is the covariance matrix of the output signal of the virtual array after expansion, where H can be considered as the steering matrix of the virtual array, whose dimension is extended to M 2 × K.
Considering the ((u − 1)M + v) th row and the k th column element of the h(θ k , ρ k ), we have From the formula above, we can see that the extended steering vector does not have the form of a distributed source steering vector. So, the traditional construction of the fourth-order cumulant is no longer applicable to the distributed source model. where

(B.2)
And then, the vector μ is denoted by where θ � (θ 1 , θ 2 , . . . , θ K ), ρ � (ρ 1 , ρ 2 , . . . , ρ K ). e vector υ is denoted by where σ � (σ 2 s 1 , σ 2 s 2 , . . . , σ 2 s K ) represents the power of signals, and σ 2 n represents the power of noise. en, we define a vector η as which contains all the unknown parameters. e CRLB of η can be calculated from [29]: Finally, we can prove that Data Availability e data, which are produced by simulations, used to support the findings of this study are available from the corresponding author upon request. Disclosure e authors claim that the data used in this article are provided by their simulations and this is developed without using any data in a published article to support their results.

Conflicts of Interest
e authors declare that they have no conflicts of interest.