A novel aliasing-free subband information fusion approach for wideband sparse spectral estimation

Wideband sparse spectral estimation is generally formulated as a multi-dictionary/multi-measurement (MD/MM) problem which can be solved by using group sparsity techniques. In this paper, the MD/MM problem is reformulated as a single sparse indicative vector (SIV) recovery problem at the cost of introducing an additional system error. Thus, the number of unknowns is reduced greatly. We show that the system error can be neglected under certain conditions. We then present a new subband information fusion (SIF) method to estimate the SIV by jointly utilizing all the frequency bins. With orthogonal matching pursuit (OMP) leveraging the binary property of SIV’s components, we develop a SIF-OMP algorithm to reconstruct the SIV. The numerical simulations demonstrate the performance of the proposed method.


Introduction
Wideband direction-of-arrival (DOA) estimation has been a popular area of research due to the various applications in radar, sonar, seismology, communications, astrophysics, and many other fields [1][2][3]. Traditional wideband array processing is to decompose the wideband signals into many narrowband signals with a filter bank or the discrete Fourier transform (DFT), and two categories, referred to as incoherent signal subspace method (ISSM) [4] and coherent signal subspace method (CSSM) [5], are utilized to realize wideband DOA estimation. The ISSM estimates the DOAs independently and average them over all the bins. The performance of ISSM may deteriorate with low signal-to-noise ratio (SNR) frequency bins and coherent sources. The CSSM align the signal subspaces by transforming the observation vectors associated with each bin into the focusing subspace and can deal with coherent sources by averaging the subspace-aligned covariance matrices. Compared with ISSM, CSSM can enhance DOA resolution and improve the accuracy of DOA estimates at low SNR. However, CSSM requires an initial DOA *Correspondence: luojian@hdu.edu.cn 1 Key Lab for IOT and Information Fusion Technology of Zhejiang, Hangzhou Dianzi University, Hangzhou 310018, China Full list of author information is available at the end of the article estimates, and the precision of DOA pre-estimates greatly influences the accuracy of DOA estimation [6,7].
Recently, a class of sparse signal representation (SSR) methods provide a new perspective for wideband DOA estimation [8][9][10][11]. The DOA estimation problem can be formulated as recovering a spatial sparse signal vector or matrix by minimizing the residual norm under sparsity constraint. Basically, the wideband SSR methods in frequency domain decompose the received signals into narrowband subbands and estimate spatial spectra in each frequency bin. One of the most successive 1 -norm-based SSR algorithms for DOA estimation is 1 -SVD (Singular Value Decomposition) [12], which reduces the computational complexity by SVD. Hyder and Mahata [13] presented a joint 2,0 -norm approximation (JLZA) method and extended it to wideband DOA estimation. Tang et al. [14] showed that the spatial ambiguity can be removed by using multiple dictionaries, each dictionary corresponding to a judiciously chosen frequency. It should be mentioned that the above SSR methods in frequency domain are generally formulated as a multi-dictionary/multi-measurement (MD/MM) joint optimization problem. These techniques usually do not use all the subband information to estimate DOAs with the aim of reducing unknown variables. Similar to ISSM, the performance of joint optimizing a MD/MM problem may deteriorate rapidly if subband signals with low SNR are chosen. Liu et al. [15] proposed a wideband covariance matrix sparse representation (WCMSR) method for DOA estimation. The WCMSR method uses time domain measurements and has its limitation for spatial nonambiguity because the spatial aliasing is frequency dependent.
In this paper, we reformulate the MD/MM problem as a sparse indicative vector (SIV) recovery problem and propose a new subband information fusion (SIF) method to estimate the SIV by jointly integrating all the frequency bins together. Thus, the number of unknown variables can be reduced dramatically. By introducing the SIV, an extra system error is also generated. We show that such error term can be ignored under certain conditions. Compared with the traditional wideband DOA estimation methods, the proposed approach does not rely on focusing technique to average the subspace matrices. By using the binary property of the SIV, we design a binary constrained orthogonal matching pursuit (OMP) algorithm to recover the SIV. The developed algorithm is called SIF-OMP algorithm.
The organization of this paper is as follows. In Section 2, we review the SSR model for wideband DOA estimation. The new SIF method is presented in Section 3. The SIF-OMP algorithm is implemented in Section 4, together with the convergence properties of the SIF-OMP algorithm. Numerical simulations are carried out in Section 5 to demonstrate the performance of our algorithm, and Section 6 gets the conclusion.

