A novel wideband DOA estimation method based on a fast sparse frame

Funding information Science and Technology Commission of the Central Military Commission, Grant/Award Number: 19H863-01-ZD-***-***-01; National Natural Science Foundation of China, Grant/Award Number: 61771015 Abstract In this study, a novel fast wideband direction of arrival (DOA) estimation algorithm is proposed to reduce the computational complexity. First, a multiple measurement vector (MMV)-based compact structure for a wideband signal is established. Combined with the focus operation, the array manifolds of different frequency bins are transformed into the dictionary of the reference frequency. Then, two efficient novel methods named adaptive step-size-based null space tuning with hard thresholding and feedback (ASNHF) and MMV-ASNHF are proposed to process single measurement vector and MMV problem, respectively. Finally, wideband DOA estimation can be achieved by MMV-ASNHF algorithm. Compared with other algorithms, the proposed algorithm has higher accuracy at a low signal-to-noise ratio and lower computational complexity. Simulation results show that the proposed estimator is effective and feasible.


INTRODUCTION
With the wide application of the 5G signal, the research on wideband signals has become a hot issue in many fields. For wideband direction of arrival (DOA) estimation, there are mainly two types of algorithms, which is based on subspace principle. One type of algorithm is to apply narrowband signal processing methods to each sub-band of the wideband signal, such as incoherent signal subspace method and test of orthogonality of projected subspaces [1]. These algorithms usually have low accuracy. The other type of algorithm is using the focus operation to transform the energy on each sub-band to the same frequency point, which can solve the problem that the energy of wideband signal disperses to each frequency bin [2]. However, there are still some common problems in subspace algorithms. For example, the estimation results are not accurate for limited snapshots and low signal-to-noise ratio (SNR).
In the wake of developments in compressed sensing (CS), the theory of sparse reconstruction has been widely used in the DOA estimation [3,4]. Generally, CS algorithms have high accuracy and good performance with deficient snapshots and low SNR. Although a single measurement vector (SMV) model is available for each sub-band, the joint sparsity of the wideband signal is underutilised. Based on the structure of multiple measurement vectors (MMVs), many researchers have proposed wideband DOA estimation methods using block sparsity approaches [5][6][7]. Although the accuracy of these algorithms is higher than that of the subspace methods, the prohibitive computational complexity restricts their development in practice. Especially when the number of sensors and frequency bins is large, the dictionary of the block sparsity model will become a huge dimension matrix.
Recently, many researchers have utilised statistical learning to solve sparse recovery problems. The most famous one is the sparse Bayesian learning (SBL). To process joint sparse signals, researchers have extended SBL to block SBL (BSBL) [8,9]. Although BSBL algorithms can be directly applied to wideband DOA estimation with high accuracy, the desired result is incapable without multiple iterations [10][11][12].
To overcome the above problems, the focusing operation is used to transform the block sparse dictionary into the single dictionary at the reference frequency bin with lower dimension. The MMV structure for the wideband signal is established and the computational complexity in wideband DOA estimation is greatly reduced. Furthermore, an efficient CS-based algorithm is considered for accelerating operation. The algorithm proposed in [13] is a fast hard thresholding (HT) algorithm (null space tuning (NST) with HT and feedback (FB); (NHF)) to process SMV problems. It is proved that the NHF algorithm is faster than most of the traditional sparse recovery algorithms, and the recovery effect is good. However, it is found that this method has a premature convergence problem in some scenarios by simulation experiments. Consequently, a new algorithm named adaptive step-size based NHF (ASNHF) is proposed. ASNHF is a fast thresholding algorithm with adaptive stepsize and double HT operators. Compared with NHF, the new algorithm may overcome the premature convergence problem and improve the probability of exact reconstruction when the dimension of the observation vector is small. In addition, the MMV-ASNHF is presented to deal with the MMV problem and wideband DOA estimation.
The rest of this study is organised as follows. In Section 2, both the traditional wideband array signal model and a wideband signal model for sparse reconstruction are studied. In Section 3, the NHF algorithm is introduced; then, the ASNHF algorithm and MMV-ASNHF algorithm are proposed. In Section 4, the performance and computational time of different algorithms are compared by simulation.

