# Mixed Far-Field and Near-Field Source Localization Algorithm via Sparse Subarrays

Authors: Jiaqi Song; Haihong Tao; Jian Xie; Chenwei Sun
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

Based on a dual-size shift invariance sparse linear array, this paper presents a novel algorithm for the localization of mixed far-field and near-field sources. First, by constructing a cumulant matrix with only direction-of-arrival (DOA) information, the proposed algorithm decouples the DOA estimation from the range estimation. The cumulant-domain quarter-wavelength invariance yields unambiguous estimates of DOAs, which are then used as coarse references to disambiguate the phase ambiguities in fine estimates induced from the larger spatial invariance. Then, based on the estimated DOAs, another cumulant matrix is derived and decoupled to generate unambiguous and cyclically ambiguous estimates of range parameter. According to the coarse range estimation, the types of sources can be identified and the unambiguous fine range estimates of NF sources are obtained after disambiguation. Compared with some existing algorithms, the proposed algorithm enjoys extended array aperture and higher estimation accuracy. Simulation results are given to validate the performance of the proposed algorithm.


---

## Body

## 1. Introduction

In recent years, passive source localization has become a key topic in array signal processing [1]. Various localization algorithms have been proposed for far-field (FF) source, whose wavefronts are plane waves, such as the multiple signal classification (MUSIC) method [2], the estimation of signal parameters via rotational invariance technique (ESPRIT) [3], and their derivatives. Nevertheless, when the radiating sources are located in near-field (NF) source, whose wavefronts are spherical waves, both the DOA and the range parameters should be determined to localize these radiating sources. As a result, traditional FF DOA estimation algorithms are no longer applicable for NF source localization. Fortunately, many advanced methods have been presented under the NF assumption, including the 2-D MUSIC algorithm [4], the high-order ESPRIT algorithm [5, 6], the covariance approximation (CA) method [7, 8], the weighted linear prediction method [9], the generalized ESPRIT algorithm [10, 11], and the signal reconstruction for near-field source localization [12]. In addition, we also proposed a novel method of passive localization for near-field noncircular sources [13].However, both FF and NF sources may coexist in many interested situations such as speaker localization using microphone arrays, seismic exploration, and electronic surveillance. Most of the algorithms, which deal with pure NF or pure FF sources, may fail in the scenarios of mixed sources.Recently, mixed source localization problem has been an important research topic in array signal processing [14–19]. A two-stage MUSIC (TSMUSIC) algorithm [14] was first advanced to localize mixed FF and NF sources. Based on fourth-order cumulant, the TSMUSIC algorithm can successfully estimate the parameters of mixed sources. However, its computational complexity is high due to the construction of high-order cumulant matrices and the spectral search. To relieve the computational burden, an efficient MUSIC-based algorithm is proposed in [15], which only utilizes second-order statistics. Unfortunately, this algorithm suffers from severe array aperture loss and 1-D spectral search in both DOA and range estimation. Based on [10], Liu and Sun presented another GESPRIT-based algorithm to alleviate the array aperture and obtain a reasonable classification result [16]. Our early works on mixed localization focused on unknown source numbers [17] and unknown mutual coupling [18].As is well known, the estimation accuracy is directly correlated with the array aperture size that a larger array would produce more precise estimates. However, most of the existing methods limit the array element spacing to be within a quarter wavelength to avoid DOA ambiguity. Recently, a mixed-order MUSIC (MOMUSIC) algorithm [19] was proposed to extend the array aperture by a special nested sparse linear array (SLA), which can improve the estimation accuracy.In this paper, a novel fourth-order cumulant-based dual-size shift invariance (CDSSI) algorithm is presented to solve the mixed source localization problem. Our technique utilizes a SLA of the dual-size spatial invariance method [20] which was designed by M. D. Zoltowski and K. T. Wong. Furthermore, they also proposed MUSIC/MODE null spectrum disambiguation algorithm [21] aiming to realize more accurate disambiguation. Unlike most of the existing algorithms, our technique utilizes a SLA of dual-size spatial invariance. As a result, the novel method enjoys significant promotion in DOA and range estimation accuracy by extending the intersubarray spacing. Furthermore, the proposed method avoids any 1-D or 2-D spectral searching and therefore has lower computational complexity.The rest of the paper is organized as follows. In Section2, the mixed FF and NF signal model based on SLA is presented. The proposed CDSSI algorithm is described in Section 3. In Section 4, we compare CDSSI with some recently developed ones, and computer simulations are conducted in Section 5 to validate the performance of the proposed algorithms. Finally, we conclude this paper in Section 6.Notation.
The complex conjugate, transpose, Hermitian transpose, and pseudoinverse are denoted by⋅∗, ⋅T, ⋅H, and ⋅#, respectively. ⊗ symbolizes the Kronecker product and ∘ represents the Khatri-Rao product (or column-wise Kronecker product), that is, a1,…,aJ∘b1,…,bJ=a1⊗b1,…,aJ⊗bJ. IJ is a J×J identity matrix.


## 2. Data Model

Suppose thatK-independent narrowband sources (FF and NF) impinge upon a symmetric SLA, as shown in Figure 1. The array has a total number of 2M+1 sensors, which are composed of 3 subarrays and each subarray contains 2Mu+1 array sensors, with M and Mu being positive integers. This SLA configuration can be considered as a particularization of the sparse rectangular dual-size spatial invariance array in [20].Figure 1
Sparse linear array geometry.In Figure1, the intersensor spacing is d, and the intersubarray spacing is ds=d⋅Δ, where ds≫d. The sensor position vector P would have the following form, if we take d as the length unit.
(1)P=p−M,p−M+1,…,p0,…,pM−1,pMT=−Δ−Mu,…,−Δ,…,−Δ+Mu⏟subarray−1,−Mu,…,0,…,Mu⏟subarray0,Δ−Mu,…,Δ,…,Δ+Mu⏟subarray1T.Let the array center be the phase reference point; the output of themth sensor can be approximated as [5]
(2)xmt=∑k=1Ksktejpmωk+pm2ϕk+nmt,t=1,2,…,L,(3)ωk=−2πdλsinθk,ϕk=πd2λrkcos2θk,where θk represents the DOA of the kth source, rk stands for the distance between the kth source and the reference sensor, skt symbolizes the kth source signal, nmt denotes the additive Gaussian noise, and L is the number of snapshots. It should be noted that, for the NF source scenario, the range parameter lies in the Fresnel region 0.62D3/λ,2D2/λ [22], with D symbolizing the array aperture. For the FF source scenario, the range parameter approaches to ∞ and the associated parameter ϕk becomes 0 [14].In a matrix form, the array data can be written as(4)Xt=ASt+Nt,where
(5)Xt=x−Mt,…,x0t,…,xMtT,(6)A=aθ1,r1,…,aθk,rk,…,aθK,rK,(7)aθk,rk=aωk,φk=ejωkP+ϕkP2,(8)St=s1t,…,skt,…,sKtT,Nt=n−Mt,…,n0t,…,nMtT.In the above equations,A represents the array steering matrix, aωk,ϕk denotes the array steering vector for the kth source, St symbolizes the source signal vector, and Nt is the noise vector.Given the array dataXt, a novel algorithm is proposed in Section 3 to obtain high performance in localizing and distinguishing the mixed sources successfully, under the following hypotheses.
(i)
The incoming signals are mutually independent, narrowband stationary, and non-Gaussian, as well as nonzero kurtosis.
(ii)
The DOAs of all source signals differ from each other.
(iii)
The noise is zero-mean, additive (white or color) Gaussian, and statistically independent from all impinging sources.
(iv)
In order to avoid the phase ambiguity, the intersensor spacingd should be within a quarter wavelength.


## 3. Algorithm Development

