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

Authors: Hao Li; Weijia Cui; Bin Ba; Haiyun Xu; Yankui Zhang
Journal: International Journal of Antennas and Propagation
 (2020)

---

## Abstract

The 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. The 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 (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. Then, 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. The 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.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation is one of the important research fields in array signal processing [1–3]. At present, the angle parameter estimation technology based on subspace method has been relatively mature, and the typical representatives include multiple signal classification (MUSIC) 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–8]. Therefore, 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. The 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–12], which averts spectral-peak search. The 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]. This property is used to construct a polynomial, and we can obtain DOA estimation by finding the root of it. The 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].The aforementioned articles [4–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. The 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 fourth-order 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 fourth-order cumulants [26]. The 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 fourth-order cumulant matrix to improve the performance of DOA estimation. The main contributions of this paper are as follows:(1)
The 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)
The 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. Then, 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 elementsThe rest of the paper is organized as follows. In Section2, 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.Table1 shows the key notations and explanations in this paper.Table 1
Key notations and explanations.
NotationsExplanation⋅∗Conjugate of matrix⋅TTranspose of matrix⋅HConjugate transpose of the matrixE⋅Statistical expectation operatordiag⋅Diagonalization operation⋅m,nm,n-th entry of a matrixtr⋅Matrix trace⊗Kornecker productℙSet of the locations of physical elementsℤSet of the locations of virtual elementsangle⋅Phase acquired operationvec⋅Matrix column vectorization


## 2. Signal Model

The coprime array proposed in the literature [17] is considered in this paper. The coprime linear array is composed of two uniform linear subarrays, whose element numbers are MA and MB, respectively. MA and MB are mutual primes, and the spacing of subarrays is MBλ/2 and MAλ/2, respectively, where λ represents the wavelength of the carrier. The topological structure of the coprime linear array is shown in Figure 1.Figure 1
Topological structure of the coprime linear array.It is assumed that there existK target distributed sources transmitting plane waves from center directions θ=θ1,θ2,…,θK to the measuring array, which contains M=MA+MB−1 array elements. ℙ represents the aggregation of the locations of all physical array elements, where ℙ=0,d1,d2,…,dM−1 and dm represents the distance between the m+1th array element and the first array element. Then, the signal output vector can be expressed as(1)xt=∑k=1K∫−π/2π/2aθskθ,t;μkdθ+nt,where aθ=1e−j2πd/λsinθ⋅⋅⋅e−j2πM−1d/λsinθT is M×1-dimensional array steering vector, skθ,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 nt represents the additive white Gaussian noise, whose mean is zero and the variance is σn2.For the coherently distributed source, the angular signal density function can be expressed as(2)skθ,t;μk=sktgkθ;μk,where skt is the kth complex random source which reflects the time characteristics of skθ,t;μk, and gkθ;μk represents the corresponding deterministic angle signal distribution function which reflects the spatial distribution characteristics of skθ,t;μk.Furthermore, formula (1) can be rewritten as(3)xt=∑k=1Ksktbμk+nt,where bμk=∫−π/2π/2aθgkθ;μkdθ defined as the generalized direction vector of the kth coherently distributed source. The angular distribution is an exponential distribution, namely,(4)gkθ−θk;ρk=11−ρkexpjθ−θk,where ρk is directly related to the angle extension parameter σk of the distributed source 0<ρk≤1. The 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(5)bθk,ρk=ρk2d1/λe−j2πd1sinθ1/λ⋯ρk2dM/λe−j2πdMsinθk/λ.Therefore, the received signal matrix can be expressed as(6)xt=Bst+nt.


## 3. Proposed Algorithm

### 3.1. Method Introduction