Classical wideband DOA model
Assuming there are K wideband far-field signals impinging on the uniform linear array (ULA) of M sensors from different directions, and the distance between the elements is d . The received signals of the mth sensor can be expressed as where s k (t ) denotes the kth signal, k denotes the direction of kth signal. m ( k ) is the propagation delay of the kth source at the mth sensor with respect to the reference sensor of the array. n m (t ) is the additive Gaussian white noise of the mth array element. Each received signal is divided into P non-overlapping groups with the same length P L , the time duration of each group is Δt = P L ∕ f s , where f s means the sampling frequency. The propagation delay between the m 1 th sensor and the m 2 th sensor can be defined as m 1 m 2 . Following the widely adopted wideband model assumption in array signal processing, it is assumed that Δt is much larger than the maximum of m 1 m 2 . A succinct definition for wideband signal in array signal processing is that the relative bandwidth should not be less than 0.2.
The frequency-domain model of the wideband signal can be obtained by the discrete Fourier transform of formula (1): which can be written as where N represents the number of frequency bins, f n represents the frequency of the nth frequency bin.

Wideband DOA model for sparse reconstruction
When the SMV algorithms are applied to wideband DOA estimation, solutions can be obtained by averaging the results of the sub-problems at each frequency bin, which will not be discussed here. Taking advantage of the joint sparsity of wideband signals, one can reconstruct the dictionary matrix and observation vector, that is, where B( k )is a block diagonal matrix, which can be written as The observation vector should be reconstructed as where , , X ( f N )] only needs one snapshot in the frequency domain, K g is the number of scanning angles. It can be seen that the dimension of the block dictionary matrix is MN × K g . With the increase of the number of array sensors, frequency bins and scanning angles, the computational complexity will become practically prohibitive. Therefore, a focusing algorithm is used to reduce the dimension of the block dictionary. The goal of the rotational signal subspace (RSS) algorithm [2] is to minimise the error between the array manifold focused by each frequency bin and the array manifold of the reference frequency point, that is, One particular solution is where T ( f n ) is the focusing matrix of frequency f n , f 0 is the reference frequency, U ( f n ) and V ( f n ) are left and right singular vectors of A( f n , )A H ( f 0 , ), respectively. Directly using the above focusing matrix to focus the observation vector, one can get A( f 0 , ) is the dictionary for the wideband DOA estimation with the MMV structure. It can be seen that the dictionary matrix of each frequency bin is transformed into the same dictionary matrix of the reference frequency point. Now, the wideband array signal model is transformed into the MMV model, and the dimension of the dictionary in Equation (4) is greatly reduced.

NHF algorithm
It can be seen from the above analysis that a large amount of calculation in the wideband DOA estimation is mainly caused by two aspects. One is that the dimension of the joint dictionary of each frequency bin is too large. The other is that the complexity of the algorithm is too high. Only the first problem is solved by transforming the block sparse dictionary into the low-dimensional dictionary at the reference frequency bin. After that, if complex algorithms are used, DOA estimation results still cannot be obtained quickly. The NHF algorithm is suitable for simplifying the calculation process. It is proved that the NHF algorithm is faster than most of the traditional sparse recovery algorithms, and the recovery effect is good [13]. By using the NST with HT and FB operations, the 'tail' of the recovery vector with less energy can be concentrated on other non-zero elements, and thus the sparse solution can be quickly obtained.
Sparse recovery can be expressed as the following optimisation problems: Suppose A is a dictionary matrix with row full rank, then NST can be expressed as the following equation NST : where N T denotes some principles to get the approximate solution u t , the corner marker t indicates the variable at the t th iteration. For example, if s t = [s 1 , s 2 ] T , u t can be obtained by setting the first term to zero, that is, u t = [0, s 2 ] T . In this case, u t does not satisfy the equation in Equation (11). Therefore, u t should be transformed by the projection operation P = I − A H (AA H ) −1 A, which represents the orthogonal projection onto the null space of matrix A. The sparsest solution is obtained through multiple iterations. The concept of NHF is similar to this, but in the first step of NHF, a 'tail' is added, which is exactly the FB operator. In a nutshell, the purpose of FB is to feed the energy of the elements with smaller values back to the entries with larger values. The NHF algorithm can be expressed as NHF : The purpose is to recover a K -sparse signal, H T is the HT operator, which keeps the largest K entries unchanged and sets all the others to zero. In this way, the elements with smaller values in the approximate solution s t can be directly set to zero because their energy can be converted into the larger ones through FB operation, that is, where Λ t K is defined as the index set of K larger values in the sparse solution s t , and Λ t K c is the complement of Λ t K . In the actual operation, the two sets can be directly obtained by sorting the elements in s t .
This step can be realised by the least square method, that is, Although the u t obtained in this way could not directly satisfy the equation, it greatly accelerates the convergence speed.

