# A Computationally Efficient Algorithm for DOA Estimation with Unfolded Coprime Linear Array

Authors: Gong Pan; Yin Huafei
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we investigate the direction of arrival (DOA) estimation problem with unfolded coprime linear array (UCLA) and propose a low computational complexity signal-subspace fitting (SF) algorithm. SF algorithm is able to achieve excellent DOA estimation performance while it requires global angular search (GAS). Especially in the several source signals situation, expensive complexity cost causes. To decrease computational complexity, we propose an initialized based SF (ISF) algorithm, which involves the several one dimensional (1D) partial angular search (PAS) instead of the multidimensional GAS. Consequently, the complexity is significantly decreased. Due to the full utilization of the array aperture, the proposed method in UCLA can attain better performance than general CLA (GCLA). In addition, as the SF is attractive in practical application, the proposed ISF algorithm lowers the computational cost, while achieving almost approximate estimation performance as traditional SF and noise subspace fitting (NF). Moreover, numerical simulations are provided and verify the effectiveness and the superiority of the proposed algorithm for the UCLA.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation is one of the fundamental issues for the array signal processing scenery and has been applied in engineering fields, including sonar, radar, navigation, and wireless communications [1–6]. In the past decades, many subspace based algorithms have been proposed [7–10], like multiple signals classification (MUSIC) based algorithms [7–10], and estimation of signal parameters via rotational invariance techniques (ESPRIT) [11–13]. These are subspace based algorithms. The propagator method (PM) [14, 15] can reduce the computational complexity by employing a linear partition operation instead of eigenvalue decomposition (EVD). These algorithms were initially designed for uniform array [16–19]. Nevertheless, for the conventional uniform arrays, the interelement spacing is required to be no larger than half-wavelength. As a result, the phase ambiguity problem can be avoided [20].Over these years, coprime array [21] attracts much attention. It can effectively increase the degrees of freedom (DOFs) [22, 23], relieve the mutual coupling (MC) effects [16, 24], and improve the angle estimation performance. Because of these advantages, the coprime array is widely used in wireless communication systems and radar location [25, 26]. Specifically, a general coprime linear array (GCLA) incorporates two sparse uniform linear subarrays with M and N sensors, where M and N are coprime integers. And, the interelement spacing of these two subarrays are Nλ/2 and Mλ/2, respectively. And, λ means the wavelength. This design concept breaks the conventional half-wavelength and can achieve the higher angle resolution compared with the classic uniform array in the same conditions.In these years, various algorithms have been proposed for DOA estimation with GCLA. Zhou proposed a total spectral search (TSS) algorithm in [27]. By combining the DOA estimates of two subarrays to attain the final DOA estimates. This algorithm results in significantly computational complexity because of the global angular search (GAS). A partial spectral searching algorithm [28], which investigates the linear relationship to obtain all estimates, was proposed. Moreover, it transforms the GAS into partial sector one. These algorithms treat the array separately, so they only employ the auto-information of two subarrays. An efficient method which can resolve the ambiguity in DOA estimation was proposed in [29]. The method offers good generalization and robustness in resolving the ambiguity problem. It achieves full degrees of freedom (DOF) with reduced complexity. An ambiguity-free algorithm via utilizing the total matrix information, such as auto-covariance information and mutual covariance information, was proposed in [30]. However, it involves high computational complexity. Along with pursuing the high resolution and DOA estimation performance, the computational complexity is also a challenging but promising task [31, 32].It is known that subspace fitting techniques [33, 34] are popular in array signal processing [35, 36]. Compared with maximum likelihood [37], signal subspace fitting (SF) and noise subspace fitting (NF) [33] algorithms obtain the similar angle estimation performance [37], while these algorithms involve high computational complexity due to the GAS, especially in the multiple signals situation. A successive scheme of SF has proposed in [38], which incorporates the coprime linear array and SF to decrease the complexity. To further expand the array aperture, we link the SF into unfolded coprime linear array (UCLA), which enlarges the array aperture and we transform the multidimensional searching into several one dimensional (1D) searching. Moreover, we replace the GAS by partial angular search (PAS). Specifically, by PM, we can attain the initial DOAs of two subarrays. And, we recover all estimates and obtain the unique initial DOA estimates according to coprime property. Then, we employ the initial estimates to reconstruct the steering matrix and transform the multidimensional search into several 1D one. Consequently, computational complexity cost can be significantly decreased. Meanwhile, via the initial estimates, we replace the GAS by PAS. The proposed ISF can acquire better DOA estimation performance with UCLA than that with GCLA due to the larger array aperture. And, it acquires similar DOA estimation performance compared with SF and NF, while ISF has the lowest complexity. Moreover, Cramer–Rao Bound (CRB) is presented as a theoretical lower bound [39]. Finally, the effectiveness and superiority of the proposed ISF algorithm for the UCLA is demonstrated by the numerical simulations.Specifically, we summarize the main contributions of this paper as follows:(1)
We integrate the UCLA with the subspace fitting method which can obtain a larger array aperture compared with GCLA. Simulations verify that the proposed algorithm with UCLA can realize more excellent estimation performance than GCLA.(2)
We propose an initialization based algorithm for DOA estimation, which can effectively decrease the complexity of the classic SF algorithm. By utilizing PM to initialize and obtain coarse estimation, and operating fine searching among a small sector, so we can achieve lower complexity.(3)
We demonstrate that the proposed algorithm can achieve the approximately the same DOA estimation performance as the classical SF and NF algorithms. And, the proposed algorithm outperforms the classic PM algorithm in DOA estimation performance.The remaining parts of this paper are organized as follows: in Section2, we elaborate the UCLA geometry and signal model. Subsequently, the proposed algorithm is introduced in Section 3. Complexity analysis and advantages are given in Section 4. Numerical simulations are provided in Sections 5 and 6 conclude this paper.Notations: we utilize lower-case (upper-case) bold characters as vectors (matrices). And, we use⋅T and ⋅H to represent the transpose and the conjugate transpose, respectively. ⊙ and ⊗ represent the Khatri–Rao product and Kronecker product, respectively. diag⋅ denotes a diagonal matrix which employs the elements of the matrix to be its diagonal elements. E⋅ represents statistical expectation. min⋅ is getting the minimum element. Dm⋅ is a diagonal matrix that the m-th row of the matrix is employed. angle⋅ and arctan⋅ denote phase operator and the arctangent function, respectively.


