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

. 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.

Over these years, coprime array [21] attracts much attention. It can e ectively increase the degrees of freedom (DOFs) [22,23], relieve the mutual coupling (MC) e ects [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]. Speci cally, 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. is 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 nal DOA estimates.
is algorithm results in signi cantly 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. ese 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]. e 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. en, 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. e 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.
e remaining parts of this paper are organized as follows: in Section 2, 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. D m (·) 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.

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 has M sensors with d 1 � (Nλ/2), where λ represents the wavelength. e other subarray is with N sensors and the interelement spacing is denoted as d 2 � (Mλ/2). So the total number of the sensors is denoted as T UCLA � M + N − 1. Figure 1 is an example of UCLA configuration where M � 3 and N � 4.
Assume that there are K uncorrelated far-field narrowband signals s k (t) impinging on the UCLA from distinct angles where t ∈ [1, L] and L represents the number of snapshots. e 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. e received signal of the array can be denoted as follows: where , a(θ 2 ), . . . , a(θ K )] is the direction matrix and the steering vector is defined by a(θ k ) � [a 1 (θ k ) T , a 2 (θ k ) T ] T , n(t) � [n 1 (t) T , n 2 (t) T ] T is the additive white Gaussian noise with zero mean and variance σ 2 n . And, the noise signal is independent of the signal resources. And represents the directional matrix and the corresponding steering vector is denoted as and the corresponding steering vector is represented as Practically, the covariance matrix is approximately computed with L snapshots [7].
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: It exists a full rank matrix Γ ∈ C (K×K) [7] to make (5) hold.

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 matrix A 1 , we can get A 11 and A 12 , which contain the first K rows and (M − K) rows, respectively.
For the subarray 1, we first partition the steering matrix A 1 as follows: where A 11 ∈ C K×K represents the matrix contains the first K rows of A 1 and A 12 ∈ C (M− K)×K stands for the matrix with the remaining (M − K) rows of A 1 , respectively. Assume that A 11 is a full rank matrix, then we can obtain A 12 by the following equation: where P 1c is the propagator method of the subarray 1. And P 1c ∈ C (M− K)×K [14]. en, we define the following equation: where I 1K is a unit matrix of I 1K ∈ C K×K . So we have the following equation: en, we partition the matrix of P 1 and can get P 1a and P 1b Where P 1a and P 1b denote the first (M − 1) rows and last (M − 1) rows of P 1 , respectively. And, Ξ ξ and Ξ ζ , respectively, represent the last row and the first row of P 1 . en, we partition A 1 by the following equation: where A 1a and A 1b denote the first (M − 1) rows and last (M − 1) rows of A 1 , respectively. Σ ξ represents the last row and Σ ζ is the first row of A 1 . en, we can get the following equation: According to (10), we have the following equation: Mathematical Problems in Engineering So it has the following equation: en, we have the following equation: where P + 1a � (P H 1a P 1a ) − 1 P H 1a gives the pseudoinverse of P 1a and Φ 1r is a diagonal matrix which is denoted as follows: We define the following equation: Because A 11 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,k (k � 1, 2, . . . , K) of subarray 1, which is denoted as follows: where angle(·) means angle function. By the similar conduction, we process the subarray 2. Separate the directional matrix A 2 into two parts and we can get A 21 and A 22 , which contain the first K rows and (N − K) rows, respectively. e steering matrix of A 2 is separated as follows: where A 21 ∈ C K×K represents the matrix contains the first K rows of A 2 and A 22 ∈ C (N− K)×K represents the matrix with the remaining (N − K) rows of A 2 , respectively. Assume that A 21 is a full rank matrix, then we can obtain A 22 by the following equation: where P 2c is the propagator method of the subarray 1. And P 2c ∈ C (N− K)×K . en, we define the following equation: where I 2K is a unit matrix of I 2K ∈ C K×K . Similar to equation 15 we have the following equation: en, we partition the matrix of P 2 and can get P 2a and P 2b where P 2a and P 2b denote the first (N − 1) rows and last (N − 1) rows of P 2 , respectively. And, Ξ 2ξ and Ξ 2ζ , respectively, represent the last row and the first row of P 2 . en, we partition A 2 by the following equation: where A 2a and A 2b denote the first (N − 1) rows and last (N − 1) rows of A 2 , respectively. Σ 2ξ represents the last row and Σ 2ζ is the first row of A 2 . en, we can get the following equation: According to (22), we have the following equation: en, we can get the following equation: en, we have the following equation: where P + 2a � (P H 2a P 2a ) − 1 P H 2a gives the pseudo-inverse of P 2a and Φ 2r is a diagonal matrix which is denoted as follows: Φ 2r � diag(e jπMsinθ 1 , e jπMsinθ 2 , . . . , e jπMsinθ K ) ∈ C K×K .
We define the following equation: Because A 21 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: where k � 1, 2, . . . , K and angle(·) is the angle function.