Problem formulation and existing methods
Consider a uniform linear array (ULA) of N omnidirectional sensors working together to estimate the spatial location parameters of Q wideband sources and Q is assumed to be unknown. The sensors of the ULA are equally placed on a line with spacing d which is not necessary to be smaller than half a wavelength.
denote the set of a sampling grid of all possible source locations, L N, L Q. We assume that the grid is fine enough that can represent the true source locations, e.g., {θ i 1 , θ i 2 , . . . , θ i Q } ∈ . The subscript i q , q = 1, 2, . . . , Q, is used to index the position of θ i q .
For each sensor, the time-samples are split into M segments, where for each segment, K frequency bins are obtained by a bank of narrowband filters or the discrete Fourier transform [7]. Let s k,q (m) denote the qth source signal at frequency ω k computed for the mth segment, k = 1, 2, . . . , K and m = 1, 2, . . . , M. The sparse representation perspective transforms the DOA estimation problem into sparse spectrum recovery problem. We use v k,i (m) to denote the kth frequency coefficient for the mth segment corresponding to the i-th grid, i = 1, 2, . . . , L. Similarly, we use y n,k (m) to represent the measurement data at the nth sensor for the mth segment at frequency ω k . By stacking all measurements into a vector, the output of the array can be expressed as where y k,m =[ y 1,k (m), y 2,k (m), . . . , y N,k (m)] T is the measurement vector, w k,m represents the N × 1 additive noise vector, a k (θ i ) is the steering vector with respect to the i-th grid and it can be written as where c is the speed of the signal propagation. We now introduce an overcomplete basis matrix . The sparse representation model in (1) can be expressed concisely where where V k is an L × M matrix, W k is an N × M matrix, and they are defined similarly as Y k . If the sources are assumed to be stationary during M snapshots, then each column of V k shares the same sparsity. As shown in the previous literature (see [12][13][14]16] and the references therein), the DOA estimation problem can be solved by reconstructing the operation to form a block diagonal matrix. Thus, the sparse representation for the MD/MM problem can be described as where the 0 norm V 0 denotes the number of non-zero rows of a matrix V and ε is an upper bound of the Frobenius norm of the residual error. However, finding such a combinatorial problem requires an enumerative search and is NP hard. Consequently, V 0 is replaced by V p,t , for 0 < p ≤ 1, t ≥ 1, given in [16] wherev i ∈ C 1×M is the ith row vector of matrix V and v i,m is the m-th element in the row vectorv i . Taking the 1,2 mixed norm minimization as an example, the matrix V k can be estimated by the following constrained optimization problem [12]: The above optimization problem can be solved with standard optimization software, i.e., second-order cone programming (SOCP). However, the optimization procedure should be repeated by K times by utilizing all frequency bins ω 1 , ω 2 , . . . , ω K , and K could be rather large in wideband signal processing. Indeed, it is not efficient to compute V using (7) for all frequency components. In this paper, instead of estimating matrix V , we solve the DOA estimation problem efficiently by recovering a SIV which is used to represent the location of sources.

Subband information fusion method
From the sparse representation model in (4), we observe that the matrices V 1 The nonzero rows of matrix V indicate the source locations. To solve the DOA estimation problem, we need to decide which row of the source matrix V is non-zero from the measurements. Once this is done, say that the ith row of V is non-zero, we can infer that θ i is one of DOA estimates for the corresponding source. Indeed, it is not necessary to estimate the whole signal matrix V to get the solution of a direction estimation problem. Actually, the DOA estimation problem can be formulated as recovering a SIV using a SIF method and the estimating of entire matrix V can be avoided.