ASNHF algorithm
In the experiment, it is found that the NHF algorithm has some problems and the probability of exact recovery is reduced if the dimension of the observation vector is small. In DOA estimation, the dimension of the observation vector is determined by the number of sensors. A large number of sensors will lead to a huge amount of computation, which is not conducive to engineering implementation. When the number of sensors is small, the accuracy of DOA algorithms based on NHF will decrease. In [14], when using NHF algorithm to narrowband DOA estimation, only one signal is measured at a time, then the received signals are projected onto the null space of the signal that has been measured. Although problems can be avoided in this way, the projection operation of each angle will undoubtedly increase the computational complexity.
To improve the NHF algorithm, the reason for its performance degradation should be discovered. s t +1 can be split as Equation (19), which gives a closed form of s t +1 show an uncertain relationship from Equation (19).
The convergence of the algorithm is proved in [13], but it does not present the iteration in which the algorithm converges. If the NHF algorithm converges in the t th iteration, new elements will not be chosen from the index set in the (t + 1)th iteration, which requires min As long as there is a larger element in s t Λ k c , the index of this element will be selected in the next iteration. A column vector z t ∈ C M is defined as can be easily expressed as ith row of s t +1 and u t +1 , respectively, that is, Since the index of minimum element in s t should be compared, that is, As the index set cannot represent the order of atomic in the dictionary A, it is impossible to determine whether the result of the above Equation (20) is positive or negative. A positive result means that the algorithm converges in the t th iteration. In the simulation, it is found that the inequality in Equation (15) holds in the first iteration when the dimension of the observation vector is small. In other words, the algorithm has the problem of premature convergence when M is small. Once the index selected in the first step is wrong, the correct result cannot be obtained in the last iteration because of the premature convergence.
A novel algorithm ASNHF is proposed to overcome this problem. Assuming that the dimension of a K -sparse vector is K g , in the NHF iterating operation, K largest elements of the sparse vector are selected in each iteration, and then the values of other K g − K smaller elements are transformed to K large elements by least square. In order to increase the robustness of the algorithm and overcome the premature convergence problem, L more elements are selected at each iteration by H T 1 . Then, another HT operator H T 2 is used to delete the previous L entries. The fault tolerance rate of the algorithm can be effectively improved by the above operations. Even if the locations of K larger entries are wrong in any iteration, as long as the correct index is in the increased step size, the correct index can be selected at the end.
An example is shown in Figure 1 to illustrate the difference between the two algorithms. Assuming that the wrong index is obtained in the first iteration of both algorithms, the NHF algorithm converges in the first step when the dimension of the observation vector is small. However, ASNHF algorithm can add and delete indexes in each iteration through step size L and double HT operators H T 1 and H T 2 until the correct result is obtained.
The H T 1 and FB operators in ASNHF algorithm can be expressed as To further improve the performance of ASNHF algorithm, the selection of step size is added. The error t can be set as

FIGURE 1 An example of two algorithms
Suppose l is a small constant, after each iteration, the step size L is updated adaptively by error judgment. If the error increases, L = L + l ; otherwise, the step size does not change. Adaptive step size L + l can overcome the problem of premature convergence more flexibly.
It is equivalent to concentrating the energy of the sparse vector into K + L larger values. Because the K -sparse result is desired, an additional delete operation is needed at the end of each iteration, that is, where Λ t K +1:K +L represents the index set corresponding to the L elements selected in Equation (21), and Λ t (K +1:K +L) C means the other indices. This step can be understood as a HT operator H T 2, which sets the L entries to zero and keeps the others unchanged. So H T 2 and projection operators can be written as The pseudocode is in Table 1.