Generally speaking, the FOC can be defined in two ways,cumxm1,xm2,xm3∗,xm4∗ and cumxm1,xm2∗,xm3,xm4∗, 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,(7)cumxm1,xm2∗,xm3,xm4∗=Exm1xm2∗xm3xm4∗−Exm1xm2∗Exm3xm4∗−Exm1xm3Exm2∗xm4∗−Exm1xm4∗Exm2∗xm3,where xm represents the received signal of the mth array element, ∀m1,m2,m3,m4∈1,2,…,M.Then, definingCx as an M2×M2-dimensional matrix whose element in m1−1M+m3th row and m2−1M+m4th column is cumxm1,xm2∗,xm3,xm4∗, thus, Cx can be calculated as follows:(8)Cx=Ex⊗xx⊗xH−Ex⊗xEx⊗xH−ExxH⊗ExxH−ExxH⊗ExxHmixed,where the symbol ⊗ represents the Kronecker product and the symbol ∗⊗∗mixed represents the mixed-Kronecker product which is defined as follows:(9)Ym×n⊗Zp×qmixed=Y⊗z1,Y⊗z2,…,Y⊗zq=yijz1i=1,j=1m,n,yijz2i=1,j=1m,n,⋯,yijzqi=1,j=1m,n=y11z11⋯y1nz11y11z12⋯y1nz1q⋮⋱⋮⋮⋱⋮y11zp1⋯y1nzp1y11zp2⋯y1nzpqy21z11⋯y2nz11y21z12⋯y2nz1q⋮⋱⋮⋮⋱⋮ym1zp1⋯ymnzp1ym1zp2⋯ymnzpqmp×nq.Since it is assumed thatnt 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(10)Cx=B⊗BCsB⊗BH,where Cs 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:(11)Cs=Es⊗ss⊗sH−Es⊗sEs⊗sH−EssH⊗EssH−EssH⊗EssHmixed,and independent sources have zero cross-cumulants only when i=j=k=l; we have cumsi,sj∗,sk,sl∗≠0. Hence,(12)Cs=cums1,s1∗,s1,s1∗⋱cumsK,sK∗,sK,sK∗,and we have(13)Cx=Hθ,ρCsHθ,ρT,where(14)H=bθ1,ρ1⊗bθ1,ρ1,…,bθK,ρK⊗bθK,ρK.H can be considered as the steering matrix of the virtual array, whose dimension is extended to M2×K.Considering theu−1M+vth row and the kth column element of the hθk,ρk, we have(15)hu,vθk,ρk=ρk2du+dv/λe−j2πdu+dvsinθk/λ,that is, the delay of distributed source signal from the center arrival angle θk at du+dv relative to the reference position. It is equivalent to the existence of a virtual array element at du+dv, where du,dv∈ℙ.We define the sum aggregate of physical array elements:(16)ℤ=du+dv1≤u,v≤M,which represents the location of all array elements in the extended virtual array.We take the coprime array in Section2; for example, when MA=2,MB=3, the set of positions of physical sensors is ℙ=0,2,3,4. Virtual array elements are constructed by FOC used in this paper. The positions of sensor elements and physical sensors are shown in Figure 2.Figure 2
The positions of virtual sensors and physical sensors. (a) The result of the traditional construction method of the fourth-order cumulant. (b) The result of the construction method adopted in this paper.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.The position of the largest continuous part of the virtual array element is denoted byℙ′, which contains M˜ elements. Assume that himθ,ρ is the imth element of hθ,ρ; then, the transformed virtual array steering vector is(17)h˜θ,ρ=hi1θ,ρ,…,hiM˜θ,ρT,im∈ℙ′,which is M˜-dimensional column vectors, and the corresponding array that received signal covariance can be expressed as follows:(18)Cvula=CXi1,j1⋯CXi1,jM˜⋮⋱⋮CXiM˜,j1⋯CXiM˜,jM˜,im,jm∈ℙ′,which is M˜×M˜-dimensional matrix. Thus, the subsequent computational complexity is reduced by taking the elements of Cvula out of Cx.Conducting the eigenvalue decomposition forCvula, we can obtain(19)CvulaE¯=ΛE¯,Cvula=UsΛsUs+UnΛnUn,where Λ=diagλ1,λ2,…,λM˜, λ1≥⋯≥λK>λK+1=⋯λM˜ and E¯=e1,e2,…,eM˜. Moreover, the eigenvectors corresponding to the former K large eigenvalues constitute the subspace of the signal Us, while the eigenvectors corresponding to the latter M˜−K small eigenvalues constitute the subspace of the noise Un.Therefore, 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:(20)Pθ,ρ=1h˜Hθ,ρUnUnHh˜θ,ρ.


### 3.2. Summary of the Method’s Steps

Through 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 fourth-order cumulant matrixCx 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 ofCx, and construct the matrix Cvula according to equation (18)(4)
Conduct the eigenvalue decomposition forCvula and obtain the noise subspace(5)
Apply equation (20) to solve the center DOA and distribution parameters of distributed source


### 3.3. 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. The complexity of the fourth-order cumulant mainly includes formula (14). The complexity is OTM4, where T represents the number of snapshots. The subspace decomposition algorithm is an alternative method. The 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. The former is OM˜3 and the latter is OM˜2+M˜2JθJρ. The total complexity is OM˜3+M˜2+M˜2JθJρ, where Jθ and Jρ represent the number of grids in a two-dimensional 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 O4M˜2T+8M˜3+KM˜2JθJρ. Table 2 shows the computation complexity elaborately.Table 2
Complexity comparison.
AlgorithmComplexityUniform-DSPEOTM2+M3+M2+M2JθJρUniform-LS-ESPRITO8M3+KM2JθJρ+4TM2CLA-FOC-DSPE (proposed)OTM4+M˜3+M˜2+M˜2JθJρ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 Figures3 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.Figure 3
Computation complexity with different numbers of array elements.Figure 4
Computation complexity with different numbers of the search grid.