Method
To develop the SIF algorithm, first, we introduce a projection matrix C k,i for the ith spatial grid at frequency ω k given by Let ψ k,i (m) = C k,i y k,m denote the projection of y k,m onto the range of a k (θ i ). k,m denotes an overcomplete dictionary whose ith column is ψ k,i (m), and it has the following expression For noiseless case, we get y k,m = Q q=1 C k,i q a k (θ i q ) s k,q (m). Inspired by this relationship, we introduce a new L × 1 binary SIV g = g 1 , g 2 , . . . , g L T whose nonzero elements indicate the location of sources, i.e., g has the same sparsity structure of source matrix V . As such, the model (4) can be rewritten as where E k = e k,1 , e k,2 , . . . , e k,M is an N × M error matrix and the mth column of E k is given by In (10), we observe that g captures the sparsity property of V in the dictionary k,m .
Using (10), we can combine all subband information to estimate a single SIV. Let X k (g) =[ k,1 g, k,2 g, . . . , k,M g] denote the sparse representation matrix of Y k . When all subbands are considered, the new observation model is given as follows where are defined similarly as V . Based on the new observation model (12), the following constrained optimization problem can be used for wideband DOA estimation, where η is the upper bound of the Frobenius norm of the residual error with respect to (12). Note that (13) generally subjects to great error bound compared with (5) for each frequency since additional error E is introduced in the measurement model (12). Note that g is the only unknown parameter appeared in (13). Thus, all the subband information can be utilized jointly to estimate g. We call it subband information fusion (SIF) method. The SIF method in (13) attempts to find a sufficiently sparse g in a manner that X(g) consistently fits Y as sparsely as possible. Instead of estimating V 1 , V 2 , . . . , V K by the MD/MM method stated in (5), the SIF approach only estimates a SIV g and therefore significantly reduces the number of unknown variables during the estimation process because the number of frequency bins K could be rather large in real applications.

Analysis
In Section 3, the system error E appeared in (12) is discarded when we solve the optimization problem (13). This leads to loss of information since E contains the DOA information. This section mainly analyzes the error term E and addresses the problem when E can be neglected.

Single source case
It is straightforward to rewrite the measurement y k,m when only a single DOA is required to be estimated, In consideration of spatial sparsity, y k,m can also be written as Combining all these measurements together, we obtain From (16), we observe that the system error E does not exist for single source case. Therefore, there is no information loss for single DOA estimation when we formulate the MD/MM problem as a SIV recovery issue.

Multi-source case
Let μ denote the parameter controlling the value of e k,m and μ is given by where κ(θ) = sin(θ i q ) − sin(θ i j ). Thus, e k,m can be written as Note from (17) and (18) that the value of e k,m is related to μ j,q . We need to thoroughly discuss the property of μ j,q . According to (17), the magnitude of μ j,q is given by When θ i q approaches close to θ i j infinitely, we obtain lim θ iq →θ i j κ(θ) = 0 and lim θ iq →θ i j |μ j,q | = 1 (20) The above limit conditions show that the error term e k,m can not be neglected if the angles of any two incident sources are close enough. However, e k,m can be very small if κ(θ) is larger than a certain value which can be determined by an upper bound of |μ j,q |, namely B |μ| . B |μ| is given by

SIF algorithm
The SIF method proposed in the previous section suggests that we can solve the wideband DOA estimation problem via single measurement vector reconstruction, as long as the system error E can be neglected. In this section, we first derive an equivalent form for SIF optimization; then, we develop a SIF algorithm to recover the SIV efficiently.
We now develop an equivalent form of (13). Let D ∈ C L×L represent a new dictionary which satisfies D H D = , where One realization of D is given by the eigenvalue decomposition (EVD) of the Hermitian matrix . Assume that can be written as = U U H , where U is a unitary matrix whose columns are composed by L orthonormal eigenvectors and is a diagonal matrix of eigenvalues. Thus, the dictionary D has the expression D = U 1 2 . Let t denote the rank of . We have rank(D) = rank( ) = t. We now define a new measurement vector ξ ∈ C L×1 given such that where h = we can conclude that the solution of (24) exists and ξ = (DD H ) −1 Dh. Based on the above definitions, the equivalent expression of (13) is presented as follows: Since the SIV g has the same sparsity structure of V , the following theorem shows that the problem of estimating g by optimizing (26) is equivalent to estimating g by using (13). (13) and (26), the problem of recovering g through (13) is equivalent to estimating g by (26).