As is shown in Figure2, the proposed algorithm has two main stages. In the first stage, by constructing a fourth-order cumulant matrix C1 containing only the DOA information, we can decouple the DOA estimation from the range estimation. After the eigendecomposition of C1, dual-size shift invariance ESPRIT [20] could be applied to generate coarse (unambiguous) and fine (ambiguous) DOA estimates of both NF and FF sources. In the second stage, another fourth-order cumulant matrix C2 needs to be constructed to overcome the rank-deficient phenomenon described in [14]. Moreover, the virtual steering matrix of C2 can be decoupled into two parts and both the coarse and fine range estimates of these sources could be obtained. In both of the two stages, disambiguation is important for the final accurate and unambiguous localization of mixed FF and NF sources. The proposed CDSSI algorithm is described in detail in the following subsections.Figure 2
Flow graph of the CDSSI algorithm.


### 3.1. DOA Estimation

#### 3.1.1. Cumulant Matrix for the DOA Estimation

Similar to the TSMUSIC algorithm in [14], we first construct a fourth-order cumulant matrix to estimate the DOAs of the radiating sources in this subsection. However, there are some differences between the two algorithms. For example, (1) the uniform linear array (ULA) configuration in TSMUSIC is generalized to a SLA in Figure 1, which promotes the estimation accuracy effectively and (2) the proposed method utilizes ESPRIT to generate the DOA estimation, which avoids the spectral search procedure in TSMUSIC and therefore reduces the computational complexity.According to the definition in [23], the fourth-order cumulant of the array outputs xmt, xnt, xit, and xjt can be written as
(9)cumxmt,xn∗t,xi∗t,xjt=∑k=1Kcskejpm−pn−pi−pjωk+pm2−pn2−pi2−pj2ϕk,where csk=cumskt,sk∗t,sk∗t,skt symbolizes the kurtosis of the kth source.We can construct a2M+1×2M+1 cumulant matrix C1, with its m,nth element being
(10)C1m,n=cumxm−M−1t,x−m+M+1∗t,xn−M−1∗t,x−n+M+1t=∑k=1Kcskej2pm−M−1ωk⋅e−j2pn−M−1ωk,m,n∈1,2M+1.Note thatC1 can be written in a matrix form
(11)C1=BC4sBH,where
(12)C4s=diagcs1,…,csk,…,csK,(13)B=bω1,…,bωk,…,bωK,bωk=ej2ωkP.From (1), B can be written as
(14)B=AsubD−1AsubAsubD1=Ds∘Asub,where
(15)D=e−j2Δω1⋱e−j2Δωk⋱e−j2ΔωK,Ds=d1,…,dk,…,dK,dk=e−j2Δωk,1,ej2ΔωkT,herein, Asub is the virtual steering matrix for each subarray in the cumulant domain, with the following form:
(16)Asub=a¯1,…,a¯k,…,a¯K,a¯K=e−j2Muωk,e−j2Mu−1ωk,…,1,…,ej2Mu−1ωk,ej2MuωkT.Therefore,C1 can be expressed as
(17)C1=Ds∘AsubC4sDs∘AsubH.The reason to rewrite it in this manner is that it has a similar form of dual-size invariance and therefore, ESPRIT could be applied to yield coarse and fine estimates of DOAs.


#### 3.1.2. Fine DOA Estimates with Ambiguities

We firstly estimate the DOAs of the mixed FF and NF sources. By taking an eigendecomposition of (17), we have
(18)C1=USΛSUSH+UNΛNUNH,where ΛS and ΛN are diagonal matrices containing K large eigenvalues and 2M+1−K small eigenvalues, respectively. UN is the 2M+1×2M+1−K eigenvector matrix spanning the noise subspace of C1. US is the 2M+1×K signal subspace eigenvector matrix of C1, which can be expressed as
(19)US=Ds∘AsubT,where T is a nonsingular matrix.To obtain high-accuracy DOA estimates, we may form two matrices that are related by the extended intersubarray spacingds along the x-axis. As is illustrated in Figure 3, by taking the first and the last 22Mu+1 rows of US, we have the following two matrices.
(20)US1=US1:4Mu+2,:=D1∘AsubT,US2=US2Mu+1:end,:=D2∘AsubT=D1∘AsubΦfT,where D1 and D2 are the first and the last 2 rows of Ds and Φf=diagej2Δω1,…,ej2ΔωK. Since US1 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray−1 and subarray0 and US2 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray0 and subarray1, fine DOA estimates can be found from Φf due to the extended invariance relationship.Figure 3
Construction of 2 subspace matrix elements for fine DOA estimation.From (20), we have US2=US1T−1ΦfT, which means that there is a rotational invariance between US1 and US2, that is, Ψf=T−1ΦfT=US1#US2. Therefore, the diagonal elements of Φf, Φfk,k=ej2Δωk=ej4πdssinθk/λ, correspond to the eigenvalues of Ψf. As ds≫λ/4 and sinθk≤1, we can find a series of ambiguous DOA estimates
(21)θ^klθ=arcsinμk+lθλ2ds,k=1,…,K,(22)2dsλ−1−μk≤lθ≤2dsλ1−μk,(23)μk=∠Φfk,k4πds/λ,where lθ is an integer and ⋅ and ⋅ represent the ceil (round toward positive infinity) and floor (round toward negative infinity) operation in MATLAB, respectively. ∠⋅ returns the phase angle of the operand, which lies between −π,π.


#### 3.1.3. Coarse DOA Estimates without Ambiguities

For the purpose of disambiguation, unambiguous but coarse DOA estimates must be obtained as references to the ambiguous fine estimates in the previous subsection. From Figure4, by extracting the cumulant domain information of the first and the last 2Mu sensors in each subarray, the quarter-wavelength spatial invariance of the array geometry can be exploited to generate unambiguous coarse DOA estimates. This can be easily done by defining a permutation matrix
(24)Γ=I3⊗e1,I3⊗e2,…,I3⊗e2Mu+1,where ei signifies the ith column of the identity matrix I2Mu+1.Figure 4
Construction of 2 subspace matrix elements for coarse DOA estimation.From (19) and (24), we have
(25)UP=ΓUS=Asub∘DsT.By taking the first and the last3×2Mu rows of UP, we obtain the following two matrices.
(26)UP1=UP1:6Mu,:=A1∘DsT,UP2=UP4:end,:=A2∘DsT=A1∘DsΦcT,where A1 and A2 are the first and the last 2Mu rows of Asub, respectively. Φc=diagej2ω1,…,ej2ωK. Since UP1 is contributed by the cumulant domain data from the first 2Mu sensors in each subarray and UP2 is contributed by the cumulant domain data from the last 2Mu sensors in each subarray, unambiguous DOA estimates with quarter-wavelength spatial invariance can be found from Φc.From (26), we have UP2=UP1T−1ΦcT, which means that there is a rotational invariance between UP1 and UP2, that is, Ψc=T−1ΦcT=UP1#UP2. Therefore, the diagonal elements of Φc, Φck,k=ej2ωk=ej4πdsinθk/λ, correspond to the eigenvalues of Ψc. The unambiguous coarse DOA estimates of K-mixed FF and NF sources are given by
(27)θ^kref=arcsin∠Φck,k4πd/λ,k=1,…,K.θ^kref and θ^klθ may be paired by the method in [20].


#### 3.1.4. Disambiguation

The coarse estimatesθ^kref could be used as references to disambiguate the cyclic ambiguities in θ^klθ. Mathematically, we obtain the disambiguated angle estimates as
(28)θ^k=arcsinμk+lθ0λ2ds,where the estimate of lθ0 is given by
(29)lθ0=argminlθθ^klθ−θ^kref.Note thatlθ0 is solved through searching over a few discrete points and the searching range is given by (22).


### 3.2. Range Estimation

#### 3.2.1. Cumulant Matrix for the Range Estimation