## 2. Signal Model

In this paper, we employ an unfolded coprime linear array (UCLA) configuration which is able to further enlarge the array aperture and promote DOA estimation performance.An UCLA configuration incorporates two uniform linear subarrays. One subarray hasM sensors with d1=Nλ/2, where λ represents the wavelength. The other subarray is with N sensors and the interelement spacing is denoted as d2=Mλ/2. So the total number of the sensors is denoted as TUCLA=M+N−1. Figure 1 is an example of UCLA configuration where M=3 and N=4.Figure 1
Structure of unfolded coprime linear array.Assume that there areK uncorrelated far-field narrow-band signals skt impinging on the UCLA from distinct angles where t∈1,L and L represents the number of snapshots. The angles are denoted as Θ=θ1,θ2,…,θK, where θk∈0,π/2, k=1,2,…,K. Here, we assume the number of sources K is known. The received signal of the array can be denoted as follows:(1)xt=Ast+nt=x1tx2t=A1A2st+n1tn2t,where A=A1T,A2TT=aθ1,aθ2,…,aθK is the direction matrix and the steering vector is defined by aθk=a1θkT,a2θkTT, nt=n1tT,n2tTT is the additive white Gaussian noise with zero mean and variance σn2. And, the noise signal is independent of the signal resources. And st=s1t,s2t,…,sKtT denotes the signal vector, where t=1,2,…,L, L means the number of snapshots. A1=a1θ1,a1θ2,…,a1θK represents the directional matrix and the corresponding steering vector is denoted as a1θk=ejM−1Nπsinθk,ejM−2Nπsinθk,…,1Tk=1,2,…,K. And, the directional matrix of subarray 2 is denoted as A2=a2θ1,a2θ2,…,a2θK and the corresponding steering vector is represented as a2θk=e−jMπsinθk,…,e−jN−1MπsinθkT.Practically, the covariance matrix is approximately computed withL snapshots [7].(2)R^=1L∑t=1LxtxHt.Then, perform eigenvalue decomposition [7].(3)R^=U^sD^sU^sH+U^nD^nU^nH,where D^s and D^n are the diagonal matrices composed of the largest K eigenvalues of R^ and the diagonal matrix containing the remaining eigenvalues, respectively. And, U^s denotes the signal subspace which consists of the eigenvectors corresponding to the largest K eigenvalues. U^n is the noise subspace including the rest eigenvectors.In the noise-free case, we can get the following equation:(4)spaceUs=spaceA.It exists a full rank matrixΓ∈ℂK×K [7] to make (5) hold.(5)U^s=AΓ.


## 3. Proposed Method for DOA Estimation

### 3.1. Initialization Processing