Proof See Appendix.
Theorem 1 suggests that the wideband DOA estimation problem can be solved by optimizing (26) using existing methods, e.g., OMP. We consider the OMP for the recovery of g based on the measurement vector ξ . OMP is an iterative greedy algorithm. For each step, it finds out the column of D which is the most correlated with the Fig. 2 The amplitudes of |μ j,q | and B |μ| with k = 50. The amplitudes of |μ j,q | and B |μ| , ω 50 = 5.9π × 10 7 rad/s, d = πc/(ω 0 + Bf /2), Bf = 0.2ω 0 , ω 0 = 7π × 10 7 rad/s, c = 3 × 10 8 m/s, N = 7 current residuals. This column is then added into the set of selected columns. After the OMP algorithm returns the column set, one can use the least squares (LS) method to further estimate the nonzero values of g. Note that the elements of g are binary by definition. Such property should be incorporated into the OMP algorithm. Once the column set is obtained, the coordinates of "1" components are also determined. Therefore, we do not need to deal with the LS process before updating the residuals.
The proposed SIF method for recovering g can be implemented in terms of algorithm 1, referred as the SIF-OMP algorithm. Let us use q to denote the iteration number. For any subset q ∈ {1, 2, . . . , L}, we denote by D q a submatrix of D consisting of the columns d γ q within γ q ∈ q . We use r q to denote the residual at the qth iteration. Based on the above notations, the SIF-OMP algorithm is given as follows.