In order to derive the range estimates, another fourth-order cumulant matrix needs to be constructed to overcome the rank-deficient phenomenon described in [14]. The virtual steering matrix for range estimation in [5] has the following form
(30)A~=11⋯1ej2φ1ej2φ2⋯ej2φK⋮⋮⋯⋮ej2N−1φ1ej2N−1φ2⋯ej2N−1φK.If thekth source is in FF, the electrical angle ϕk will become 0 and the kth column of A~ will be 1,…,1T. When the number of FF sources is more than one, A~ will drop rank and thus the algorithm will fail to localize radiating sources.In this section, a new fourth-order cumulant matrixC2 is defined, with its m,nth element being
(31)C2m,n=cumx0t,x0∗t,xm−M−1t,xn−M−1∗t=∑k=1Kcskejpm−M−1ωk+pm−M−12ϕk⋅ejpn−M−1ωk+pn−M−12ϕk∗,m,n∈1,2M+1.Note thatC2 can be expressed as
(32)C2=cumx0t,x0∗t,Xt,X∗t=AC4sA∗,where A is the steering matrix in (6) and C4s is the kurtosis diagonal matrix of sources as defined in (12). When the kth source is in FF, the corresponding column vector aωk,ϕk becomes
(33)aωk,ϕk=ejωkp−M,…,ejωkp0,…,ejωkpMT.Therefore,A will still be of full column rank when there are multiple FF sources.


#### 3.2.2. Unambiguous and Ambiguous Range Estimations

From (1) and (7), aωk,ϕk can be decoupled as
(34)aωk,ϕk=Vωkgϕk,Vωk=ejωk−Δ−Mu⋱ejωk−Δ+Muejωk−Mu⋱1⋰ejωkMuejωkΔ−Mu⋰ejωkΔ+Mu,(35)gϕk=ejϕkΔ+Mu2,…,ejϕkΔ2,…,ejϕkΔ−Mu2,ejϕkMu2,…,ejϕk12,ejϕk02T,where Vωk is a 2M+1×M+1 matrix which only contains the DOA information and gϕk is a M+1×1 column vector which only depends on ϕk.By taking an eigendecomposition ofC2, we can obtain a 2M+1×K signal subspace matrix ES and a 2M+1×2M+1−K noise subspace matrix EN, where ES is composed of the K principle eigenvectors of C2 and EN contains the remaining 2M+1−K eigenvectors.According to [4], the 2-D MUSIC spectrum for DOA and range estimation is given by
(36)Pω,ϕ=aHω,ϕENENHaω,ϕ.Substituting the disambiguated DOA estimatesθ^kk=1,2,…,K into the above equation, the 2-D spectrum search can be reduced to the K 1-D ones. Therefore, the estimates of ϕkk=1,2,…K are given by
(37)ϕ^k=minϕaHω^,ϕENENHaω^,ϕ=minϕgHϕVHω^ENENHVω^gϕ.Actually, (37) implies that gϕk is the eigenvector associated with the minimum eigenvalue of the Hermitian matrix VHω^kENENHVω^k. Through the eigendecomposition of VHω^kENENHVω^k, we can obtain an estimation of gϕk.In order to generate the coarse range estimates unambiguously,ϕk can be estimated as
(38)ϕ^k=∠gϕkMgϕkM+1,where gϕkn represents the nth element in gϕk. From (3), the coarse range estimate is given by
(39)r^kref=πd2λϕ^kcos2θ^k,where θ^k is the disambiguated DOA estimation. According to the definition of the Fresnel region, NF sources are located in the range of 0.62D3/λ,2D2/λ. As a result, when r^kref is in this region, the corresponding source is classified as a NF one; otherwise, it is regarded as a FF one.Similarly, through the extended aperture size between each subarray, fine range estimates with ambiguity can be generated. From (35), we have
(40)gϕkMu+1gϕkM+1=ejΔ2ϕk=ejπds2/λrkcos2θk.Therefore, we can find a series of ambiguous range estimates(41)r^klr=πds2νk+2lrπλcos2θ^k,(42)−νk2π≤lr≤ds22λr^kref−νk2π,(43)νk=∠gϕkMu+1gϕkM+1.


#### 3.2.3. Disambiguation

In an analogous manner to that of Section3.1.4, the coarse estimates r^kref could be used as references to disambiguate the cyclic ambiguities in r^klr. Mathematically, the disambiguated range estimates are
(44)r^k=πds2νk+2lr0πλcos2θ^k,where the estimate of lr0 is given by
(45)lr0=argminlθr^klr−r^kref.Note thatlr0 is solved through searching over a few discrete points, and the searching range is given by (42).


## 3.1. DOA Estimation

### 3.1.1. Cumulant Matrix for the DOA Estimation

Similar to the TSMUSIC algorithm in [14], we first construct a fourth-order cumulant matrix to estimate the DOAs of the radiating sources in this subsection. However, there are some differences between the two algorithms. For example, (1) the uniform linear array (ULA) configuration in TSMUSIC is generalized to a SLA in Figure 1, which promotes the estimation accuracy effectively and (2) the proposed method utilizes ESPRIT to generate the DOA estimation, which avoids the spectral search procedure in TSMUSIC and therefore reduces the computational complexity.According to the definition in [23], the fourth-order cumulant of the array outputs xmt, xnt, xit, and xjt can be written as
(9)cumxmt,xn∗t,xi∗t,xjt=∑k=1Kcskejpm−pn−pi−pjωk+pm2−pn2−pi2−pj2ϕk,where csk=cumskt,sk∗t,sk∗t,skt symbolizes the kurtosis of the kth source.We can construct a2M+1×2M+1 cumulant matrix C1, with its m,nth element being
(10)C1m,n=cumxm−M−1t,x−m+M+1∗t,xn−M−1∗t,x−n+M+1t=∑k=1Kcskej2pm−M−1ωk⋅e−j2pn−M−1ωk,m,n∈1,2M+1.Note thatC1 can be written in a matrix form
(11)C1=BC4sBH,where
(12)C4s=diagcs1,…,csk,…,csK,(13)B=bω1,…,bωk,…,bωK,bωk=ej2ωkP.From (1), B can be written as
(14)B=AsubD−1AsubAsubD1=Ds∘Asub,where
(15)D=e−j2Δω1⋱e−j2Δωk⋱e−j2ΔωK,Ds=d1,…,dk,…,dK,dk=e−j2Δωk,1,ej2ΔωkT,herein, Asub is the virtual steering matrix for each subarray in the cumulant domain, with the following form:
(16)Asub=a¯1,…,a¯k,…,a¯K,a¯K=e−j2Muωk,e−j2Mu−1ωk,…,1,…,ej2Mu−1ωk,ej2MuωkT.Therefore,C1 can be expressed as
(17)C1=Ds∘AsubC4sDs∘AsubH.The reason to rewrite it in this manner is that it has a similar form of dual-size invariance and therefore, ESPRIT could be applied to yield coarse and fine estimates of DOAs.


### 3.1.2. Fine DOA Estimates with Ambiguities

We firstly estimate the DOAs of the mixed FF and NF sources. By taking an eigendecomposition of (17), we have
(18)C1=USΛSUSH+UNΛNUNH,where ΛS and ΛN are diagonal matrices containing K large eigenvalues and 2M+1−K small eigenvalues, respectively. UN is the 2M+1×2M+1−K eigenvector matrix spanning the noise subspace of C1. US is the 2M+1×K signal subspace eigenvector matrix of C1, which can be expressed as
(19)US=Ds∘AsubT,where T is a nonsingular matrix.To obtain high-accuracy DOA estimates, we may form two matrices that are related by the extended intersubarray spacingds along the x-axis. As is illustrated in Figure 3, by taking the first and the last 22Mu+1 rows of US, we have the following two matrices.
(20)US1=US1:4Mu+2,:=D1∘AsubT,US2=US2Mu+1:end,:=D2∘AsubT=D1∘AsubΦfT,where D1 and D2 are the first and the last 2 rows of Ds and Φf=diagej2Δω1,…,ej2ΔωK. Since US1 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray−1 and subarray0 and US2 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray0 and subarray1, fine DOA estimates can be found from Φf due to the extended invariance relationship.Figure 3
Construction of 2 subspace matrix elements for fine DOA estimation.From (20), we have US2=US1T−1ΦfT, which means that there is a rotational invariance between US1 and US2, that is, Ψf=T−1ΦfT=US1#US2. Therefore, the diagonal elements of Φf, Φfk,k=ej2Δωk=ej4πdssinθk/λ, correspond to the eigenvalues of Ψf. As ds≫λ/4 and sinθk≤1, we can find a series of ambiguous DOA estimates
(21)θ^klθ=arcsinμk+lθλ2ds,k=1,…,K,(22)2dsλ−1−μk≤lθ≤2dsλ1−μk,(23)μk=∠Φfk,k4πds/λ,where lθ is an integer and ⋅ and ⋅ represent the ceil (round toward positive infinity) and floor (round toward negative infinity) operation in MATLAB, respectively. ∠⋅ returns the phase angle of the operand, which lies between −π,π.