MMV-ASNHF algorithm
ASNHF algorithm can deal with the DOA estimation of narrowband signals or each sub-band in the wideband signals.
When it is used in the MMV wideband signal model mentioned in Section 2, ASNHF should be extended to the MMV case. In this section, the MMV-ASNHF algorithm is presented. The difference between MMV-ASNHF and ASNHF algorithm is that the selection method of index set is different. Elements in ASNHF are chosen in comparison, while in MMV-ASNHF, the rows are selected. MH T 1 is defined as an operator that keeps the largest K rows of a matrix and sets the elements in other rows equal to zero. To sort rows of the matrix, a norm constraint is imposed on each row to get the largest K rows. The MH T 1 and FB operators in MMV-ASNHF algorithm can be expressed as where ‖Λ t K +L ‖ denotes the index set of K + L largest rows in S t , and ‖Λ t The other steps are the same as ASNHF. MH T 2 is defined as an operator that sets the L rows of S t to zero and keeps all the other rows unchanged, that is, where ‖Λ t K +1:K +L ‖ represents the index set corresponding to the L more rows we selected, and ‖Λ t (K +1:K +L) C ‖ means the other indexes.
One can use Equation (9) as the wideband signal model for MMV-ASNHF. The result of wideband DOA estimation can be calculated by imposing a norm constraint to each row of the returned value S t . The pseudocode of MMV-ASNHF is in Table 2.
It is worth noting that in [12], a method of further simplifying the operation is proposed in which the term is written as