Simulation results
In this section, we illustrate numerically the performance of the proposed algorithm via various simulations. We consider a similar example given in [15] for ease of comparison. Assume that two BPSK signals with central frequency of 70 MHz and bandwidth of 20% impinge on a ULA with 7 sensors. The code-rate of incident BPSK signals is assumed unknown. First, we show the DOA estimation results of two sources obtained by the MUSIC, MD/MM [14], WCMSR [15], and SIF-OMP, respectively. Within all approaches, the spatial range [ −90 • , 90 • ] is split into 180 grids with an interval of 1 • . The total number of frequency bins is K = 256. The number of snapshots or segments is M = 1 and the SNR is set by SNR = 20 dB in this example. The parameters ε is set to 100 for the MD/MM algorithm. The weight factor is selected by 1 for the WCMSR method.
In first part of simulations, we consider the comparison of the DOA estimates of two sources from the directions of −10 • and 10 • using the proposed algorithm and previous methods. First, we fix the inter-spacing of the ULA d at half-wavelength corresponding to the highest signal frequency. We then enlarge d to 100 times to demonstrate the aliasing-free property of the proposed algorithm. In Fig. 3, d equals to the half-wavelength with respect to the highest signal frequency and d = πc/ ω 0 + B f 2 , where c = 3 × 10 8 is the propagation speed of the signal, ω 0 = 1.4π × 10 8 rad/s is the central frequency and B f = 0.2ω 0 is the bandwidth. The MUSIC method decomposes the wideband signals into narrowband subbands first and obtains DOA estimates for each frequency bin. Thus, we plot 3D view of spatial spectra with respect to various frequency components and angle grids; see Fig. 3. In Fig. 4, we compare the results of the proposed method with MUSIC, MD/MM, and WCMSR algorithms. The 3D spectrum of the MUSIC method is averaged in angle domain for ease of comparison. All these methods can detect the DOAs successfully for this example.
When the inter-spacing of the ULA is expanded 100 times, the phenomenon of spatial aliasing appears for MUSIC and WCMSR algorithms, yet the MD/MM and SIF-OMP methods do not suffer spatial aliasing problem; see Figs. 5 and 6. Figure 5 plots the 3D spectrum of the MUSIC method. The WCMSR also has the spatial ambiguity problem because it works in time domain. Both MD/MM and SIF-OMP algorithms can overcome the spatial ambiguity for this example, and they detect the spectra correctly. With the interval of two adjacent frequencies ω = 2π, the spatial nonambiguity of the proposed algorithm is guaranteed if d < c 2 according to Theorem 1 in [14].
We compare the root mean-squared errors (RMSEs) of the proposed approach with MD/MM and WCMSR algorithms. Please refer to Figs. 7,8,9,and 10. The RMSE is defined as follows: where E{·} represents the expectation operation. First, we fix the inter-spacing of the ULA at half-wavelength with respect to the highest frequency. The RMSEs of the three methods are summarized in Fig. 7. All methods have satisfactory performance for high SNR. However, the performance of the SIF-OMP method gets worse at low SNR region due to the additional errors. We also examine the RMSEs of the three methods against the number of utilized frequencies when the SNR is fixed at 0 dB; see Fig. 8. Note that the RMSE of the SIF-OMP algorithm does not reduce as the number of frequencies increases. That is because the total value of errors can not be decreased with the increase of frequency bins. We then enlarge the Fig. 3 The 3D spectral-spatial spectrum of MUSIC with half-wavelength. The 3D spectral-spatial spectrum of MUSIC with two BPSK signals from −10 • and 10 • impinging on a ULA interspaced by half-wavelength relative to the highest frequency, SNR: 20 dB, K = 256 and M = 1 Fig. 4 The spatial spectra of two sources with half-wavelength. The spatial spectra of two BPSK signals from −10 • and 10 • on a ULA interspaced by half-wavelength relative to the highest frequency, SNR: 20 dB, K = 256 and M = 1 Fig. 5 The 3D spectral-spatial spectrum of MUSIC with 100 times half-wavelength. The 3D spectral-spatial spectrum of MUSIC with two BPSK signals from −10 • and 10 • impinging on a ULA interspaced by 100 times half-wavelength relative to the highest frequency, SNR: 20 dB, K = 256 and M = 1 Fig. 6 The spatial spectra of two sources with 100 times half-wavelength. The spatial spectra of two BPSK signals from −10 • and 10 • on a ULA interspaced by 100 times half-wavelength relative to the highest frequency, SNR: 20 dB, K = 256 and M = 1 inter-spacing of the ULA to 100 times half-wavelength with respect to the highest frequency and keep the other settings unchanged. Figure 9 plots the RMSEs of the MD/MM and SIF-OMP algorithms with respect to various SNR. Again, we observe that the additional errors affect the performance of the SIF-OMP method for low SNR. Figure 10 shows the RMSEs of the MD/MM and SIF-OMP algorithms against the number of frequency bins when SNR equals 0 dB. These two algorithms both have good performance. We then derive the separation probabilities of those three algorithms in Figs. 11, 12, and 13. A successful separation is defined similarly with [15] when two conditions are satisfied. (1) The amplitude of the highest pseudopeak is lower than half of that of the lowest signal peak and (2) the bias of the DOA estimates is less than 3 • . In the first group of examples, the first source signal is fixed at −10 • , and the second one varies from −4 • to 10 • . For each angle separation, the experiments are done by 100 trials. Figure 11 illustrates the separation probabilities of MD/MM, WCMSR, and SIF-OMP at different angle separations. We observe from Fig. 11 that WCMSR gets better results than MD/MM and SIF-OMP when d equals to half-wavelength relative to the highest frequency, whereas the MD/MM and SIF-OMP algorithms can separate two BPSK signals with almost 100% probability when d is expanded to 100 times half-wavelength relative to the highest frequency. We then fix the directions of two signals at −10 • and 0 • and compare the separation probabilities with respect to various SNRs and different number of frequency bins. First, K is fixed at 256 and SNR varies from −10 to 10 dB. The separation probabilities of MD/MM, WCMSR, and SIF-OMP are drawn in Fig. 12. Again, we can conclude from Fig. 12 that it is more potential to obtain 100% separation probability when the inter-spacing of the ULA is much larger than half-wavelength. If the SNR is fixed at 0 dB and the number of frequency bins varies from 128 to 1280, their probabilities of separation are plotted in Fig. 13. The angular separation performance of SIF-OMP can be improved when the number of utilized frequencies increases.

Conclusions
In this paper, we present a new aliasing-free SIF algorithm to solve the wideband DOA estimation problem. The SIF algorithm utilizes all frequency bin information to recover a SIV at the expense of introducing an additional system error. By doing this, we reformulate the MD/MM problem as a single measurement vector recovery problem and therefore reduce the unknown variables greatly. We show that the additional system error can be neglected under certain conditions. After using the binary property of the SIV, we develop a SIF-OMP algorithm to estimate the SIV. We compare the performances of the SIF-OMP algorithm