In this subsection, we first utilize subarray 1 to introduce the proposed algorithm. And we can operate the subarray 2 by the similar method.By partitioning the directional matrixA1, we can get A11 and A12, which contain the first K rows and M−K rows, respectively.For the subarray 1, we first partition the steering matrixA1 as follows:(6)A1=A11A12,where A11∈ℂK×K represents the matrix contains the first K rows of A1 and A12∈ℂM−K×K stands for the matrix with the remaining M−K rows of A1, respectively.Assume thatA11 is a full rank matrix, then we can obtain A12 by the following equation:(7)A12=P1cA11,where P1c is the propagator method of the subarray 1. And P1c∈ℂM−K×K [14].Then, we define the following equation:(8)P1=I1KP1c,where I1K is a unit matrix of I1K∈ℂK×K.So we have the following equation:(9)P1A11=I1KA11P1cA11=A11A12=A1.Then, we partition the matrix ofP1 and can get P1a and P1b(10)P1=P1aΞξ=ΞζP1b,Where P1a and P1b denote the first M−1 rows and last M−1 rows of P1, respectively. And, Ξξ and Ξζ, respectively, represent the last row and the first row of P1.Then, we partitionA1 by the following equation:(11)A1=A1aΣξ=ΣζA1b,where A1a and A1b denote the first M−1 rows and last M−1 rows of A1, respectively. Σξ represents the last row and Σζ is the first row of A1.Then, we can get the following equation:(12)P1A11=P1aΞξA11=ΞζP1bA11,A1=A11A12=A1aΣξ=ΣζA1b.According to (10), we have the following equation:(13)P1aA11=A1aP1bA11=A1b.So it has the following equation:(14)P1aP1bA11=A1aA1b=A1aA1aΦ1r.Then, we have the following equation:(15)P1a+P1b=A11Φ1rA11−1,where P1a+=P1aHP1a−1P1aH gives the pseudoinverse of P1a and Φ1r is a diagonal matrix which is denoted as follows:Φ1r=diagejπNsinθ1,ejπNsinθ2,…,ejπNsinθK∈ℂK×K.We define the following equation:(16)Ψ1r=P1a+P1b.BecauseA11 is a full rank matrix, so Ψ1r is the similar transformation of Φ1r.AsΦ1r is a diagonal matrix of eigenvalues, Ψ1r and Φ1r possess the same eigenvalues. As a result, operate eigenvalues decomposition of Ψ1r, then we can obtain the diagonal elements δ1,k. And, we can get the initial DOA estimates sinθ⟶1,kk=1,2,…,K of subarray 1, which is denoted as follows:(17)sinθ⟶1,k=angleδ1,kNπ,where angle⋅ means angle function.By the similar conduction, we process the subarray 2.Separate the directional matrixA2 into two parts and we can get A21 and A22, which contain the first K rows and N−K rows, respectively.The steering matrix ofA2 is separated as follows:(18)A2=A21A22,where A21∈ℂK×K represents the matrix contains the first K rows of A2 and A22∈ℂN−K×K represents the matrix with the remaining N−K rows of A2, respectively.Assume thatA21 is a full rank matrix, then we can obtain A22 by the following equation:(19)A22=P2cA21,where P2c is the propagator method of the subarray 1. And P2c∈ℂN−K×K.Then, we define the following equation:(20)P2=I2KP2c,where I2K is a unit matrix of I2K∈ℂK×K.Similar to equation15 we have the following equation:(21)P2A21=I2KA21P2cA21=A21A22=A2.Then, we partition the matrix ofP2 and can get P2a and P2b(22)P2=P2aΞ2ξ=Ξ2ζP2b,where P2a and P2b denote the first N−1 rows and last N−1 rows of P2, respectively. And, Ξ2ξ and Ξ2ζ, respectively, represent the last row and the first row of P2.Then, we partitionA2 by the following equation:(23)A2=A2aΣ2ξ=Σ2ζA2b,where A2a and A2b denote the first N−1 rows and last N−1 rows of A2, respectively. Σ2ξ represents the last row and Σ2ζ is the first row of A2.Then, we can get the following equation:(24)P2A21=P2aΞξA11=ΞζP2bA21,A2=A21A22=A2aΣξ=ΣζA2b.According to (22), we have the following equation:(25)P1aA11=A1aP1bA11=A1b.Then, we can get the following equation:(26)P2aP2bA21=A2aA2b=A2aA2aΦ2r.Then, we have the following equation:(27)P2a+P2b=A21Φ2rA21−1,where P2a+=P2aHP2a−1P2aH gives the pseudo-inverse of P2a and Φ2r is a diagonal matrix which is denoted as follows:Φ2r=diagejπMsinθ1,ejπMsinθ2,…,ejπMsinθK∈ℂK×K.We define the following equation:(28)Ψ2r=P2a+P2b.BecauseA21 is a full rank matrix, so Ψ2r is the similar transformation of Φ2r.AsΦ2r is a diagonal matrix of eigenvalues, Ψ2r and Φ2r possess the same eigenvalues. As a result, operate eigenvalues decomposition of Ψ2r, then we can obtain the diagonal elements δ2,k. And, we can get the initial DOA estimates sinθ^2,k of subarray 2, which is denoted as follows:(29)sinθ^2,k=angleδ2,k/Mπ,where k=1,2,…,K and angle⋅ is the angle function.


### 3.2. Ambiguity Elimination Based on Coprime Property

In this part, according to the obtained angles, we first recover all the estimates. Then, we eliminate the ambiguity problem based on the coprime property.It is known that there exists2kπk∈ℤ between the real and ambiguous angles for the sinusoid function [28].(30)2πdisinθ^iλ−2πdisinθ^i,amλ=2Qiπ,where Qi∈ℤ, d1=Nd,d2=Md, θi,am means the ambiguous angle of the subarray i. It has the following equation:(31)sinθ⟶1,kλ−sinθ^1,am=2Q1N,sinθ⟶2,kλ−sinθ^2,am=2Q2M.According to the variation range ofθ, it is indicated that Q1∈−N−1,N−1 and Q2∈−M−1,M−1, where M and N are integers [27].Then, we have the following equation:(32)2Q1N=2Q2M.It is known that the interelement spacing of a uniform linear array is no larger than half wave length to avoid the phase ambiguity. As a result, no phase ambiguity problem results in. But the coprime array, due to the element spacing larger than half wavelength, arises phase ambiguity.To illustrate the phase ambiguity problem, we provide the simulation about the coprime array. Figure2 depicts the DOA estimation with the three different element spacing, where there is one signal θ=25° arrives at the array. And it can be noticed that there are ambiguous angles when d=3λ/2 and d=5λ/2.Figure 2
DOA estimation with the varying element spacing.Due to the coprime property ofM and N, there only exists Q1=Q2=0 which makes the equation (32) satisfied.Via equations (33) and (34), all the DOA estimates are obtained.(33)θ^M,k=arcsinsinθ⟶1,k−2Q1N,(34)θ^N,k=arcsinsinθ⟶2,k−2Q2M,where k=1,2,…,K, θ^M=θ^M,1,θ^M,2,…,θ^M,k and θ^N=θ^N,1,θ^N,2,…,θ^N,k.Practically, considering that noise exists, to attain the overlapped angle estimation is always difficult. Consequently, we replace searching the overlap by finding the nearest angles fromθ^Mξ and θ^Nζ , which contain all the estimates of two subarrays, respectively.(35)minθ^Mξ,θ^Nζθ^Mξ−θ^Nζξ=1,2,…,N,ζ=1,2,…,M.By equation (36), we can get the initial DOA estimates.(36)θ^kini=θ^Mξ+θ^Nζ2k=1,2,…,K.


### 3.3. Initialization Based Algorithm