### 3.1.3. Coarse DOA Estimates without Ambiguities

For the purpose of disambiguation, unambiguous but coarse DOA estimates must be obtained as references to the ambiguous fine estimates in the previous subsection. From Figure4, by extracting the cumulant domain information of the first and the last 2Mu sensors in each subarray, the quarter-wavelength spatial invariance of the array geometry can be exploited to generate unambiguous coarse DOA estimates. This can be easily done by defining a permutation matrix
(24)Γ=I3⊗e1,I3⊗e2,…,I3⊗e2Mu+1,where ei signifies the ith column of the identity matrix I2Mu+1.Figure 4
Construction of 2 subspace matrix elements for coarse DOA estimation.From (19) and (24), we have
(25)UP=ΓUS=Asub∘DsT.By taking the first and the last3×2Mu rows of UP, we obtain the following two matrices.
(26)UP1=UP1:6Mu,:=A1∘DsT,UP2=UP4:end,:=A2∘DsT=A1∘DsΦcT,where A1 and A2 are the first and the last 2Mu rows of Asub, respectively. Φc=diagej2ω1,…,ej2ωK. Since UP1 is contributed by the cumulant domain data from the first 2Mu sensors in each subarray and UP2 is contributed by the cumulant domain data from the last 2Mu sensors in each subarray, unambiguous DOA estimates with quarter-wavelength spatial invariance can be found from Φc.From (26), we have UP2=UP1T−1ΦcT, which means that there is a rotational invariance between UP1 and UP2, that is, Ψc=T−1ΦcT=UP1#UP2. Therefore, the diagonal elements of Φc, Φck,k=ej2ωk=ej4πdsinθk/λ, correspond to the eigenvalues of Ψc. The unambiguous coarse DOA estimates of K-mixed FF and NF sources are given by
(27)θ^kref=arcsin∠Φck,k4πd/λ,k=1,…,K.θ^kref and θ^klθ may be paired by the method in [20].


### 3.1.4. Disambiguation

The coarse estimatesθ^kref could be used as references to disambiguate the cyclic ambiguities in θ^klθ. Mathematically, we obtain the disambiguated angle estimates as
(28)θ^k=arcsinμk+lθ0λ2ds,where the estimate of lθ0 is given by
(29)lθ0=argminlθθ^klθ−θ^kref.Note thatlθ0 is solved through searching over a few discrete points and the searching range is given by (22).


## 3.1.1. Cumulant Matrix for the DOA Estimation

Similar to the TSMUSIC algorithm in [14], we first construct a fourth-order cumulant matrix to estimate the DOAs of the radiating sources in this subsection. However, there are some differences between the two algorithms. For example, (1) the uniform linear array (ULA) configuration in TSMUSIC is generalized to a SLA in Figure 1, which promotes the estimation accuracy effectively and (2) the proposed method utilizes ESPRIT to generate the DOA estimation, which avoids the spectral search procedure in TSMUSIC and therefore reduces the computational complexity.According to the definition in [23], the fourth-order cumulant of the array outputs xmt, xnt, xit, and xjt can be written as
(9)cumxmt,xn∗t,xi∗t,xjt=∑k=1Kcskejpm−pn−pi−pjωk+pm2−pn2−pi2−pj2ϕk,where csk=cumskt,sk∗t,sk∗t,skt symbolizes the kurtosis of the kth source.We can construct a2M+1×2M+1 cumulant matrix C1, with its m,nth element being
(10)C1m,n=cumxm−M−1t,x−m+M+1∗t,xn−M−1∗t,x−n+M+1t=∑k=1Kcskej2pm−M−1ωk⋅e−j2pn−M−1ωk,m,n∈1,2M+1.Note thatC1 can be written in a matrix form
(11)C1=BC4sBH,where
(12)C4s=diagcs1,…,csk,…,csK,(13)B=bω1,…,bωk,…,bωK,bωk=ej2ωkP.From (1), B can be written as
(14)B=AsubD−1AsubAsubD1=Ds∘Asub,where
(15)D=e−j2Δω1⋱e−j2Δωk⋱e−j2ΔωK,Ds=d1,…,dk,…,dK,dk=e−j2Δωk,1,ej2ΔωkT,herein, Asub is the virtual steering matrix for each subarray in the cumulant domain, with the following form:
(16)Asub=a¯1,…,a¯k,…,a¯K,a¯K=e−j2Muωk,e−j2Mu−1ωk,…,1,…,ej2Mu−1ωk,ej2MuωkT.Therefore,C1 can be expressed as
(17)C1=Ds∘AsubC4sDs∘AsubH.The reason to rewrite it in this manner is that it has a similar form of dual-size invariance and therefore, ESPRIT could be applied to yield coarse and fine estimates of DOAs.


## 3.1.2. Fine DOA Estimates with Ambiguities

We firstly estimate the DOAs of the mixed FF and NF sources. By taking an eigendecomposition of (17), we have
(18)C1=USΛSUSH+UNΛNUNH,where ΛS and ΛN are diagonal matrices containing K large eigenvalues and 2M+1−K small eigenvalues, respectively. UN is the 2M+1×2M+1−K eigenvector matrix spanning the noise subspace of C1. US is the 2M+1×K signal subspace eigenvector matrix of C1, which can be expressed as
(19)US=Ds∘AsubT,where T is a nonsingular matrix.To obtain high-accuracy DOA estimates, we may form two matrices that are related by the extended intersubarray spacingds along the x-axis. As is illustrated in Figure 3, by taking the first and the last 22Mu+1 rows of US, we have the following two matrices.
(20)US1=US1:4Mu+2,:=D1∘AsubT,US2=US2Mu+1:end,:=D2∘AsubT=D1∘AsubΦfT,where D1 and D2 are the first and the last 2 rows of Ds and Φf=diagej2Δω1,…,ej2ΔωK. Since US1 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray−1 and subarray0 and US2 is contributed by the cumulant domain data from all 22Mu+1 sensors at subarray0 and subarray1, fine DOA estimates can be found from Φf due to the extended invariance relationship.Figure 3
Construction of 2 subspace matrix elements for fine DOA estimation.From (20), we have US2=US1T−1ΦfT, which means that there is a rotational invariance between US1 and US2, that is, Ψf=T−1ΦfT=US1#US2. Therefore, the diagonal elements of Φf, Φfk,k=ej2Δωk=ej4πdssinθk/λ, correspond to the eigenvalues of Ψf. As ds≫λ/4 and sinθk≤1, we can find a series of ambiguous DOA estimates
(21)θ^klθ=arcsinμk+lθλ2ds,k=1,…,K,(22)2dsλ−1−μk≤lθ≤2dsλ1−μk,(23)μk=∠Φfk,k4πds/λ,where lθ is an integer and ⋅ and ⋅ represent the ceil (round toward positive infinity) and floor (round toward negative infinity) operation in MATLAB, respectively. ∠⋅ returns the phase angle of the operand, which lies between −π,π.


## 3.1.3. Coarse DOA Estimates without Ambiguities