SIMULATION
In this section, simulations results are given to show the effectiveness and feasibility of the proposed algorithm. All the simulation results are obtained using the same computer with an Intel i7-8700 CPU and 32 GB of RAM, running MATLAB R2017a on 64-bit Windows 10. Case 1: In this case, the proposed algorithm is compared with NHF algorithm, Orthogonal Matching Pursuit (OMP) algorithm [15], Subspace Pursuit (SP) algorithm [16], Basis Pursuit (BP) algorithm [17], Generalized Orthogonal Matching Pursuit (gOMP) algorithm [18] and Extended Orthogonal Match- ing Pursuit (eOMP) algorithm [15] in performance and running time. Suppose that the observation vector is a Gaussian random signal, then the zero position is randomly generated, and the dictionary matrix is a random matrix. When the dimension of original signal s is set to 200, the dimension of observation vector x is 80, which means that the   Figures 2(a) and (b). When the dimension of original signal s is set to 100, the dimension of observation vector x is 20, which means that the dictionary A ∈ ℂ 20 × 100 , and the simulated results are shown in Figures 2(c) and (d). Assume that the   Figure 2(a) shows the curve of the operation time of various algorithms with the increase of the sparsity, and Figure 2(c) shows the curve of the probability of exact recovery of various algorithms with the increase of sparsity. Figures 2(b) and (d) show the change of the step size L with the number of iterations, and it updates automatically according to the error t . As shown in Figure 2(a), NHF algorithm has the least operation time, followed by the proposed algorithm. This is because the proposed algorithm has double HT operators and selects more entries in each iteration. Because the running time of BP algorithm is too long, it is not plotted in the figure. In Figure 2(c), it is found that the performance of the NHF algorithm is not ideal when the number of rows in the dictionary matrix is much smaller compared with the number of columns. In this case, the performance of the proposed algorithm is better than that of NHF algorithm, which is close to other sparse recovery algorithms; moreover, the running time is less than other algorithms.
Case 2: In this case, the performance of the NHF algorithm is compared with that of ASNHF algorithm in narrowband DOA estimation. Parameters used in the simulation are given in Table 3. The element spacing equals the half wavelength at the centre frequency. The ULA is assumed to be ideal without the impact of array imperfection and mutual coupling. The simulation results are shown in Figure 3. It can be seen that when NHF algorithm is used for narrowband DOA estimation, only one target can be detected accurately without projection operation. This is the problem mentioned at the beginning of Section 3. ASNHF algorithm can accurately detect multiple targets.
Then, the MMV-ASNHF algorithm is applied to wideband DOA estimation, and signal model is defined as Equation (9). As mentioned in Section 1, wideband DOA estimator can be established by subspace methods or block sparsity-based algorithm, including BSBL algorithms. Therefore, three algorithms are studied for comparison: BOMP algorithm [7], BSBL  [9] and the traditional RSS algorithm [2], which can represent most algorithms for wideband DOA estimation in the frequency domain. Furthermore, two wideband DOA estimators presented most recently [19,20] are taken for comparison.
The two sources at (−2 • , 1 • ) are coherent and the last source is uncorrelated to them. The simulation parameters are shown in Table 4. The number of snapshots in the frequency domain is 1. [−10 • : 10 • ] is the prediction angle range of focusing operation, and the step is set as 0.1 • within (−90 • , 90 • ). The array spacing is the half wavelength at the centre frequency. The simulation results are shown in Figure 4. It can be seen that the performance of RSS is poor with limited snapshots, as the two signals are coherent. The other three algorithms can measure three signals.
Case3: To further illustrate the performance of these algorithms, the resolution probabilities of the algorithms are studied for comparing. Suppose there are two coherent wideband signals impinging from (0 • , 2 • ). The number of elements and subbands are 40 and 11, respectively. The other simulation parameters are the same as that in Table 4. The resolution result with respect to SNR is shown in Figure 5(a). Then, the SNR is set as −10 dB, and results versus the number of elements and subbands are displayed in Figures 5(b) and (c), respectively, when the source separation is varied from 1 to 6 • . Figure 5(d) depicts the resolution probability with respect to angle separation. The two signals are resolved if both |̂1 − 1 | and |̂2 − 2 | are smaller than 1 • . Monte Carlo trial number T is 200. As shown in the results, our algorithm is only a litter weaker than BSBL algorithm but better than the other algorithms in resolution probability.
Case 4: In this case, the Root Mean Squard Error (RMSE) of the six algorithms is compared for wideband DOA estimation, which is mentioned in Case 2. The RMSE can be defined as Monte Carlo trial number T is 200, K is the signal number. Suppose that two signals are incident from (−10 • , 10 • ), then the SNR increases from −5 to 10 dB gradually, and the prediction angle range of RSS algorithm is (−15 • , 15 • ). The number of sensors is 20 and the other simulation parameters are the same as that in Table 4; CRB is the wideband Cramer-Rao lower bound [21] and the simulation results are shown in Figure 6(a). Next, the number of frequency bins is changed from 11 to 20, SNR = −5 dB, and the results are shown in Figure 6(b). Then, the number of elements is changed from 16 to 30, the number of sub-bands is 11, and the results are displayed in Figure 6(c),  The impact of multipath propagation mode is also considered, and two coherent signals are used in this case. Simulation results are shown in Figures 6(e) and (f). The RSS algorithm has poor performance with limited snapshot, and the two algorithms proposed in [19,20] are inapplicable to coherent signals; thus, they are not plotted in Figures 6(e) and (f). It can be seen that the BSBL algorithm has the highest estimation accuracy but with the longest running time. The estimation accuracy of MMV-ASNHF algorithm is better than the other algorithms.
Case 5: In this case, the computational complexity of the above four algorithms is compared. The main computational complexity of the RSS algorithm includes calculation of SVD, EVD and spectrum peak searching, which requires The main computational complexity of the BOMP algorithm in each iteration includes calculation of the inner product, the leastsquares problem and the residual, which requires K g MN 2 + Block sparse Bayesian learning (BSBL)-FM (one iteration) algorithm in each iteration includes calculation of FB and P, which requires (M + N )(K + L) 2 + MN (K g + L) + K 2 g N flops. P can be calculated and stored in the memory before the first iteration. Table 5 shows the computational complexity comparison of the above algorithms.
Then, their computational time is compared. The influence of matrix dimension on algorithms of block sparse is reflected in the number of blocks, the number of elements in each block, and the dimension of observation vector. These three items  correspond to the number of scan points, frequency bins and array elements in wideband DOA estimation, respectively. So four variables are set in simulations. First, let the SNR increase from −5 to 10 dB with step 1 dB, while other conditions remain unchanged. The simulation results are shown in Table 6. Next, the number of scanning points is changed, SNR = 5 dB, and the frequency bins are 11, other conditions remain unchanged. The results are shown in Table 7. Then, the number of frequency bins is changed from 4 to 64, other conditions are unchanged. The results are shown in Table 8. Finally, the number of array elements is changed from 5 to 50 with step 10, and the results are shown in Table 9. From the simulation results, one can see that the operation time of the proposed algorithm is much shorter than that of the other three algorithms. Although the estimation accuracy of BSBL algorithm is very high, the computational complexity is too large, which is not suitable for practical application.

CONCLUSION
In this study, a fast wideband DOA estimator is presented. Specifically, a novel fast HT sparse recovery algorithm ASNHF is proposed, and then it is extended to MMV-ASNHF for processing MMV problem. Finally, MMV-ASNHF is applied to wideband DOA estimation. Simulation results show that the greatest advantage of the new algorithm is the incredible efficiency. In addition, MMV-ASNHF has an excellent performance in accuracy and resolution. This work contributes to extending the application of CS theory in wideband DOA estimation.