Via equation (37), we get the covariance matrix [7].(37)R=ARsAH+σn2IM+N×M+N,where Rs is the covariance matrix of the signals and σn2IM+N×M+N denotes the power of noise.From equations (3) and (37), it has the following equation:(38)ARsAH+σn2IM+N×M+N=UsDsUsH+σn2UnUnH.Due to the orthogonality of the signal and noise subspace, it existsI=UsUsH+UnUnH. So equation (38) can be rewritten as follows:(39)ARsAH+σn2IM+N×M+N=UsDsUsH+σn2IM+N×M+N−UsUsH.Then, we can get the following equation:(40)ARsAH+σn2UsUsH=UsDsUsH.AsUs=AΓ and UsHUs=IM+N×M+N, we have the following equation:(41)Γ=RsAHUsDs−σn2IM+N×M+N−1.However, the noise exists. To solve this problem, establish a fitting relationship to compute the matrixΓ.(42)θ,Γ^=minU^s−A^Γ^F2,which can make the equation (5) hold.By utilizing the least square (LS) criterion, we can get the following equation:(43)Γ^=A^HA^−1A^HU^s=A^+U^s.Incorporate (36) and (37), then we have the following equation:(44)θ=minU^s−A^A^+U^sF2=mintrI−A^A^HA^−1A^HU^sU^sH=maxtrA^A^HA^−1A^HU^sU^sH.When there are numerical signals, the problem of equation (44) is becoming a multidimensional SF problem. Consequently, it will have a higher computational cost. In view of this, we utilize the initialization based method to reconstruct the steering matrix and search within a small sector. In this way, complexity gets significantly decreased.According to the obtained initial DOA estimatesθ^in=θ^1in,θ^2in,…,θ^Kin, the new manifold matrix A^1 is obtained.(45)A^1=aθ,aθ^2in,…,aθ^Kin.Then, the angleθ1 can be computed by the following equation:(46)θ^1=argminθ∈θ^1in−Δθ,θ^1in+ΔθU^s−A^1A^1U^sF2.It can be noted the searching region isθ∈θ^1in−Δθ,θ^1in+Δθ , where Δθ is a tiny value. In this way, we can get the more accurate DOA estimate of θ1.From equation (46), we can obtained θ^1. Then, we keep θ^1,θ^3in,…,θ^Kin unchanged and elaborate a new directional matrix A^2.(47)A^2=aθ^1,aθ,aθ^3in,…,aθ^Kin.Here,θ is the angle that we will estimate in the following step.By equation (48), we can obtain the estimate of θ2 by PAS within θ∈θ^2in−Δθ,θ^2in+Δθ.(48)θ^2=argminθ∈θ^2in−Δθ,θ^2in+ΔθU^s−A^2A^2+U^sF2.Similarly, keepθ^1,θ^2,θ^4in,…,θ^Kin unchanged. And we employ θ^2 to establish A^3,(49)A^3=aθ^1,aθ^2,aθ,aθ^4in,…,aθ^Kin.It is noted thatθ^1 and θ^2 is estimated via equations (46) and (48), and θ is the goal that we are to estimate in the next step.Via equation (50), we can get the more accurate DOA estimate of θ3 within a small searching region θ∈θ^3in−Δθ,θ^3in+Δθ.(50)θ^3=argminθ∈θ^3in−Δθ,θ^3in+ΔθU^s−A^3A^3U^sF2.By the similar method, we reconstruct the new directional matrixA^K via using θ^,1θ^2,…,θ^K−1.(51)A^K=aθ^1,aθ^2,…,aθ^K−1,aθ.Then. we can attain the estimate ofθK by the following equation:(52)θ^K=argminθ∈θ^Kini−Δθ,θ^Kini+ΔθU^s−A^KA^K+U^sF2.Here, the angle searches within a small region θ∈θ^Kini−Δθ,θ^Kini+Δθ.Due to transforming the multi-dimensional GAS of SF into initialization based 1D PAS, the computational complexity is significantly reduced.


### 3.4. Detailed Steps

The detailed steps of the proposed method are (Algorithm1) as follows:Algorithm 1: The details of the proposed ISF algorithm.
Step 1: compute the covariance matrixR^ according to equation (2)Step 2: operate the EVD ofR^ and get the signal subspace by equation (3)Step 3: via propagator method, obtain the initial angle estimation and recover all the DOA estimates by equations (33) and (34)Step 4: ambiguity problem is solved via the coprime property and initial estimatesθ^kin,k=1,2,…,K are achievedStep 5: compute fine DOA estimates according to equation (52)


## 3.1. Initialization Processing