For the purpose of disambiguation, unambiguous but coarse DOA estimates must be obtained as references to the ambiguous fine estimates in the previous subsection. From Figure4, by extracting the cumulant domain information of the first and the last 2Mu sensors in each subarray, the quarter-wavelength spatial invariance of the array geometry can be exploited to generate unambiguous coarse DOA estimates. This can be easily done by defining a permutation matrix
(24)Γ=I3⊗e1,I3⊗e2,…,I3⊗e2Mu+1,where ei signifies the ith column of the identity matrix I2Mu+1.Figure 4
Construction of 2 subspace matrix elements for coarse DOA estimation.From (19) and (24), we have
(25)UP=ΓUS=Asub∘DsT.By taking the first and the last3×2Mu rows of UP, we obtain the following two matrices.
(26)UP1=UP1:6Mu,:=A1∘DsT,UP2=UP4:end,:=A2∘DsT=A1∘DsΦcT,where A1 and A2 are the first and the last 2Mu rows of Asub, respectively. Φc=diagej2ω1,…,ej2ωK. Since UP1 is contributed by the cumulant domain data from the first 2Mu sensors in each subarray and UP2 is contributed by the cumulant domain data from the last 2Mu sensors in each subarray, unambiguous DOA estimates with quarter-wavelength spatial invariance can be found from Φc.From (26), we have UP2=UP1T−1ΦcT, which means that there is a rotational invariance between UP1 and UP2, that is, Ψc=T−1ΦcT=UP1#UP2. Therefore, the diagonal elements of Φc, Φck,k=ej2ωk=ej4πdsinθk/λ, correspond to the eigenvalues of Ψc. The unambiguous coarse DOA estimates of K-mixed FF and NF sources are given by
(27)θ^kref=arcsin∠Φck,k4πd/λ,k=1,…,K.θ^kref and θ^klθ may be paired by the method in [20].


## 3.1.4. Disambiguation

The coarse estimatesθ^kref could be used as references to disambiguate the cyclic ambiguities in θ^klθ. Mathematically, we obtain the disambiguated angle estimates as
(28)θ^k=arcsinμk+lθ0λ2ds,where the estimate of lθ0 is given by
(29)lθ0=argminlθθ^klθ−θ^kref.Note thatlθ0 is solved through searching over a few discrete points and the searching range is given by (22).


## 3.2. Range Estimation

### 3.2.1. Cumulant Matrix for the Range Estimation

In order to derive the range estimates, another fourth-order cumulant matrix needs to be constructed to overcome the rank-deficient phenomenon described in [14]. The virtual steering matrix for range estimation in [5] has the following form
(30)A~=11⋯1ej2φ1ej2φ2⋯ej2φK⋮⋮⋯⋮ej2N−1φ1ej2N−1φ2⋯ej2N−1φK.If thekth source is in FF, the electrical angle ϕk will become 0 and the kth column of A~ will be 1,…,1T. When the number of FF sources is more than one, A~ will drop rank and thus the algorithm will fail to localize radiating sources.In this section, a new fourth-order cumulant matrixC2 is defined, with its m,nth element being
(31)C2m,n=cumx0t,x0∗t,xm−M−1t,xn−M−1∗t=∑k=1Kcskejpm−M−1ωk+pm−M−12ϕk⋅ejpn−M−1ωk+pn−M−12ϕk∗,m,n∈1,2M+1.Note thatC2 can be expressed as
(32)C2=cumx0t,x0∗t,Xt,X∗t=AC4sA∗,where A is the steering matrix in (6) and C4s is the kurtosis diagonal matrix of sources as defined in (12). When the kth source is in FF, the corresponding column vector aωk,ϕk becomes
(33)aωk,ϕk=ejωkp−M,…,ejωkp0,…,ejωkpMT.Therefore,A will still be of full column rank when there are multiple FF sources.


### 3.2.2. Unambiguous and Ambiguous Range Estimations

From (1) and (7), aωk,ϕk can be decoupled as
(34)aωk,ϕk=Vωkgϕk,Vωk=ejωk−Δ−Mu⋱ejωk−Δ+Muejωk−Mu⋱1⋰ejωkMuejωkΔ−Mu⋰ejωkΔ+Mu,(35)gϕk=ejϕkΔ+Mu2,…,ejϕkΔ2,…,ejϕkΔ−Mu2,ejϕkMu2,…,ejϕk12,ejϕk02T,where Vωk is a 2M+1×M+1 matrix which only contains the DOA information and gϕk is a M+1×1 column vector which only depends on ϕk.By taking an eigendecomposition ofC2, we can obtain a 2M+1×K signal subspace matrix ES and a 2M+1×2M+1−K noise subspace matrix EN, where ES is composed of the K principle eigenvectors of C2 and EN contains the remaining 2M+1−K eigenvectors.According to [4], the 2-D MUSIC spectrum for DOA and range estimation is given by
(36)Pω,ϕ=aHω,ϕENENHaω,ϕ.Substituting the disambiguated DOA estimatesθ^kk=1,2,…,K into the above equation, the 2-D spectrum search can be reduced to the K 1-D ones. Therefore, the estimates of ϕkk=1,2,…K are given by
(37)ϕ^k=minϕaHω^,ϕENENHaω^,ϕ=minϕgHϕVHω^ENENHVω^gϕ.Actually, (37) implies that gϕk is the eigenvector associated with the minimum eigenvalue of the Hermitian matrix VHω^kENENHVω^k. Through the eigendecomposition of VHω^kENENHVω^k, we can obtain an estimation of gϕk.In order to generate the coarse range estimates unambiguously,ϕk can be estimated as
(38)ϕ^k=∠gϕkMgϕkM+1,where gϕkn represents the nth element in gϕk. From (3), the coarse range estimate is given by
(39)r^kref=πd2λϕ^kcos2θ^k,where θ^k is the disambiguated DOA estimation. According to the definition of the Fresnel region, NF sources are located in the range of 0.62D3/λ,2D2/λ. As a result, when r^kref is in this region, the corresponding source is classified as a NF one; otherwise, it is regarded as a FF one.Similarly, through the extended aperture size between each subarray, fine range estimates with ambiguity can be generated. From (35), we have
(40)gϕkMu+1gϕkM+1=ejΔ2ϕk=ejπds2/λrkcos2θk.Therefore, we can find a series of ambiguous range estimates(41)r^klr=πds2νk+2lrπλcos2θ^k,(42)−νk2π≤lr≤ds22λr^kref−νk2π,(43)νk=∠gϕkMu+1gϕkM+1.


### 3.2.3. Disambiguation

In an analogous manner to that of Section3.1.4, the coarse estimates r^kref could be used as references to disambiguate the cyclic ambiguities in r^klr. Mathematically, the disambiguated range estimates are
(44)r^k=πds2νk+2lr0πλcos2θ^k,where the estimate of lr0 is given by
(45)lr0=argminlθr^klr−r^kref.Note thatlr0 is solved through searching over a few discrete points, and the searching range is given by (42).


## 3.2.1. Cumulant Matrix for the Range Estimation

In order to derive the range estimates, another fourth-order cumulant matrix needs to be constructed to overcome the rank-deficient phenomenon described in [14]. The virtual steering matrix for range estimation in [5] has the following form
(30)A~=11⋯1ej2φ1ej2φ2⋯ej2φK⋮⋮⋯⋮ej2N−1φ1ej2N−1φ2⋯ej2N−1φK.If thekth source is in FF, the electrical angle ϕk will become 0 and the kth column of A~ will be 1,…,1T. When the number of FF sources is more than one, A~ will drop rank and thus the algorithm will fail to localize radiating sources.In this section, a new fourth-order cumulant matrixC2 is defined, with its m,nth element being
(31)C2m,n=cumx0t,x0∗t,xm−M−1t,xn−M−1∗t=∑k=1Kcskejpm−M−1ωk+pm−M−12ϕk⋅ejpn−M−1ωk+pn−M−12ϕk∗,m,n∈1,2M+1.Note thatC2 can be expressed as
(32)C2=cumx0t,x0∗t,Xt,X∗t=AC4sA∗,where A is the steering matrix in (6) and C4s is the kurtosis diagonal matrix of sources as defined in (12). When the kth source is in FF, the corresponding column vector aωk,ϕk becomes
(33)aωk,ϕk=ejωkp−M,…,ejωkp0,…,ejωkpMT.Therefore,A will still be of full column rank when there are multiple FF sources.


