Fast DOA Estimation Based on the Transform Domain Weighted Noise Subspace Fitting Algorithm for Generalized Sparse Array

Compared with uniform arrays, a generalized sparse array (GSA) can obtain larger array aperture because of its larger element spacing, which improves the accuracy of DOA estimation. At present, most DOA estimation algorithms are only suitable for the uniform arrays, while a few DOA estimate algorithms that can be applied to the GSA are unsatisfactory in terms of computational speed and accuracy. To compensate this deﬁciency, an improved DOA estimation algorithm which can be applied to the GSA is proposed in this paper. First, the received signal model of the GSA is established. Then, a fast DOA estimation method is derived by combining the weighted noise subspace algorithm (WNSF) with the concept of “transform domain” (TD). Theoretical analysis and simulation results show that compared with the traditional multiple signal classiﬁcation (MUSIC) algorithm and the traditional WNSF algorithm, the proposed algorithm has higher accuracy and lower computational complexity.


Introduction
With the rapid development of information technology, wireless transmission systems, such as mobile communications, radar, and unmanned aerial vehicles (UAV), have become widely used. As such, antenna technology has become more important [1][2][3][4][5], as it has array antenna direction finding, which is an important branch of antenna technology. Most of the existing arrays are uniform arrays, whose element spacing is less than half wavelength of the incident signal, so as to avoid the ambiguity of angle estimation [6][7][8][9]. However, when the number of array elements is limited, the aperture of the uniform array will also be affected, which will lead to lower accuracy of DOA estimation [10,11]. To solve this problem, scholars have proposed the GSA-a nonequidistant array system-in which the distance between adjacent array elements is more than half of the wavelength of the incident signal [12][13][14][15]. Compared with the uniform arrays, the GSA can obtain larger array aperture. erefore, the estimation accuracy and resolution of the algorithm are effectively improved. e GSA is derived from a uniform array, and it can be obtained by redeploying the array elements of uniform array according to some optimization algorithms (such as simulated annealing algorithm [16,17] and genetic algorithm [18,19]). For instance, Chen et al. [20] has proposed an effective method to construct GSA based on the improved genetic algorithm. e most significant feature of GSA is that the spacing between adjacent array elements is unequal. In this case, most conventional DOA estimation algorithms which have strict requirements on array structure, such as ESPRIT algorithm [21][22][23], will lose their effect. On the other hand, a few algorithms that are less dependent on the array structure, such as the MUSIC algorithm [24][25][26], do not perform well in the condition of low SNR. At present, researchers generally use the WNSF algorithm for DOA estimation [27].
e WNSF algorithm has no requirement on array structure and can accurately estimate signal DOA under bad conditions such as low SNR or small number of snapshots [28][29][30]. However, the WNSF algorithm needs spectral peak search. In accurate search or two-dimensional DOA estimation, the total number of points of the spatial spectrum is very huge, which incurs vast amount of computation and hinders the application of the algorithm in practical engineering [31,32].
To reduce the computation, a reduced-dimension search algorithm was proposed in [33,34], which could reduce the dimensionality of the spectral function and then performed multiple low-dimensional searches. However, this algorithm not only reduced the amount of computation but also reduced the accuracy of parameter estimation. To solve this problem, a concept of transform domain (TD) was proposed [35,36]. e authors in [35] transformed the received signal into the coarray domain and then iteratively corrected the phase offset between the coarray data and presumed model caused by angle biases according to a closed-form formula. On the other hand, Zhang et al. [36] used spherical Fourier domain to construct the array signal model. However, these methods are only applicable to uniform liner arrays, which means that the method is not universal. Nevertheless, it is a concept worth learning and developing.
Motivated by the above facts, in this paper, we combine the concept of "transform domain" with the WNSF algorithm to obtain a fast DOA estimation method which is suitable for the GSA. e procedure of the proposed method is as follows. Firstly, we establish the received data model of the GSA. Secondly, the model is transformed into the TD. en, the intersection of the noise subspace and its conjugate space in the TD is used to replace the noise subspace in the spectral search function. In this manner, the calculation range of the search function can be reduced by half, and the "rapidity" of the algorithm can be achieved. e theoretical analysis and computer simulation results show that the proposed algorithm has higher accuracy and lower computational complexity than the MUSIC algorithm and traditional WNSF algorithm.