In this subsection, we first utilize subarray 1 to introduce the proposed algorithm. And we can operate the subarray 2 by the similar method.By partitioning the directional matrixA1, we can get A11 and A12, which contain the first K rows and M−K rows, respectively.For the subarray 1, we first partition the steering matrixA1 as follows:(6)A1=A11A12,where A11∈ℂK×K represents the matrix contains the first K rows of A1 and A12∈ℂM−K×K stands for the matrix with the remaining M−K rows of A1, respectively.Assume thatA11 is a full rank matrix, then we can obtain A12 by the following equation:(7)A12=P1cA11,where P1c is the propagator method of the subarray 1. And P1c∈ℂM−K×K [14].Then, we define the following equation:(8)P1=I1KP1c,where I1K is a unit matrix of I1K∈ℂK×K.So we have the following equation:(9)P1A11=I1KA11P1cA11=A11A12=A1.Then, we partition the matrix ofP1 and can get P1a and P1b(10)P1=P1aΞξ=ΞζP1b,Where P1a and P1b denote the first M−1 rows and last M−1 rows of P1, respectively. And, Ξξ and Ξζ, respectively, represent the last row and the first row of P1.Then, we partitionA1 by the following equation:(11)A1=A1aΣξ=ΣζA1b,where A1a and A1b denote the first M−1 rows and last M−1 rows of A1, respectively. Σξ represents the last row and Σζ is the first row of A1.Then, we can get the following equation:(12)P1A11=P1aΞξA11=ΞζP1bA11,A1=A11A12=A1aΣξ=ΣζA1b.According to (10), we have the following equation:(13)P1aA11=A1aP1bA11=A1b.So it has the following equation:(14)P1aP1bA11=A1aA1b=A1aA1aΦ1r.Then, we have the following equation:(15)P1a+P1b=A11Φ1rA11−1,where P1a+=P1aHP1a−1P1aH gives the pseudoinverse of P1a and Φ1r is a diagonal matrix which is denoted as follows:Φ1r=diagejπNsinθ1,ejπNsinθ2,…,ejπNsinθK∈ℂK×K.We define the following equation:(16)Ψ1r=P1a+P1b.BecauseA11 is a full rank matrix, so Ψ1r is the similar transformation of Φ1r.AsΦ1r is a diagonal matrix of eigenvalues, Ψ1r and Φ1r possess the same eigenvalues. As a result, operate eigenvalues decomposition of Ψ1r, then we can obtain the diagonal elements δ1,k. And, we can get the initial DOA estimates sinθ⟶1,kk=1,2,…,K of subarray 1, which is denoted as follows:(17)sinθ⟶1,k=angleδ1,kNπ,where angle⋅ means angle function.By the similar conduction, we process the subarray 2.Separate the directional matrixA2 into two parts and we can get A21 and A22, which contain the first K rows and N−K rows, respectively.The steering matrix ofA2 is separated as follows:(18)A2=A21A22,where A21∈ℂK×K represents the matrix contains the first K rows of A2 and A22∈ℂN−K×K represents the matrix with the remaining N−K rows of A2, respectively.Assume thatA21 is a full rank matrix, then we can obtain A22 by the following equation:(19)A22=P2cA21,where P2c is the propagator method of the subarray 1. And P2c∈ℂN−K×K.Then, we define the following equation:(20)P2=I2KP2c,where I2K is a unit matrix of I2K∈ℂK×K.Similar to equation15 we have the following equation:(21)P2A21=I2KA21P2cA21=A21A22=A2.Then, we partition the matrix ofP2 and can get P2a and P2b(22)P2=P2aΞ2ξ=Ξ2ζP2b,where P2a and P2b denote the first N−1 rows and last N−1 rows of P2, respectively. And, Ξ2ξ and Ξ2ζ, respectively, represent the last row and the first row of P2.Then, we partitionA2 by the following equation:(23)A2=A2aΣ2ξ=Σ2ζA2b,where A2a and A2b denote the first N−1 rows and last N−1 rows of A2, respectively. Σ2ξ represents the last row and Σ2ζ is the first row of A2.Then, we can get the following equation:(24)P2A21=P2aΞξA11=ΞζP2bA21,A2=A21A22=A2aΣξ=ΣζA2b.According to (22), we have the following equation:(25)P1aA11=A1aP1bA11=A1b.Then, we can get the following equation:(26)P2aP2bA21=A2aA2b=A2aA2aΦ2r.Then, we have the following equation:(27)P2a+P2b=A21Φ2rA21−1,where P2a+=P2aHP2a−1P2aH gives the pseudo-inverse of P2a and Φ2r is a diagonal matrix which is denoted as follows:Φ2r=diagejπMsinθ1,ejπMsinθ2,…,ejπMsinθK∈ℂK×K.We define the following equation:(28)Ψ2r=P2a+P2b.BecauseA21 is a full rank matrix, so Ψ2r is the similar transformation of Φ2r.AsΦ2r is a diagonal matrix of eigenvalues, Ψ2r and Φ2r possess the same eigenvalues. As a result, operate eigenvalues decomposition of Ψ2r, then we can obtain the diagonal elements δ2,k. And, we can get the initial DOA estimates sinθ^2,k of subarray 2, which is denoted as follows:(29)sinθ^2,k=angleδ2,k/Mπ,where k=1,2,…,K and angle⋅ is the angle function.


## 3.2. Ambiguity Elimination Based on Coprime Property

In this part, according to the obtained angles, we first recover all the estimates. Then, we eliminate the ambiguity problem based on the coprime property.It is known that there exists2kπk∈ℤ between the real and ambiguous angles for the sinusoid function [28].(30)2πdisinθ^iλ−2πdisinθ^i,amλ=2Qiπ,where Qi∈ℤ, d1=Nd,d2=Md, θi,am means the ambiguous angle of the subarray i. It has the following equation:(31)sinθ⟶1,kλ−sinθ^1,am=2Q1N,sinθ⟶2,kλ−sinθ^2,am=2Q2M.According to the variation range ofθ, it is indicated that Q1∈−N−1,N−1 and Q2∈−M−1,M−1, where M and N are integers [27].Then, we have the following equation:(32)2Q1N=2Q2M.It is known that the interelement spacing of a uniform linear array is no larger than half wave length to avoid the phase ambiguity. As a result, no phase ambiguity problem results in. But the coprime array, due to the element spacing larger than half wavelength, arises phase ambiguity.To illustrate the phase ambiguity problem, we provide the simulation about the coprime array. Figure2 depicts the DOA estimation with the three different element spacing, where there is one signal θ=25° arrives at the array. And it can be noticed that there are ambiguous angles when d=3λ/2 and d=5λ/2.Figure 2
DOA estimation with the varying element spacing.Due to the coprime property ofM and N, there only exists Q1=Q2=0 which makes the equation (32) satisfied.Via equations (33) and (34), all the DOA estimates are obtained.(33)θ^M,k=arcsinsinθ⟶1,k−2Q1N,(34)θ^N,k=arcsinsinθ⟶2,k−2Q2M,where k=1,2,…,K, θ^M=θ^M,1,θ^M,2,…,θ^M,k and θ^N=θ^N,1,θ^N,2,…,θ^N,k.Practically, considering that noise exists, to attain the overlapped angle estimation is always difficult. Consequently, we replace searching the overlap by finding the nearest angles fromθ^Mξ and θ^Nζ , which contain all the estimates of two subarrays, respectively.(35)minθ^Mξ,θ^Nζθ^Mξ−θ^Nζξ=1,2,…,N,ζ=1,2,…,M.By equation (36), we can get the initial DOA estimates.(36)θ^kini=θ^Mξ+θ^Nζ2k=1,2,…,K.