## 3.2.2. Unambiguous and Ambiguous Range Estimations

From (1) and (7), aωk,ϕk can be decoupled as
(34)aωk,ϕk=Vωkgϕk,Vωk=ejωk−Δ−Mu⋱ejωk−Δ+Muejωk−Mu⋱1⋰ejωkMuejωkΔ−Mu⋰ejωkΔ+Mu,(35)gϕk=ejϕkΔ+Mu2,…,ejϕkΔ2,…,ejϕkΔ−Mu2,ejϕkMu2,…,ejϕk12,ejϕk02T,where Vωk is a 2M+1×M+1 matrix which only contains the DOA information and gϕk is a M+1×1 column vector which only depends on ϕk.By taking an eigendecomposition ofC2, we can obtain a 2M+1×K signal subspace matrix ES and a 2M+1×2M+1−K noise subspace matrix EN, where ES is composed of the K principle eigenvectors of C2 and EN contains the remaining 2M+1−K eigenvectors.According to [4], the 2-D MUSIC spectrum for DOA and range estimation is given by
(36)Pω,ϕ=aHω,ϕENENHaω,ϕ.Substituting the disambiguated DOA estimatesθ^kk=1,2,…,K into the above equation, the 2-D spectrum search can be reduced to the K 1-D ones. Therefore, the estimates of ϕkk=1,2,…K are given by
(37)ϕ^k=minϕaHω^,ϕENENHaω^,ϕ=minϕgHϕVHω^ENENHVω^gϕ.Actually, (37) implies that gϕk is the eigenvector associated with the minimum eigenvalue of the Hermitian matrix VHω^kENENHVω^k. Through the eigendecomposition of VHω^kENENHVω^k, we can obtain an estimation of gϕk.In order to generate the coarse range estimates unambiguously,ϕk can be estimated as
(38)ϕ^k=∠gϕkMgϕkM+1,where gϕkn represents the nth element in gϕk. From (3), the coarse range estimate is given by
(39)r^kref=πd2λϕ^kcos2θ^k,where θ^k is the disambiguated DOA estimation. According to the definition of the Fresnel region, NF sources are located in the range of 0.62D3/λ,2D2/λ. As a result, when r^kref is in this region, the corresponding source is classified as a NF one; otherwise, it is regarded as a FF one.Similarly, through the extended aperture size between each subarray, fine range estimates with ambiguity can be generated. From (35), we have
(40)gϕkMu+1gϕkM+1=ejΔ2ϕk=ejπds2/λrkcos2θk.Therefore, we can find a series of ambiguous range estimates(41)r^klr=πds2νk+2lrπλcos2θ^k,(42)−νk2π≤lr≤ds22λr^kref−νk2π,(43)νk=∠gϕkMu+1gϕkM+1.


## 3.2.3. Disambiguation

In an analogous manner to that of Section3.1.4, the coarse estimates r^kref could be used as references to disambiguate the cyclic ambiguities in r^klr. Mathematically, the disambiguated range estimates are
(44)r^k=πds2νk+2lr0πλcos2θ^k,where the estimate of lr0 is given by
(45)lr0=argminlθr^klr−r^kref.Note thatlr0 is solved through searching over a few discrete points, and the searching range is given by (42).


## 4. Discussion

In this section, we compare CDSSI with some recently developed mixed source localization algorithms, including TSMUSIC [14], the generalized ESPRIT-based (GESPRIT) algorithm [16], and MOMUSIC [19]. All of the four methods are analysed from the following aspects:


### 4.1. Array Aperture

Both TSMUSIC and GESPRIT employ an ULA which requires the intersensor spacing to be within a quarter wavelength. Therefore, with the same sensor number, SLA enjoys extended array aperture size, producing better estimation accuracy for the proposed algorithm. MOMUSIC also utilizes a special nested SLA in the development of the algorithm. However, according to the array model described in [19], the array aperture of MOMUSIC equals to when M1=2 and M2=2 (9 sensors in all), while the array aperture of CDSSI, with the same number of sensors, could be extended by spacing subarrays apart at tens of quarter wavelengths or more [20].


### 4.2. Maximum Number of Sources That Can Be Resolved

In CDSSI, the coarse and fine DOA estimates are derived from two6Mu+K matrices and two 22Mu+1×K matrices, respectively; therefore, the maximum number of sources that can be resolved is 4Mu+1. For MOMUSIC algorithm, a virtual ULA with 2M1+1M2+1+1 sensors can be constructed. However, it has a half aperture loss in the formation of a Toeplitz matrix, and therefore, it can resolve M1+1M2+1 sources at most. With a 2M+1 element ULA, TSMUSIC is capable of resolving up to 2M sources, while GESPRIT can resolve a maximum number of 2M−1 sources.


### 4.3. Computational Complexity

Only the major computation load is considered in this comparison, including construction of the cumulant matrices, eigenvalue decomposition (EVD), and spectral search. The searching steps for the angle parameterθ∈−90°,90° and the range parameter r∈0.62D3/λ,2D2/λ are denoted as θΔ and rΔ. Let L and 2M+1 symbolize the snapshot number and sensor number, respectively. For MOMUSIC, we assume that M1=M2=M, and therefore, the sensor number is 4M+1. TSMUSIC involves computing two fourth-order matrices with dimensions 2M+1×2M+1 and 4M+1×4M+1, two EVDs, and spectral search for DOA estimation. The GESPRIT algorithm requires the construction of two second-order covariance matrices, two EVDs, and spectral search for DOA and range estimation. MOMUSIC requires the construction of a fourth-order cumulant matrix and a second-order covariance matrix, two EVDs, and spectral search for DOA and range estimation. CDSSI does not need any spectral search operation but requires two fourth-order cumulant matrices with dimension 2M+1×2M+1, EVDs on C1, C2, Ψf, Ψc, and VHω^kENENHVω^k (K times). The computational complexity of the four methods is listed in Table 1.Table 1
Computational complexity of different algorithms.
Algorithm
Cumulant matrix
EVDs
Spectral search
CDSSSI
182M+12L
4/322M+13+2K3+KM+13
182M+12L
TSMUSIC
92M+12L+94M+12L
4/32M+13+4/34M+13
2⋅1802M+12/θΔ
GESPRIT
22M+12L
8/32M+13
K×2D2/λ−0.62D3/λ/rΔ×2M+12+2⋅1802M+12/θΔ
MOMUSIC
92M+12+1L+4M+12L
4/3M+12+13+4/34M+13
K×2D2/λ−0.62D3/λ/rΔ×4M+12+2⋅180M+12+12/θΔCompared with the second-order-based methods, CDSSI requires more computations in constructing the fourth-order cumulant matrix and EVDs. However, it avoids any 1-D or 2-D complicated spectral search. Since the search steps need to be dense enough for the spectral search-based algorithms to approach their theoretical bounds, the computation load of these algorithms will in turn increase dramatically.


### 4.4. Performance in Correlated Noise

Both TSMUSIC and the proposed algorithm are capable of suppressing additive Gaussian colored noise since they apply fourth-order cumulant in the whole estimation procedure, while GESPRIT and MOMUSIC, which rely on the second-order statistics, will degrade in the presence of spatially correlated noise.


### 4.5. Arc Length and First-Order Curvature