Signal Model of GSA
Consider a GSA composed of M elements, as shown in Figure 1. e array elements are arbitrarily distributed in space, and the position coordinates of the array elements are given by (x m , y m , z m ), m � 1, 2, . . . , M. e black dots in Figure 1 represent the array elements. e elevation angle θ is defined as the intersection angle between the signal and the positive direction of Z-axis, while the azimuth angle φ is defined as the angle between the projection of the signal onto the XOY plane and the positive direction of the X-axis. Here, the value ranges of θ and φ are both (− π/2, π/2). Assume that K narrow-band waves impinge upon the GSA from the elevation angle θ k and the azimuth angle φ k , the data received by a snapshot of the array can be expressed as where , s 2 (t), · · · , s K (t)] T is the signal sampling data; and n(t) is the additive white Gaussian noise matrix with the same dimension as x(t). e array steering matrix A(θ, φ) ≜ [a(θ 1 , φ 1 ), a(θ 2 , φ 2 ), . . . , a(θ K , φ K )] and the steering vector a(θ k , φ k ) can be expressed as a θ k , φ k ≜ e jβ k,1 , e jβ k,2 , . . . , e jβ k,M , where j � �� � − 1 √ is the plural unit. Define λ as the signal wavelength. en, β k,m is given by e process of estimating the 2D spatial angle of the K signals based on the signal model constructed by (1) is called spatial spectral estimation.

DOA Estimation Algorithm Based on TD-WNSF
3.1. Principle of the Traditional WNSF Algorithm. In practical application, according to the received data from P snapshots, the spatial correlation is estimated using a time average, and the covariance matrix of the array output is obtained by where X is the received data matrix of P snapshots. After performing an eigenvalue decomposition, we can derive the following result: where U S is the signal subspace composed of eigenvectors corresponding to K larger eigenvalues and U N is the noise subspace composed of eigenvectors corresponding to M − K smaller eigenvalues. Accordingly, we can deduce that the signal subspace and the space formed by the steering vector of the array are the same, and the steering vector space of the array and the noise subspace are orthogonal to each other [37]. is orthogonal relation can be expressed as where the symbol 'H' represents the conjugate transpose operation, a(θ, φ) is the M × 1 steering vector, and O is the Considering that the length of the actual received data matrix is limited and noise is mixed in the actual received data matrix, a(θ, φ) and U N are not completely orthogonal. In other words, equation (6) is not completely valid. erefore, the following formula of noise subspace fitting (NSF) is considered [38]: Equation (7) can be further extended to a weighted form. e relationship between the noise subspace and the steering vector space of the array is given by where W is the weight. en, the fitting formula given by equation (7) can be transformed into eoretically, the DOA estimation in equation (9) can be estimated by bringing a 2D search to bear on the ranges of all parameters; however, this is computationally exhaustive.