## 3.3. Initialization Based Algorithm

Via equation (37), we get the covariance matrix [7].(37)R=ARsAH+σn2IM+N×M+N,where Rs is the covariance matrix of the signals and σn2IM+N×M+N denotes the power of noise.From equations (3) and (37), it has the following equation:(38)ARsAH+σn2IM+N×M+N=UsDsUsH+σn2UnUnH.Due to the orthogonality of the signal and noise subspace, it existsI=UsUsH+UnUnH. So equation (38) can be rewritten as follows:(39)ARsAH+σn2IM+N×M+N=UsDsUsH+σn2IM+N×M+N−UsUsH.Then, we can get the following equation:(40)ARsAH+σn2UsUsH=UsDsUsH.AsUs=AΓ and UsHUs=IM+N×M+N, we have the following equation:(41)Γ=RsAHUsDs−σn2IM+N×M+N−1.However, the noise exists. To solve this problem, establish a fitting relationship to compute the matrixΓ.(42)θ,Γ^=minU^s−A^Γ^F2,which can make the equation (5) hold.By utilizing the least square (LS) criterion, we can get the following equation:(43)Γ^=A^HA^−1A^HU^s=A^+U^s.Incorporate (36) and (37), then we have the following equation:(44)θ=minU^s−A^A^+U^sF2=mintrI−A^A^HA^−1A^HU^sU^sH=maxtrA^A^HA^−1A^HU^sU^sH.When there are numerical signals, the problem of equation (44) is becoming a multidimensional SF problem. Consequently, it will have a higher computational cost. In view of this, we utilize the initialization based method to reconstruct the steering matrix and search within a small sector. In this way, complexity gets significantly decreased.According to the obtained initial DOA estimatesθ^in=θ^1in,θ^2in,…,θ^Kin, the new manifold matrix A^1 is obtained.(45)A^1=aθ,aθ^2in,…,aθ^Kin.Then, the angleθ1 can be computed by the following equation:(46)θ^1=argminθ∈θ^1in−Δθ,θ^1in+ΔθU^s−A^1A^1U^sF2.It can be noted the searching region isθ∈θ^1in−Δθ,θ^1in+Δθ , where Δθ is a tiny value. In this way, we can get the more accurate DOA estimate of θ1.From equation (46), we can obtained θ^1. Then, we keep θ^1,θ^3in,…,θ^Kin unchanged and elaborate a new directional matrix A^2.(47)A^2=aθ^1,aθ,aθ^3in,…,aθ^Kin.Here,θ is the angle that we will estimate in the following step.By equation (48), we can obtain the estimate of θ2 by PAS within θ∈θ^2in−Δθ,θ^2in+Δθ.(48)θ^2=argminθ∈θ^2in−Δθ,θ^2in+ΔθU^s−A^2A^2+U^sF2.Similarly, keepθ^1,θ^2,θ^4in,…,θ^Kin unchanged. And we employ θ^2 to establish A^3,(49)A^3=aθ^1,aθ^2,aθ,aθ^4in,…,aθ^Kin.It is noted thatθ^1 and θ^2 is estimated via equations (46) and (48), and θ is the goal that we are to estimate in the next step.Via equation (50), we can get the more accurate DOA estimate of θ3 within a small searching region θ∈θ^3in−Δθ,θ^3in+Δθ.(50)θ^3=argminθ∈θ^3in−Δθ,θ^3in+ΔθU^s−A^3A^3U^sF2.By the similar method, we reconstruct the new directional matrixA^K via using θ^,1θ^2,…,θ^K−1.(51)A^K=aθ^1,aθ^2,…,aθ^K−1,aθ.Then. we can attain the estimate ofθK by the following equation:(52)θ^K=argminθ∈θ^Kini−Δθ,θ^Kini+ΔθU^s−A^KA^K+U^sF2.Here, the angle searches within a small region θ∈θ^Kini−Δθ,θ^Kini+Δθ.Due to transforming the multi-dimensional GAS of SF into initialization based 1D PAS, the computational complexity is significantly reduced.


## 3.4. Detailed Steps

The detailed steps of the proposed method are (Algorithm1) as follows:Algorithm 1: The details of the proposed ISF algorithm.
Step 1: compute the covariance matrixR^ according to equation (2)Step 2: operate the EVD ofR^ and get the signal subspace by equation (3)Step 3: via propagator method, obtain the initial angle estimation and recover all the DOA estimates by equations (33) and (34)Step 4: ambiguity problem is solved via the coprime property and initial estimatesθ^kin,k=1,2,…,K are achievedStep 5: compute fine DOA estimates according to equation (52)


## 4. Discussions

### 4.1. Complexity

In this part, we give the computational complexity comparison results of the proposed ISF algorithm, SF [2], NF [2], and TSS [27]. For ISF, it has the complexity of M2+N2L+M3+N3+2M−12+N−12K+M−13+N−13+M−1+N−1K2+Π14KM2+N2+M3+N3 where Π1=K⋅2Δ/ds means the search times and ds=0.001 is the search step, Δ denotes a tiny search value. Moreover, we provide the computational complexity comparison of the different algorithms in Table 1, including SF, NF, and TSS. The comparison of the computational complexity versus number of snapshots and sensors are illustrated in Figures 3 and 4, where M=3,N=4,K=3,ds=0.001 and N=4,5,7,8, respectively. As the proposed method transforms the GAS into PAS and searches over a small sector, it shows clearly that its complexity is much lower than SF, NF, and TSS. Figure 5 depicts the complexity comparison versus the search step. It is seen that ISF can significantly relieve the computational complexity burden.Table 1
Comparison of the computational complexity.
AlgorithmsComputational complexityISFOM2+N2L+M3+N3+2M−12+N−12K+M−13+N−13+M−1+N−1K2+G14KM2+N2+M3+N3SFOM2+N2L+M3+N3+G24KM2+N2+M3+N3NFOM2+N2L+M3+N3+G2MM−KK+NN−KKTSSOM2+N2L+M3+N3+G3MM−K+NN−KFigure 3
Complexity versus the number of snapshots.Figure 4
Complexity versus element number.Figure 5
Complexity versus search step number.


