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.
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 and sensors, where and are coprime integers. And, the interelement spacing of these two subarrays are and , 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 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 and to represent the transpose and the conjugate transpose, respectively. and represent the Khatri–Rao product and Kronecker product, respectively. denotes a diagonal matrix which employs the elements of the matrix to be its diagonal elements. represents statistical expectation. is getting the minimum element. is a diagonal matrix that the m-th row of the matrix is employed. and 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 has sensors with , where represents the wavelength. The other subarray is with sensors and the interelement spacing is denoted as . So the total number of the sensors is denoted as . Figure 1 is an example of UCLA configuration where and .

Assume that there are uncorrelated far-field narrow-band signals impinging on the UCLA from distinct angles where and represents the number of snapshots. The angles are denoted as , where , . Here, we assume the number of sources is known. The received signal of the array can be denoted as follows:where is the direction matrix and the steering vector is defined by , is the additive white Gaussian noise with zero mean and variance . And, the noise signal is independent of the signal resources. And denotes the signal vector, where , L means the number of snapshots. represents the directional matrix and the corresponding steering vector is denoted as . And, the directional matrix of subarray 2 is denoted as and the corresponding steering vector is represented as .
Practically, the covariance matrix is approximately computed with snapshots [7].
Then, perform eigenvalue decomposition [7].where and are the diagonal matrices composed of the largest eigenvalues of and the diagonal matrix containing the remaining eigenvalues, respectively. And, denotes the signal subspace which consists of the eigenvectors corresponding to the largest K eigenvalues. 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 [7] to make (5) hold.
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 matrix , we can get and , which contain the first rows and rows, respectively.
For the subarray 1, we first partition the steering matrix as follows:where represents the matrix contains the first rows of and stands for the matrix with the remaining rows of , respectively.
Assume that is a full rank matrix, then we can obtain by the following equation:where is the propagator method of the subarray 1. And [14].
Then, we define the following equation:where is a unit matrix of .
So we have the following equation:
Then, we partition the matrix of and can get and Where and denote the first rows and last rows of , respectively. And, and , respectively, represent the last row and the first row of .
Then, we partition by the following equation:where and denote the first rows and last rows of , respectively. represents the last row and is the first row of .
Then, we can get the following equation:
According to (10), we have the following equation:
So it has the following equation:
Then, we have the following equation:where gives the pseudoinverse of and is a diagonal matrix which is denoted as follows:.
We define the following equation:
Because is a full rank matrix, so is the similar transformation of .
As is a diagonal matrix of eigenvalues, and possess the same eigenvalues. As a result, operate eigenvalues decomposition of , then we can obtain the diagonal elements . And, we can get the initial DOA estimates of subarray 1, which is denoted as follows:where means angle function.
By the similar conduction, we process the subarray 2.
Separate the directional matrix into two parts and we can get and , which contain the first rows and rows, respectively.
The steering matrix of is separated as follows:where represents the matrix contains the first rows of and represents the matrix with the remaining rows of , respectively.
Assume that is a full rank matrix, then we can obtain by the following equation:where is the propagator method of the subarray 1. And .
Then, we define the following equation:where is a unit matrix of .
Similar to equation 15 we have the following equation:
Then, we partition the matrix of and can get and where and denote the first rows and last rows of , respectively. And, and , respectively, represent the last row and the first row of .
Then, we partition by the following equation:where and denote the first rows and last rows of , respectively. represents the last row and is the first row of .
Then, we can get the following equation:
According to (22), we have the following equation:
Then, we can get the following equation:
Then, we have the following equation:where gives the pseudo-inverse of and is a diagonal matrix which is denoted as follows:.
We define the following equation:
Because is a full rank matrix, so is the similar transformation of .
As is a diagonal matrix of eigenvalues, and possess the same eigenvalues. As a result, operate eigenvalues decomposition of , then we can obtain the diagonal elements . And, we can get the initial DOA estimates of subarray 2, which is denoted as follows:where and 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 exists between the real and ambiguous angles for the sinusoid function [28].where , , means the ambiguous angle of the subarray i. It has the following equation:
According to the variation range of , it is indicated that and , where and are integers [27].
Then, 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 arrives at the array. And it can be noticed that there are ambiguous angles when and .

Due to the coprime property of and , there only exists which makes the equation (32) satisfied.
Via equations (33) and (34), all the DOA estimates are obtained.where , and .
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 and , which contain all the estimates of two subarrays, respectively.
By equation (36), we can get the initial DOA estimates.
3.3. Initialization Based Algorithm
Via equation (37), we get the covariance matrix [7].where is the covariance matrix of the signals and 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 . So equation (38) can be rewritten as follows:
Then, we can get the following equation:
As and , 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 reconstruct the steering matrix and search within a small sector. In this way, complexity gets significantly decreased.
According to the obtained initial DOA estimates , the new manifold matrix is obtained.
Then, the angle 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 .
From equation (46), we can obtained . Then, we keep unchanged and elaborate a new directional matrix .
Here, is the angle that we will estimate in the following step.
By equation (48), we can obtain the estimate of by PAS within .
Similarly, keep unchanged. And we employ to establish ,
It is noted that and 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 within a small searching region .
By the similar method, we reconstruct the new directional matrix via using .
Then. we can attain the estimate of 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.
3.4. Detailed Steps
The detailed steps of the proposed method are (Algorithm 1) as follows:
|
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 where means the search times and 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 and 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.



4.2. Cramer-Rao Bound
Here, we derive the CRB [37] of the UCLA.
Elaborate the manifold matrix of the UCLA as follows:where denotes the rows from the second one to the last one of the .where , , means the Hadamard operation. And is the Kth column of .
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 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 complexity has an effective reduction, which can be seen in Section 4.(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 Section 5.
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 is the number of Monte Carlo simulations, stands for the estimate of the -th trial for the -th theoretical angle . And, in this paper, we set = 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 is presented in Figure 6, where , , SNR = 5 dB. And, we define the search step and the tiny searching restrain as and . It is shown clearly that the proposed ISF algorithm detects the source signals successfully.

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 is given in Figures 7 and 8 by the same algorithm. It is defined that 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.


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 and It is defined that 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.


5.4. RMSE with Different Snapshots and SNR
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.


5.5. Estimation Probability Comparison of Different Algorithms
Figures 13 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 . The two sources can be resolved if where , [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.


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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by China NSF under Grant nos. 61371169, 61601167, 61601504, and 62004100, Jiangsu NSF under Grant no. BK20161489, the open research fund of National Mobile Communications Research Laboratory, Southeast University under Grant no. 2015D03, Natural Science Research Project of Higher Education in Jiangsu Province under Grant no. 20KJD430010, and Jiangsu Wind Power Engineering Technology Center of Research on wind farm data acquisition system based on 5G and University level scientific research under Grant nos. YK20-02-11 and YK19-02-06.