Furthermore, knowledge of the “manifold shape” not only is essential for the investigation of ambiguities and assessment of the detection-resolution capabilities of an array but it may also prove useful in developing new and more effective methods for its search process. We studied on the two array manifold properties, namely, arc length and first-order curvature, and analyse the accuracy and resolution capabilities of mixed sources [24].It has been shown that the response of sparse array towards a far-field/near-field source emitting narrowband spherical wavefront from azimuthθ and range r can be written as
(46)aθ,r=ejωP+ϕP2,where the associated parameter is ωk=−2πd/λsinθ, ϕk=πd2/λrkcos2θ, and P is the sensor position vector.
(47)P=p−M,p−M+1,…,p0,…,pM−1,pMT=−Δ−Mu,…,−Δ,…,−Δ+Mu⏟subarray−1,−Mu,…,0,…,Mu⏟subarray0,Δ−Mu,…,Δ,…,Δ+Mu⏟subarray1T.For the far-field source scenario, the range parameter approaches to∞ and the associated parameter ϕk becomes 0.The rate of change of arc lengthds/dθ and first-order curvature k1 of the θ parameter curves are as follows:
(48)ṡθ,r=dsdθ=−j2πdλcosθP+drsinθ⋅P2,ȧθ,r=−j2πdλcosθP+drsinθ⋅P2⋅aθ,r,k1θ,r=1ṡ2a¨−s¨ṡȧ,where ȧ implies differentiation with respect to θ and the arc length s is defined as
(49)s=∫0θdaθdθdθ.The Cramer-Rao lower bound under the assumption of spatially and temporally uncorrelated large number of snapshotsL≫1 may be written as
(50)CRBθ=σ22Lṡ2θ,r1−u1Hθ,raθ,r2aθ,r2−1,where the unit-norm tangent vector u1θ,r to the array manifold has been substituted for ȧθ,r/ṡθ,r.The expression of the asymptotic (L≫1) variance for a single emitter of unit power can be written as
(51)Varθ=1+σ2aθ,r2CRBθ.The variance asymptotically approaches the CRB for high SNR or largeN and has a similar dependence on ṡθ,r.Figures5, 6, and 7 show the Cramer-Rao lower bound capabilities of our proposed sparse array. Figures 5 and 6 indicate the variation of theoretical CRB with respect to the change of r. Figure 7 shows the theoretical CRB with respect to the change of snapshot. From Figures 5 and 6, we can see that the smaller the angle with the normal direction of array, the less the estimation error and when the range between the array and source is enlarged, the theoretical estimation performance becomes worse. And we can conclude that more sample data can improve estimation performance from Figure 7.Figure 5
Theoretical CRB with respect to the change ofθ and r.Figure 6
Asymptotic variance with respect to the change ofθ and r.Figure 7
Theoretical CRB with respect to the change of snapshot.In addition, the angle and range estimation accuracy is not only related to the intersubarray spacing of sparse but also related to the possibility of disambiguation. As the intersubarray spacing increases, the possibility of wrong disambiguation increases simultaneously. The MIE method [25] can be used to predict MSE performance; the angle and range estimation can be represented as
(52)Varu=Eu−u^2∣intervalerrorPintervalerror+Eu−u^2∣nointervalerrorPnointervalerror,where u and u^ denote the estimation value and true value of angle and range parameters, Pintervalerror is the possibility of right disambiguation, and Pnointervalerror is the possibility of wrong disambiguation.


## 4.1. Array Aperture

Both TSMUSIC and GESPRIT employ an ULA which requires the intersensor spacing to be within a quarter wavelength. Therefore, with the same sensor number, SLA enjoys extended array aperture size, producing better estimation accuracy for the proposed algorithm. MOMUSIC also utilizes a special nested SLA in the development of the algorithm. However, according to the array model described in [19], the array aperture of MOMUSIC equals to when M1=2 and M2=2 (9 sensors in all), while the array aperture of CDSSI, with the same number of sensors, could be extended by spacing subarrays apart at tens of quarter wavelengths or more [20].


## 4.2. Maximum Number of Sources That Can Be Resolved

In CDSSI, the coarse and fine DOA estimates are derived from two6Mu+K matrices and two 22Mu+1×K matrices, respectively; therefore, the maximum number of sources that can be resolved is 4Mu+1. For MOMUSIC algorithm, a virtual ULA with 2M1+1M2+1+1 sensors can be constructed. However, it has a half aperture loss in the formation of a Toeplitz matrix, and therefore, it can resolve M1+1M2+1 sources at most. With a 2M+1 element ULA, TSMUSIC is capable of resolving up to 2M sources, while GESPRIT can resolve a maximum number of 2M−1 sources.


## 4.3. Computational Complexity

Only the major computation load is considered in this comparison, including construction of the cumulant matrices, eigenvalue decomposition (EVD), and spectral search. The searching steps for the angle parameterθ∈−90°,90° and the range parameter r∈0.62D3/λ,2D2/λ are denoted as θΔ and rΔ. Let L and 2M+1 symbolize the snapshot number and sensor number, respectively. For MOMUSIC, we assume that M1=M2=M, and therefore, the sensor number is 4M+1. TSMUSIC involves computing two fourth-order matrices with dimensions 2M+1×2M+1 and 4M+1×4M+1, two EVDs, and spectral search for DOA estimation. The GESPRIT algorithm requires the construction of two second-order covariance matrices, two EVDs, and spectral search for DOA and range estimation. MOMUSIC requires the construction of a fourth-order cumulant matrix and a second-order covariance matrix, two EVDs, and spectral search for DOA and range estimation. CDSSI does not need any spectral search operation but requires two fourth-order cumulant matrices with dimension 2M+1×2M+1, EVDs on C1, C2, Ψf, Ψc, and VHω^kENENHVω^k (K times). The computational complexity of the four methods is listed in Table 1.Table 1
Computational complexity of different algorithms.
Algorithm
Cumulant matrix
EVDs
Spectral search
CDSSSI
182M+12L
4/322M+13+2K3+KM+13
182M+12L
TSMUSIC
92M+12L+94M+12L
4/32M+13+4/34M+13
2⋅1802M+12/θΔ
GESPRIT
22M+12L
8/32M+13
K×2D2/λ−0.62D3/λ/rΔ×2M+12+2⋅1802M+12/θΔ
MOMUSIC
92M+12+1L+4M+12L
4/3M+12+13+4/34M+13
K×2D2/λ−0.62D3/λ/rΔ×4M+12+2⋅180M+12+12/θΔCompared with the second-order-based methods, CDSSI requires more computations in constructing the fourth-order cumulant matrix and EVDs. However, it avoids any 1-D or 2-D complicated spectral search. Since the search steps need to be dense enough for the spectral search-based algorithms to approach their theoretical bounds, the computation load of these algorithms will in turn increase dramatically.


## 4.4. Performance in Correlated Noise

Both TSMUSIC and the proposed algorithm are capable of suppressing additive Gaussian colored noise since they apply fourth-order cumulant in the whole estimation procedure, while GESPRIT and MOMUSIC, which rely on the second-order statistics, will degrade in the presence of spatially correlated noise.


## 4.5. Arc Length and First-Order Curvature