### 4.2. Cramer-Rao Bound

Here, we derive the CRB [37] of the UCLA.Elaborate the manifold matrix of the UCLA as follows:(53)At=A1A2p,where A2p denotes the rows from the second one to the last one of the A2.(54)CRB=σn22LReDHI−AtAtHAt−1AtHD⊕Rs−1,where Rs=1/L∑t=1LstsHt, D=∂at,1/∂θ1,∂at,2/∂θ2,…,∂at,K/∂θK, ⊕ means the Hadamard operation. And at,k is the Kth column of At.


### 4.3. Advantages

We give the advantages of the proposed ISF algorithm in the following:(1)
We incorporate the signal subspace fitting method into UCLA, which can achieve the more superior performance than GCLA due to the larger array aperture. It is seen in Figure5.(2)
When there are multiple signals, the proposed ISF transforms the conventional multi-dimensional search into several 1D search, which can remarkably decrease the computational complexity. It is seen in Figure2.(3)
By employing the obtained initial DOA estimates, the GAS is transformed into PAS. In this way, the complexity has an effective reduction, which can be seen in Section4.(4)
The proposed ISF is able to attain similar DOA estimation performance as traditional SF and NF algorithms and outperform ESPRIT and PM, which is seen in Section5.


## 4.1. Complexity

In this part, we give the computational complexity comparison results of the proposed ISF algorithm, SF [2], NF [2], and TSS [27]. For ISF, it has the complexity of M2+N2L+M3+N3+2M−12+N−12K+M−13+N−13+M−1+N−1K2+Π14KM2+N2+M3+N3 where Π1=K⋅2Δ/ds means the search times and ds=0.001 is the search step, Δ denotes a tiny search value. Moreover, we provide the computational complexity comparison of the different algorithms in Table 1, including SF, NF, and TSS. The comparison of the computational complexity versus number of snapshots and sensors are illustrated in Figures 3 and 4, where M=3,N=4,K=3,ds=0.001 and N=4,5,7,8, respectively. As the proposed method transforms the GAS into PAS and searches over a small sector, it shows clearly that its complexity is much lower than SF, NF, and TSS. Figure 5 depicts the complexity comparison versus the search step. It is seen that ISF can significantly relieve the computational complexity burden.Table 1
Comparison of the computational complexity.
AlgorithmsComputational complexityISFOM2+N2L+M3+N3+2M−12+N−12K+M−13+N−13+M−1+N−1K2+G14KM2+N2+M3+N3SFOM2+N2L+M3+N3+G24KM2+N2+M3+N3NFOM2+N2L+M3+N3+G2MM−KK+NN−KKTSSOM2+N2L+M3+N3+G3MM−K+NN−KFigure 3
Complexity versus the number of snapshots.Figure 4
Complexity versus element number.Figure 5
Complexity versus search step number.


## 4.2. Cramer-Rao Bound

Here, we derive the CRB [37] of the UCLA.Elaborate the manifold matrix of the UCLA as follows:(53)At=A1A2p,where A2p denotes the rows from the second one to the last one of the A2.(54)CRB=σn22LReDHI−AtAtHAt−1AtHD⊕Rs−1,where Rs=1/L∑t=1LstsHt, D=∂at,1/∂θ1,∂at,2/∂θ2,…,∂at,K/∂θK, ⊕ means the Hadamard operation. And at,k is the Kth column of At.


## 4.3. Advantages

We give the advantages of the proposed ISF algorithm in the following:(1)
We incorporate the signal subspace fitting method into UCLA, which can achieve the more superior performance than GCLA due to the larger array aperture. It is seen in Figure5.(2)
When there are multiple signals, the proposed ISF transforms the conventional multi-dimensional search into several 1D search, which can remarkably decrease the computational complexity. It is seen in Figure2.(3)
By employing the obtained initial DOA estimates, the GAS is transformed into PAS. In this way, the complexity has an effective reduction, which can be seen in Section4.(4)
The proposed ISF is able to attain similar DOA estimation performance as traditional SF and NF algorithms and outperform ESPRIT and PM, which is seen in Section5.


## 5. Simulations

In the simulation section, the root mean square error (RMSE) is used as the performance comparison metric, which is defined as follows:(55)RMSE=∑p=1Q∑k=1Kθ^k,p−θk2PK,where P is the number of Monte Carlo simulations, θ^k,p stands for the estimate of the p-th trial for the k-th theoretical angle θk. And, in this paper, we set P = 1000.


### 5.1. Scattering Figure of the Proposed ISF with UCLA

The scattering figure of the proposed ISF algorithm with UCLA for three distant sourcesθ=10∘,30∘,50∘ is presented in Figure 6, where M=3,N=4, L=200, SNR = 5 dB. And, we define the search step and the tiny searching restrain as ds=0.001 and Δ=0.5. It is shown clearly that the proposed ISF algorithm detects the source signals successfully.Figure 6
Scattering figure via ISF algorithm.


### 5.2. Comparison of Different Arrays with the Same Algorithm

The RMSE comparison versus SNR and snapshots with different configurations, including UCLA and GCLA, for two sourcesθ1,θ2=25∘,45∘ is given in Figures 7 and 8 by the same algorithm. It is defined that L=200 and SNR = 5 dB, respectively. From these two figures, we can notice that the UCLA is able to obtain the lower CRB and better DOA estimation performance than the GCLA. Moreover, the proposed ISF algorithm can attain the better DOA estimation performance with the UCLA than that with the GCLA because of the extension of the array aperture.Figure 7
RMSE performance versus SNR for different configurations.Figure 8
RMSE performance versus the number of snapshots for different configurations.