TD-WNSF Algorithm.
When running the traditional WNSF algorithm, an extreme value test should be conducted for each point in the spatial spectrum. With the improvement in search accuracy, the number of spectral points is increased further, leading to a sharp increase in the running time of the algorithm. If we can find a way to compress the extremum search range, then the speed of DOA estimation can also be improved. From this analysis, the following transformation is considered: Combining equations (3) and (10) leads to Moreover, combining equations (2) and (11) leads to where * represents the conjugate operation. e steering vector a(θ, ϕ) and the transformation domain steering vector a(u, v) are clearly equivalent in a physical sense. en, (12) can be depicted as a virtual mirror signal source symmetry to (u, v). Assume that the weight is given by W � I, where I represents the identity matrix. Moreover, combining equations (6) and (12) leads to In equation (13), the steering vector corresponding to the real signal source (u, v) is orthogonal to the noise subspace U N , whereas the steering vector corresponding to the virtual mirror signal source (− u, − v) is orthogonal to the conjugate noise subspace U * N . If we replace the noise space in the WNSF algorithm with the intersection space of U N and U * N , as the intersection space is orthogonal to the real steering vector and the virtual steering vector, then the WNSF algorithm can generate an extremum at the real signal source and virtual signal source simultaneously. is characteristic means that the DOA estimation only needs to search half of the (u, v) domain. erefore, the purpose of "rapidity" can be achieved. is algorithm of constructing the WNSF spectral function in the TD is called the TD-WNSF algorithm.
Solving the intersection space of U N and U * N is essential to constructing the TD-WNSF algorithm. e steps discussed briefly describe the method to find the intersection space of the two subspaces.
Suppose that the noise subspace U N is represented by [α 1 , α 2 , . . . , α M− K ] and the conjugate noise subspace U * N is represented by [β 1 , β 2 en, the intersection space U inter of U N and U * N can be determined by the adjoint solution of which is specifically expressed by Definition 1 in this study.
. . , k s )} is the adjoint solution of the solution of equation (14).
By contrast, if (k 1 , k 2 , . . . , k s ) is the adjoint solution of a solution to that of equation (14), then (k 1 , k 2 , . . . , k s , l 1 , l 2 , . . . , l s ) is the solution of equation (14), which is also the expression of equation (17). Given that the left side of equation (17) belongs to U N , whereas the right side of equation (17) In summary, U inter � U N ∩ U * N � k 1 α 1 + k 2 α 2 + · · · + k s α s |(k 1 , k 2 , · · · , k s )} is the adjoint solution of the solution of equation (14). At this stage, the proof has been completed.
On the basis of on the above algorithm, the TD-WNSF spectrum can be defined as As mentioned, this TD-WNSF spectral function generates extreme values at the real and mirror positions of the signal at the same time. us, only a half-domain search in the (u, v) domain is needed to estimate the signal direction.
en, an accurate search is performed at the vicinity of the spectral peak position and its mirror image position to obtain the precise arrival angle. On the basis of the previous analysis, the proposed algorithm does not entail any array structure requirement and can be applied to the GSA.

Description of Algorithm
Step. e implementation steps of the proposed method are summarized as follows: (i) Step 1: perform eigenvalue decomposition of the array-received data in order to obtain the noise subspace. (ii) Step 2: calculate the intersection space based on the algorithm of Definition 1 and then construct the TD-WNSF spectrum based on equation (19). (iii) Step 3: search the positive half-spectrum of equation (18) to obtain the estimated value (u i , v i ) of the DOA parameter in the (u, v) domain, where i � 1, 2, . . . , K. (iv) Step 4: replace a(θ, φ) in equation (6) with a(u i , v i ) and for the extreme value test. Among the two elements, the one that satisfies a H (u, v)U N � 0 is the real TD-DOA. (v) Step 5: substitute the TD-DOA, which was obtained in step 4, into equation (10) to calculate the spatial DOA (θ i , φ i ). (vi) Step 6: perform an accurate searching in a small area near (θ i , φ i ). e element that satisfies a H (θ, φ)U N � 0 is the real spatial domain DOA.
As shown in these steps, the algorithm first implements a rough search process to obtain the TD-DOA (u i , v i ). en, the inverse trigonometric function is transformed to calculate the rough estimation of the angle (θ i , φ i ). Finally, an accurate search is conducted in the small neighborhood of (θ i , φ i ) to obtain the fine estimate (θ ⌢ i , φ ⌢ i ). As such, no angle measurement blurring will occur.

Algorithm Complexity Analysis.
e traditional MUSIC algorithm, the traditional WNSF algorithm, and the TD-WNSF algorithm are compared in this study. e array structure used by the three algorithms is shown in Figure 1. Consider that K uncorrelated signals impinge upon an GSA of M elements, the number of snapshots is given by L, and the number of search points is given by q. Assume that all the three algorithms perform a rough search, which means that the step size is 1°.
For the traditional MUSIC algorithm, the modulus ‖a H (θ, φ)U N ‖ 2 needs to be calculated for each spectral point, and the dimension of U N is M × (M − K). erefore, the computation of the spectral search of the traditional MUSIC algorithm is q(M − K)(M + 1). e computation of the eigenvalue decomposition in the M × M dimensional autocovariance matrix is M(K + 2) 2 . erefore, the total computation of the traditional MUSIC algorithm is For the traditional WNSF algorithm, the weight W is usually taken as a unit matrix I in the calculation. e operation used to find the trace tr(‖a H (θ, φ)U N ‖ 2 ) is needed to calculate each spectral point. e dimension of U N is M × (M − K), and the computation involved in finding the trace is q(M − K)(M + M). e computation of the matrix's eigenvalue decomposition is the same as that of the traditional MUSIC algorithm, which is M(K + 2) 2 . erefore, the total computation of the traditional MUSIC algorithm is For the TD-WNSF algorithm proposed in this study, the operation used to find the trace tr(‖a H (θ, φ)U inter ‖ 2 ) is needed to calculate each spectral value point. e intersection space U inter has a lower dimension denoted by M × (M − 2K) compared with U N , and the search range of the TD-WNSF algorithm is reduced by half.
us, the computation of the peak search is q(M − 2K)(M + 1)/2, while the computation of the eigenvalue decomposition for the autocovariance matrix in the TD-MUSIC is erefore, the total computation of the TD-MUSIC algorithm is q(M − 2K)(M + 1)/ 2 + 2(M − 1)(K − 1) 2 . Table 1 shows a comparison of the computations of the three algorithms with the number of array elements.
We can easily see that, with the increase of the number of array elements, the computational complexity of our algorithm is far lower than the other two algorithms, which reflects the rapidity of our algorithm.