Furthermore, knowledge of the “manifold shape” not only is essential for the investigation of ambiguities and assessment of the detection-resolution capabilities of an array but it may also prove useful in developing new and more effective methods for its search process. We studied on the two array manifold properties, namely, arc length and first-order curvature, and analyse the accuracy and resolution capabilities of mixed sources [24].It has been shown that the response of sparse array towards a far-field/near-field source emitting narrowband spherical wavefront from azimuthθ and range r can be written as
(46)aθ,r=ejωP+ϕP2,where the associated parameter is ωk=−2πd/λsinθ, ϕk=πd2/λrkcos2θ, and P is the sensor position vector.
(47)P=p−M,p−M+1,…,p0,…,pM−1,pMT=−Δ−Mu,…,−Δ,…,−Δ+Mu⏟subarray−1,−Mu,…,0,…,Mu⏟subarray0,Δ−Mu,…,Δ,…,Δ+Mu⏟subarray1T.For the far-field source scenario, the range parameter approaches to∞ and the associated parameter ϕk becomes 0.The rate of change of arc lengthds/dθ and first-order curvature k1 of the θ parameter curves are as follows:
(48)ṡθ,r=dsdθ=−j2πdλcosθP+drsinθ⋅P2,ȧθ,r=−j2πdλcosθP+drsinθ⋅P2⋅aθ,r,k1θ,r=1ṡ2a¨−s¨ṡȧ,where ȧ implies differentiation with respect to θ and the arc length s is defined as
(49)s=∫0θdaθdθdθ.The Cramer-Rao lower bound under the assumption of spatially and temporally uncorrelated large number of snapshotsL≫1 may be written as
(50)CRBθ=σ22Lṡ2θ,r1−u1Hθ,raθ,r2aθ,r2−1,where the unit-norm tangent vector u1θ,r to the array manifold has been substituted for ȧθ,r/ṡθ,r.The expression of the asymptotic (L≫1) variance for a single emitter of unit power can be written as
(51)Varθ=1+σ2aθ,r2CRBθ.The variance asymptotically approaches the CRB for high SNR or largeN and has a similar dependence on ṡθ,r.Figures5, 6, and 7 show the Cramer-Rao lower bound capabilities of our proposed sparse array. Figures 5 and 6 indicate the variation of theoretical CRB with respect to the change of r. Figure 7 shows the theoretical CRB with respect to the change of snapshot. From Figures 5 and 6, we can see that the smaller the angle with the normal direction of array, the less the estimation error and when the range between the array and source is enlarged, the theoretical estimation performance becomes worse. And we can conclude that more sample data can improve estimation performance from Figure 7.Figure 5
Theoretical CRB with respect to the change ofθ and r.Figure 6
Asymptotic variance with respect to the change ofθ and r.Figure 7
Theoretical CRB with respect to the change of snapshot.In addition, the angle and range estimation accuracy is not only related to the intersubarray spacing of sparse but also related to the possibility of disambiguation. As the intersubarray spacing increases, the possibility of wrong disambiguation increases simultaneously. The MIE method [25] can be used to predict MSE performance; the angle and range estimation can be represented as
(52)Varu=Eu−u^2∣intervalerrorPintervalerror+Eu−u^2∣nointervalerrorPnointervalerror,where u and u^ denote the estimation value and true value of angle and range parameters, Pintervalerror is the possibility of right disambiguation, and Pnointervalerror is the possibility of wrong disambiguation.


## 5. Simulation Results

In this section, numerical simulations are conducted to validate the performance of the proposed algorithm relative to TSMUSIC, GESPRIT, and MOMUSIC. In the following experiments, we consider a SLA composed of2M+1=9Mu=1 elements with ds=20d for the proposed algorithm, a quarter-wavelength-spaced ULA consisting of 15 elements for TSMUSIC and GESPRIT, and a sparse nested array with 9 elements M1=M2=2 for the MOMUSIC algorithm. NF sources are set to lie in a common Fresnel region of these algorithms. Moreover, the source signals are equipower, statistically independent, and of the form ejφk, where the phase φk is uniformly distributed between 0,2π. The performance is measured by the root mean squared error (RMSE) of 500 independent Monte Carlo trials. The RMSE is defined as
(53)RMSE=1500∑n=1500y^n,k−yk2,where yk stands for the DOA θk or the range rk and y^n,k denotes the estimation of yk in the nth trial.In the first experiment, we consider two equipower sources that are located atθ1=10°,r1=15λ and θ2=−10°,r2=+∞, that is, a mixed FF and NF scenario. The number of snapshots equals to 500, and the SNR varies from −10 dB to 20 dB in steps of 3 dB. The RMSEs of the four algorithms as a function of SNR are plotted in Figure 8. From these figures, one can observe that the coarse DOA and range estimates of the proposed algorithm have higher RMSEs due to the quarter-wavelength spatial shift invariance. However, after disambiguation, the fine estimates have superior estimation accuracy than those of the other three algorithms. Moreover, RMSEs of the DOA and range estimates decrease as the SNR increases.Figure 8
RMSEs of DOA and range estimates versus SNR.
(a)
(b)In the second experiment, we investigate the RMSEs of the four algorithms with the variation of the number of snapshots. The parameter settings are the same as those of the first experiment except that SNR is set equal to 10 dB and the number of snapshots varies from 100 to 10,000. From Figure9, it is obvious that, as a result of the extended aperture, CDSSI outperforms the other three algorithms in DOA and range estimation accuracy for all snapshot numbers. In addition, both the DOA and range estimation performance of all four algorithms improves as the snapshot number increases. This is because that larger sample support will produce better estimate of the covariance matrix for stationary data.Figure 9
RMSEs of DOA and range estimates versus snapshot number.
(a)
(b)In the third experiment, the scenario of two NF sources is investigated, with the source location parameters beingθ1=−10°,r1=10λ and θ1=10°,r1=20λ. The snapshot number is fixed at 500 and the SNR varies from −10 dB to 30 dB in steps of 5 dB. Figure 10 leads to a similar conclusion as in the first experiment that the proposed algorithm achieves the best performance owing to its extended aperture. Additionally, from the second figure, one can observe that the range estimation accuracy of the first source, which is closer to the array, is better than that of the second source. This result is consistent with the theoretical analysis developed in [6].Figure 10
RMSEs of DOA and range estimates versus SNR.
(a)
(b)In the fourth experiment, we study the dependence of DOA and range estimation accuracy upon the angular gap between two NF sources.θ1 varies form −80° to 80° in steps of 10°, while θ2=0°, r1=10λ, and r2=20λ, with the snapshot number and the SNR equal to 500 and 10 dB, respectively. The RMSEs of the angle and range estimation are shown in Figure 11. From these figures, one can see that CDSSI is superior to the other three algorithms in the performance of both DOA and range estimation for all θ values. In addition, there is a similarity in Figure 11 that the range estimation accuracy of the first source decreases with the growth of θ1, while that of the second source remains almost constant. This phenomenon can be explained as follows: in all the four algorithms, range estimates are based on the estimation of DOAs and the DOA estimation errors are propagated to the subsequent range estimation.Figure 11
RMSEs of DOA and range estimates versus DOA gap of two NF sources.
(a)
(b)In the fifth experiment, the range parameter of the first sourcer1 varies from 10λ to 20λ in steps of λ, while the range of the second source r2 is fixed at 20λ. The DOAs of the two sources are θ1=−10° and θ2=10°. Let the snapshot number and SNR be 500 and 10 dB, respectively. The influence of r1 on the DOA and range estimation is shown in Figure 12. From the first figure, it is obvious that the DOA estimation performance is insensitive to the change of the range parameter, since the DOA estimation is decoupled with the range estimation in all these algorithms. However, the GESPRIT algorithm behaves abnormally when the range of the first source approaches to 20λ (the range parameter of the second source). This is because when the two sources are symmetrical with respect to the broadside, that is, θ1=−θ2 and r1=r2, the GESPRIT algorithm will generate image sources, which are misidentified as real ones. Moreover, the second figure reveals that (1) the range estimation performance of the second source is hardly affected by the range variation of the first source and (2) the range estimation accuracy of the first source (closer to the array) is superior to that of the second one, which corroborates the theoretical analysis in [6]. Finally, from Figure 12, it is clear that the proposed algorithm outperforms the other three methods in both DOA and range estimation for all range parameters.Figure 12
RMSEs of DOA and range estimates versus range gap of two NF sources.
(a)
(b)In the last experiment, two FF sources are considered with their location parameters beingθ1=10°,r1=+∞ and θ2=60°,r2=+∞. The snapshot number is fixed at 500 and the SNR varies from −10 dB to 20 dB in steps of 3 dB. Note that the intersubarray spacing ds is extended from 20d to 40d without introducing additional sensors. The RMSEs of the DOA estimates with the variation of SNR are shown in Figure 13. From this figure, it is seen that all of the four algorithms are still effective in the multiple FF source scenarios and CDSSI outperforms the other three in the performance of DOAs estimation. Also, one can observe that as the intersubarray spacing becomes larger, the accuracy of the proposed method turns better compared with the results demonstrated in Figure 2.Figure 13
RMSEs of DOA estimates versus SNR.


## 6. Conclusion

In this paper, an efficient and high-performance algorithm is proposed for the mixed far-field and near-field source localization problems. Based on a sparse linear array of dual-size spatial invariance, the proposed algorithm can offer enhanced accuracy due to the extended aperture size. Moreover, the proposed method has lower computational complexity because it does not require any 1-D or 2-D spectral search. According to the simulations, the proposed algorithm outperforms the conventional ones in the performance of both angle and range estimation.



---

*Source: 3237167-2018-03-26.xml*