### 5.3. Comparison of Different Algorithms with the UCLA

In this subsection, the RMSE comparison of the proposed ISF algorithm, SF [33], NF [33], TSS [27], S-SF [38], ESPRIT [11], and PM [14] versus SNR and the number of snapshots is given in Figures 9 and 10, where M=3,N=4 and θ1,θ2=25∘,45∘ It is defined that L=200 and SNR = 5 dB, respectively. From these two figures, we can notice that ISF can achieve nearly similar estimation performance as SF, NF, and TSS but with the lower complexity due to the initialization operation to decrease the complexity which is verified in Figure 2. What’s more, ISF performs the better DOA estimation than ESPRIT and PM.Figure 9
RMSE performance versus SNR.Figure 10
RMSE performance versus the number of snapshots.


### 5.4. RMSE with Different Snapshots and SNR

Figures11 and 12 compare the estimation performance with a different number of snapshots and SNR, respectively. It shows clearly that the performance of angle estimation becomes better with the number of snapshots and SNR increasing.Figure 11
RMSE performance with different (L).Figure 12
RMSE performance with different SNR.


### 5.5. Estimation Probability Comparison of Different Algorithms

Figures13 and 14 depict the estimation probability versus the number of SNR and snapshots of the proposed ISF algorithm, SF [33], NF [33], TSS [27], S-SF [37], ESPRIT [11], and PM [14]. Suppose two closely located targets impinging on the arrays, where SNR=5dB,K=2,θ1,θ2=20°,21°. The two sources can be resolved if θ−θ^<θ1−θ2/2 where θ=θ1,θ2, θ^=θ^1,θ^2 [40]. We can clearly see that the proposed ISF algorithm performs the almost the same estimation probability than SF, NF, and TSS. It can be also inferred that ISF outperforms the ESPRIT and PM algorithms.Figure 13
Estimation probability versus SNR.Figure 14
Estimation probability versus snapshots.


## 5.1. Scattering Figure of the Proposed ISF with UCLA

The scattering figure of the proposed ISF algorithm with UCLA for three distant sourcesθ=10∘,30∘,50∘ is presented in Figure 6, where M=3,N=4, L=200, SNR = 5 dB. And, we define the search step and the tiny searching restrain as ds=0.001 and Δ=0.5. It is shown clearly that the proposed ISF algorithm detects the source signals successfully.Figure 6
Scattering figure via ISF algorithm.


## 5.2. Comparison of Different Arrays with the Same Algorithm

The RMSE comparison versus SNR and snapshots with different configurations, including UCLA and GCLA, for two sourcesθ1,θ2=25∘,45∘ is given in Figures 7 and 8 by the same algorithm. It is defined that L=200 and SNR = 5 dB, respectively. From these two figures, we can notice that the UCLA is able to obtain the lower CRB and better DOA estimation performance than the GCLA. Moreover, the proposed ISF algorithm can attain the better DOA estimation performance with the UCLA than that with the GCLA because of the extension of the array aperture.Figure 7
RMSE performance versus SNR for different configurations.Figure 8
RMSE performance versus the number of snapshots for different configurations.


## 5.3. Comparison of Different Algorithms with the UCLA

In this subsection, the RMSE comparison of the proposed ISF algorithm, SF [33], NF [33], TSS [27], S-SF [38], ESPRIT [11], and PM [14] versus SNR and the number of snapshots is given in Figures 9 and 10, where M=3,N=4 and θ1,θ2=25∘,45∘ It is defined that L=200 and SNR = 5 dB, respectively. From these two figures, we can notice that ISF can achieve nearly similar estimation performance as SF, NF, and TSS but with the lower complexity due to the initialization operation to decrease the complexity which is verified in Figure 2. What’s more, ISF performs the better DOA estimation than ESPRIT and PM.Figure 9
RMSE performance versus SNR.Figure 10
RMSE performance versus the number of snapshots.


## 5.4. RMSE with Different Snapshots and SNR

Figures11 and 12 compare the estimation performance with a different number of snapshots and SNR, respectively. It shows clearly that the performance of angle estimation becomes better with the number of snapshots and SNR increasing.Figure 11
RMSE performance with different (L).Figure 12
RMSE performance with different SNR.


## 5.5. Estimation Probability Comparison of Different Algorithms

Figures13 and 14 depict the estimation probability versus the number of SNR and snapshots of the proposed ISF algorithm, SF [33], NF [33], TSS [27], S-SF [37], ESPRIT [11], and PM [14]. Suppose two closely located targets impinging on the arrays, where SNR=5dB,K=2,θ1,θ2=20°,21°. The two sources can be resolved if θ−θ^<θ1−θ2/2 where θ=θ1,θ2, θ^=θ^1,θ^2 [40]. We can clearly see that the proposed ISF algorithm performs the almost the same estimation probability than SF, NF, and TSS. It can be also inferred that ISF outperforms the ESPRIT and PM algorithms.Figure 13
Estimation probability versus SNR.Figure 14
Estimation probability versus snapshots.


## 6. Conclusions

In this paper, we propose an ISF algorithm for DOA estimation with UCLA and verify that UCLA behaves the better DOA estimation performance than GCLA due to the larger array aperture. In the multiple signals scenery, the classic SF needs severe computational complexity cost due to the multidimensional GAS. To solve this problem, we transform the multi-dimensional search into several 1D one. In addition, GAS is changed to be PAS. Specifically, the propagator method is employed to obtain the initial DOA estimation. By initialization, we can transform the multidimensional GAS into several 1D partial one. As a result, the complexity is significantly reduced. CRB is presented and the simulations verify the effectiveness of the proposed algorithm.



---

*Source: 5429488-2022-09-05.xml*