Simulation and Analysis
Consider a generalized sparse array with 16 array elements, the position distribution of its array elements is shown in Table 2.

Simulation 1: Signal DOA Parameter Estimation.
is simulation was implemented to evaluate the effectiveness of the proposed TD-WNSF algorithm. Considering a signal with the following values, θ � 60 o and φ � 20 o , impinging upon the GSA. SNR is taken to be 0 dB. e running result of our algorithm is shown in Figure 2. From Figure 2, we can see that the estimation of TD-DOA is (0.17, 0.865). According to equation (10), the spatial angle (60°, 20°) and the TD angle (0.17, 0.865) are equivalent, which illustrates that our algorithm can accurately measure the DOA of the signal.

Simulation 2: Relationship between Algorithm Performance and SNR.
is simulation is implemented to compare the DOA estimation performances of the MUSIC, WNSF, and TD-WNSF algorithms.
e Monte Carlo number is L � 500, and the SNR shifts from − 10 to 20 dB. e rootmean-squared error (RMSE) is defined as where θ k,l and φ k,l represent the estimated values of the Kth signal in the Lth Monte Carlo simulation, respectively. Figure 3 shows the relationship between the RMSE of the azimuth angle and SNR, while Figure 4 shows the relationship between the RMSE of the elevation angle and SNR. From Figures 3 and 4, we can see that the RMSE curves of three methods decrease as SNR increases. However, the RMSE curve of our TD-WNSF algorithm locates below that of MUSIC and WNSF all along, which indicates that the TD-WNSF is superior to MUSIC and RD-MUSIC in a different SNR.
is is because our method performs many finegrained searches in the neighborhood of the TD estimation.    International Journal of Antennas and Propagation

Simulation 3: Relationship between Algorithm Performance and the Angle of Incident Signal.
To observe the performance of the new approach more clearly, we compared the algorithm performance with the angle of incident signal. e simulation consists of two parts. e first part entails fixing the elevation angles to 15°, 30°, 45°, 60°, and 75°. e change in azimuth angle is conducted to explore the relationship between RMSE and azimuth angle. e second part involves fixing the azimuth angles to 15°, 30°, 45°, 60°, and 75°. e change in elevation angle is performed to explore the relationship between RMSE and elevation angle. e simulation result of the first part is shown in Figure 5. Figure 5 shows the relationship between the performance of algorithms and the azimuth angle of the incident signal under five typical elevation angles. e comparison of the five subgraphs indicates that the array has high sensitivity to the signals with azimuth angle between 0°and 40°. e simulation result of the second part is shown in Figure 6. Figure 6 shows the relationship between the performance of algorithms and the incident elevation angle of the incident signal under five typical azimuth angles. e comparison of the five subgraphs indicates that the array has high sensitivity to the signals with an elevation angle between 10°and 35°.

Conclusions
In this paper, a novel and efficient TD-WNSF algorithm for estimating the DOA was proposed, which can be applied to the GSA. We first established the received data model of the GSA. And then, we exploited the symmetry of the fitting estimator in the TD to reduce the spectral searching range. Compared with the traditional MUSIC and WNSF algorithms, the computational complexity of the proposed algorithm is significantly reduced. e simulation results indicate that the TD-WNSF algorithm has high accuracy and efficiency.
Data Availability e data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.