Ambiguity Elimination Based on Coprime Property.
In this part, according to the obtained angles, we first recover all the estimates. en, we eliminate the ambiguity problem based on the coprime property.
It is known that there exists 2kπ(k ∈ Z) between the real and ambiguous angles for the sinusoid function [28].
where Q i ∈ Z, d 1 � N d, d 2 � M d, θ i,am means the ambiguous angle of the subarray i. It has the following equation: According to the variation range of θ, it is indicated that where M and N are integers [27]. en, we have the following equation: 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. Figure 2 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.
Due to the coprime property of M and N, there only exists Q 1 � Q 2 � 0 which makes the equation (32) satisfied.
By equation (36), we can get the initial DOA estimates.
where R s is the covariance matrix of the signals and σ 2 n I (M+N)×(M+N) denotes the power of noise. From equations (3) and (37), it has the following equation: Due to the orthogonality of the signal and noise subspace, it exists I � U s U H s + U n U H n . So equation (38) can be rewritten as follows: en, we can get the following equation: As U s � AΓ and U H s U s � I (M+N)×(M+N) , we have the following equation: However, the noise exists. To solve this problem, establish a fitting relationship to compute the matrix Γ.
which can make the equation (5) hold. By utilizing the least square (LS) criterion, we can get the following equation: Incorporate (36) and (37), then we have the following equation: 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 Mathematical Problems in Engineering reconstruct the steering matrix and search within a small sector. In this way, complexity gets significantly decreased.
According to the obtained initial DOA estimates (45) en, the angle θ 1 can be computed by the following equation: It can be noted the searching region is 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 . en, we keep in K ] unchanged and elaborate a new directional matrix A (2) .
Here, θ is the angle that we will estimate in the following step.
By equation (48), we can obtain the estimate of θ 2 by PAS And we employ θ 2 to establish A (3) , a θ 1 , a θ 2 , a(θ) 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 By the similar method, we reconstruct the new direc- a θ 1 , a θ 2 , . . . , a θ K− 1 , a(θ) .
en. we can attain the estimate of θ K by the following equation: Here, the angle searches within a small region . Due to transforming the multi-dimensional GAS of SF into initialization based 1D PAS, the computational complexity is significantly reduced.  Mathematical Problems in Engineering

Detailed
Steps.
e detailed steps of the proposed method are (Algorithm 1) as follows:

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 ( 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. e 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.

Cramer-Rao Bound.
Here, we derive the CRB [37] of the UCLA.
Elaborate the manifold matrix of the UCLA as follows: where A 2p denotes the rows from the second one to the last one of the A 2 .

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 Figure 5. (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 Figure 2. (3) By employing the obtained initial DOA estimates, the GAS is transformed into PAS. In this way, the Step 1: compute the covariance matrix R according to equation (2) Step 2: operate the EVD of R 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 θ complexity has an effective reduction, which can be seen in Section 4. (4) e 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 Section 5.

Simulations
In the simulation section, the root mean square error (RMSE) is used as the performance comparison metric, which is defined as follows: 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.  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.

Comparison of Different Arrays with the Same Algorithm.
e 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.

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. Figures 11 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.

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.

Data Availability
e data used to support the findings of this study are included within the article.

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