## 3.1. Method Introduction

Generally speaking, the FOC can be defined in two ways,cumxm1,xm2,xm3∗,xm4∗ and cumxm1,xm2∗,xm3,xm4∗, 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,(7)cumxm1,xm2∗,xm3,xm4∗=Exm1xm2∗xm3xm4∗−Exm1xm2∗Exm3xm4∗−Exm1xm3Exm2∗xm4∗−Exm1xm4∗Exm2∗xm3,where xm represents the received signal of the mth array element, ∀m1,m2,m3,m4∈1,2,…,M.Then, definingCx as an M2×M2-dimensional matrix whose element in m1−1M+m3th row and m2−1M+m4th column is cumxm1,xm2∗,xm3,xm4∗, thus, Cx can be calculated as follows:(8)Cx=Ex⊗xx⊗xH−Ex⊗xEx⊗xH−ExxH⊗ExxH−ExxH⊗ExxHmixed,where the symbol ⊗ represents the Kronecker product and the symbol ∗⊗∗mixed represents the mixed-Kronecker product which is defined as follows:(9)Ym×n⊗Zp×qmixed=Y⊗z1,Y⊗z2,…,Y⊗zq=yijz1i=1,j=1m,n,yijz2i=1,j=1m,n,⋯,yijzqi=1,j=1m,n=y11z11⋯y1nz11y11z12⋯y1nz1q⋮⋱⋮⋮⋱⋮y11zp1⋯y1nzp1y11zp2⋯y1nzpqy21z11⋯y2nz11y21z12⋯y2nz1q⋮⋱⋮⋮⋱⋮ym1zp1⋯ymnzp1ym1zp2⋯ymnzpqmp×nq.Since it is assumed thatnt 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(10)Cx=B⊗BCsB⊗BH,where Cs 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:(11)Cs=Es⊗ss⊗sH−Es⊗sEs⊗sH−EssH⊗EssH−EssH⊗EssHmixed,and independent sources have zero cross-cumulants only when i=j=k=l; we have cumsi,sj∗,sk,sl∗≠0. Hence,(12)Cs=cums1,s1∗,s1,s1∗⋱cumsK,sK∗,sK,sK∗,and we have(13)Cx=Hθ,ρCsHθ,ρT,where(14)H=bθ1,ρ1⊗bθ1,ρ1,…,bθK,ρK⊗bθK,ρK.H can be considered as the steering matrix of the virtual array, whose dimension is extended to M2×K.Considering theu−1M+vth row and the kth column element of the hθk,ρk, we have(15)hu,vθk,ρk=ρk2du+dv/λe−j2πdu+dvsinθk/λ,that is, the delay of distributed source signal from the center arrival angle θk at du+dv relative to the reference position. It is equivalent to the existence of a virtual array element at du+dv, where du,dv∈ℙ.We define the sum aggregate of physical array elements:(16)ℤ=du+dv1≤u,v≤M,which represents the location of all array elements in the extended virtual array.We take the coprime array in Section2; for example, when MA=2,MB=3, the set of positions of physical sensors is ℙ=0,2,3,4. Virtual array elements are constructed by FOC used in this paper. The positions of sensor elements and physical sensors are shown in Figure 2.Figure 2
The positions of virtual sensors and physical sensors. (a) The result of the traditional construction method of the fourth-order cumulant. (b) The result of the construction method adopted in this paper.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.The position of the largest continuous part of the virtual array element is denoted byℙ′, which contains M˜ elements. Assume that himθ,ρ is the imth element of hθ,ρ; then, the transformed virtual array steering vector is(17)h˜θ,ρ=hi1θ,ρ,…,hiM˜θ,ρT,im∈ℙ′,which is M˜-dimensional column vectors, and the corresponding array that received signal covariance can be expressed as follows:(18)Cvula=CXi1,j1⋯CXi1,jM˜⋮⋱⋮CXiM˜,j1⋯CXiM˜,jM˜,im,jm∈ℙ′,which is M˜×M˜-dimensional matrix. Thus, the subsequent computational complexity is reduced by taking the elements of Cvula out of Cx.Conducting the eigenvalue decomposition forCvula, we can obtain(19)CvulaE¯=ΛE¯,Cvula=UsΛsUs+UnΛnUn,where Λ=diagλ1,λ2,…,λM˜, λ1≥⋯≥λK>λK+1=⋯λM˜ and E¯=e1,e2,…,eM˜. Moreover, the eigenvectors corresponding to the former K large eigenvalues constitute the subspace of the signal Us, while the eigenvectors corresponding to the latter M˜−K small eigenvalues constitute the subspace of the noise Un.Therefore, 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:(20)Pθ,ρ=1h˜Hθ,ρUnUnHh˜θ,ρ.


## 3.2. Summary of the Method’s Steps

Through 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 fourth-order cumulant matrixCx 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 ofCx, and construct the matrix Cvula according to equation (18)(4)
Conduct the eigenvalue decomposition forCvula and obtain the noise subspace(5)
Apply equation (20) to solve the center DOA and distribution parameters of distributed source


## 3.3. 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. The complexity of the fourth-order cumulant mainly includes formula (14). The complexity is OTM4, where T represents the number of snapshots. The subspace decomposition algorithm is an alternative method. The 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. The former is OM˜3 and the latter is OM˜2+M˜2JθJρ. The total complexity is OM˜3+M˜2+M˜2JθJρ, where Jθ and Jρ represent the number of grids in a two-dimensional 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 O4M˜2T+8M˜3+KM˜2JθJρ. Table 2 shows the computation complexity elaborately.Table 2
Complexity comparison.
AlgorithmComplexityUniform-DSPEOTM2+M3+M2+M2JθJρUniform-LS-ESPRITO8M3+KM2JθJρ+4TM2CLA-FOC-DSPE (proposed)OTM4+M˜3+M˜2+M˜2JθJρ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 Figures3 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.Figure 3
Computation complexity with different numbers of array elements.Figure 4
Computation complexity with different numbers of the search grid.


## 4. Simulation Results

In this section, the simulation experiments of the proposed method are given. The 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:(21)RMSE=1NK∑j=1N∑i=1Kθi−θ∧i,j2,where θi is the true DOA of the distributed source, θ∧i,j is the estimated DOA coordinate of the ith distributed source by the jth Monte Carlo experiment, and N is the number of Monte Carlo. The 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(22)CRLBμ=1TUHΠV⊥U−1,where ΠV⊥=I−VVHV−1VH. Refer to Appendix B for the detailed derivation of (22). The CRLB derived from the maximum-likelihood criterion will be plotted in the results of simulation experiments to measure the performance of the proposed method.Simulation 1.
Estimate ability of different distribution parameters.
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. The subarraysMA=2MB=3. The DOA of the three distributed sources is evenly distributed from −20° to 20°, and the parameter ρ=0.7, ρ=0.8, and ρ=0.9, respectively. The 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. The 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.Figure 5
Scatter diagram of parameter estimation (SNR = 0 dB).Figure 6
Scatter diagram of parameter estimation (SNR = 10 dB).Figure 7
Scatter diagram of parameter estimation (SNR = 20 dB).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. The proposed algorithm improves the degree of freedom by constructing a virtual array, which can estimate more sources than the number of array elements. The subarraysMA=2MB=3. The DOA of the seven distributed sources is evenly distributed from −50° to 50°, and the parameter ρ=0.8. The 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. The 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.Figure 8
Normalized spectral peak search (SNR = 10 dB).Figure 9
Normalized spectral peak search (SNR = 20 dB).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. The DOA of the three distributed sources is evenly distributed from−20° to 20°, and the parameter ρ=0.8. The number of snapshots is set as T = 2000. The subarrays MA=2MB=3 in the proposed algorithm. The 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 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.Figure 10
RMSE of DOA estimates versus SNRs.Simulation 4.
RMSE of DOA estimates versus the number of snapshots.
In this part, we experiment on the effect of different snapshots on the performance of the algorithm. The SNR is 20 dB, and the other parameters are identical to those in Simulation3. 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.Figure 11
RMSE of DOA estimates versus the number of snapshots.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. The number of snapshots is fixed to 2000 and the SNR is set as 10 dB. The other parameters are the same as those in Simulation 3. Based on 500 Monte Carlo experiments, the simulation results are shown 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 12
RMSE of DOA estimates versus different parametersρ.Simulation 6.
RMSE of DOA estimates versus different step lengthsΔθ.
In order to research the performance of the proposed algorithm under different step lengths of searching, we conduct 500 Monte Carlo experiments under the condition of different step lengthsΔθ with the number of snapshots being fixed to 2000. The SNR ranges from −4 dB to 20 dB and Δρ=0.01. The other parameters are identical to those in Simulation 3. The simulation results are shown in Figure 13. It is obvious that the smaller the step length Δθ, the better the performance, but at the same time, the greater the complexity. Therefore, in practical applications, it is necessary to balance the weight of complexity and accuracy to obtain the best estimation effect.Figure 13
RMSE of DOA estimates versus different step lengths of searching.


## 5. Conclusion

In order to achieve high precision and high DOF estimation of multiple distributed sources, this paper introduces fourth-order 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. Then, 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.



---

*Source: 8968136-2020-07-11.xml*