DOA estimation of wideband signals based on slice-sparse representation

In this article, the direction-of-arrival (DOA) estimation problem of wideband signal sources is studied. We pass the incident signals through a bank of narrowband filters to split the array outputs into several narrowband components. Then, a novel slice-sparse representation model of the joint narrowband array covariance data is proposed in the frequency domain to enforce joint sparsity in the concatenated covariance matrix of all frequencies. Based on the greed matching pursuit algorithm, a multiple measurement slices orthogonal matching pursuit algorithm is proposed to exploit the joint frequency processing in the case of wideband scenarios. The DOA estimation is achieved by joint processing of the array covariance data at different frequency bins. The estimated performance is compared with the representative DOA estimation methods. Simulation experiments are conducted to validate the effectiveness of the proposed method.


Introduction
Due to the increase use of wideband signals in the fields of wireless communication system and radar, the problem of direction of arrival (DOA) estimation of wideband signals has been of considerable interest to the array signal processing in recent years. Many methods have been proposed to estimate the DOAs of wideband signals, among which the maximum-likelihood (ML) methods and subspace methods are most studied. The ML estimators show excellent performance but it needs multidimensional nonlinear global search [1]. The subspace methods, although not optimal, are computationally more attractive than ML methods. The subspace methods can be classified into two major categories: incoherent signal-subspace method (ISSM) [2] and coherent signal-subspace method (CSSM) [3][4][5][6]. These methods decompose the incident wideband signals into narrowband components by passing them through a bank of narrowband filters, and then obtain the DOA estimation incoherently or coherently. The ISSM incoherently constructs the final result by taking an average of different frequency bins. Although it *Correspondence: ganlu@uestc.edu.cn School of Electronic Engineering, University of Electronic Science and Technology of China, Chengdu, Sichuan, 611731, China works well at high signal-to-noise ratio (SNR), its performance degrades greatly at low SNR or the SNR varies at different frequency bins because even a single outlier from one frequency bin can potentially lead to inaccurate estimate through the averaging processing. The CSSM surpasses the ISSM because of the capability to process the coherent sources by employing the spectral focusing technique [3]. However, the accuracy of preliminary DOA estimation which is required for the spectral focusing largely influences the performance of CSSM [6]. Therefore, we must explore better solutions for the problem of wideband DOA estimation.
In recent years, the sparse signal representation has attracted enormous attention. By now, there has been several articles addressed sparse representation (SR)-based DOA estimation [7][8][9]. Among them, the l1-SVD method proposed by Malioutov et al. [8] is a successful attempt. This method has also been extended to the wideband signals by treating each frequency band independently [8]. But, repeated using convex optimization technique leads to heavy computational burden in the multiple measurement vectors (MMV) scenario. The wideband covariance matrix sparse representation (W-CMSR) method presented by Liu and Haung [10] forms a new measurement vector by aligning the lower left triangular elements of the array output covariance matrix and reconstructs this http://asp.eurasipjournals.com/content/2013/1/18 vector on an over-complete dictionary to realize wideband DOA estimation. But, it relies on some prior information such as signal modulation types and the pre-estimated signal power spectrum. These give rise to performance degradation when the incident signals share different types of modulations or signal power spectrums. Therefore, further efforts are required to derive integrative wideband DOA estimates method.
In this article, we still address the wideband DOA estimation problem by dividing the wideband data into narrowband counterparts. To work in the frequency domain, we propose a slice-SR model of the joint array covariance data that are stacked as a tensor matrix. Based on the greedy matching pursuit algorithm [11,12], we propose a multiple measurement slices (MMS) orthogonal matching pursuit (MMS-OMP) algorithm. This algorithm processes the narrowband covariance data jointly to obtain the spatial-frequency spectrum. The proposed method and other representative methods are compared. The simulation results show that this algorithm has better noise robustness, higher resolution, and can resolve coherent sources.
This article contributes to the field of wideband DOA estimation in the following four aspects, (1) We use the idea of joint processing of the narrowband components, so the spectral focusing is not introduced. (2) The proposed method is based on the greedy pursuit algorithms so that, it is a low-complexity and high resolution estimator for the wideband signal sources. (3) An important benefit that comes with our algorithm is the ability to incorporate prior information on the frequency spectra of the sources. (4) Furthermore, we can allow the distance between adjacent array elements larger than the smallest half-wavelength with respect to the wideband source frequency, just like what W-CMSR has achieved.
The remainder of the article mainly consists of four sections. In Section 2, we describe the frequency domain model of wideband DOA estimation. In Section 3, we present the viewpoint of the SR of array output covariance vectors. Then, a concept of joint K-slice sparse is proposed for the stack operation of the MMS in the form of tensor matrix. Finally, we formulate the MMS-OMP algorithm in detail. Some numerical results are provided to demonstrate the performance and computational efficiency of MMS-OMP in Section 4. In Section 5, we make a conclusion.

Data model
For wideband signal sources, suppose that K far-field wideband signals impinge onto an N-element array from the directions of θ 1 , θ 2 , . . . , θ K respectively, which are corrupted by additive complex Gaussian white noise. At time t, the data of the n th array element are collected by where s k (t) is the waveform of the kth source, τ nk is the propagation delay of the nth element of the array with respect to the reference element, w n (t) is the zeromean complex additive Gaussian white noise. We assume that the incident signals and the additive noise are mutually independent. The number of sources K can be estimated by Akaike information criterion (AIC) or minimum description length (MDL) criterion [13].
Then, the observation time T is divided into L subsegments. Each sub-segment has an observation time T d , apparently, T = LT d . We transform the T d received data of the lth segment into the frequency domain resulting in J non-overlapping narrowband components. We denote the frequency domain data at frequency f j by where are derived from the discrete Fourier transform (DFT) of the received data, signals and noise, respectively. The array manifold matrix with respect to the above wideband model is Considering a uniform linear array (ULA), the array steering vector a k (f j ) can be expressed as follows where τ nk = 1 c (n − 1)d sin(θ k ), d is the distance between adjacent array element, c is the velocity of wave propagation.
Equation (2) is the frequency domain model of wideband signal sources. It is similar to time domain model of narrowband sources in the form. The estimated array covariance matrix at frequency f j is given bŷ whereR(f j ) ∈ C N×N , [ ·] H denotes the complex conjugate transpose operator. The wideband DOA estimation method proposed in next section can be derived from SR of the array output covariance vectors inR(f j ). http://asp.eurasipjournals.com/content/2013/1/18 3 Proposed method

SR of the array covariance vectors
For the frequency f j , (5) can be reformulated as followŝ is the noise power, I is an N × N identity matrix. In order to a convenient derivation, we temporarily neglect the noise item. Then we havê According to SR theory [14], it is obviously that (7) can be rewritten asR We representR(f j ) in the form of column vectors aŝ . So each column vector inR(f j ) can be regarded as the linear combination of the array steering vectors in A(f j ), i.e., where g kn is the kth row nth column element of G. To cast this problem as a SR problem, we introduce an overcomplete representation j in terms of all possible DOAs. Then, let {θ 1 ,θ 2 , . . . ,θ N s } be a sampling grid of all DOAs of interest. The number of potential DOAs N s will typically be much greater than the number of sources K. Then, we construct the over-complete basic matrix composed of steering vectors corresponding to each potential DOA as its columns So (8) can be represented aŝ where the linear combination coefficient b n s satisfy Sor n (f j ) is the linear combination of the range space of j . We reformulate (10) as the form of matrices, then we obtain where b n,j is the N s × 1 representation coefficient vector with respect to the above over-complete basis at frequency f j . Since j has a nontrivial nullspace, many solutions of b n,j can fitr n (f j ) well. To make b n,j unique Ksparse solution, we should impose the sparsity of the signals relative to the whole spatial domain. If {θ 1 ,θ 2 , . . . ,θ N s } are dense enough, some K angles in {θ 1 ,θ 2 , . . . ,θ N s } can be expected to be very close (or even equal) to the actual {θ 1 , θ 2 , . . . , θ K }, so an ideal b n,j should be the vector with all elements zero except for the K elements associated with these K angles, which are related to the DOAs of the signals [14]. Then, we rewrite the model in (7) in the form of matrix as followŝ where . . , N should share the same sparse structure, and the nonzero elements of each ideal b n,j should occur in the same rows of B j . In other words, only K rows of b n,j are nonzero. Such a matrix is called joint K-rows sparse [12]. We callr n (f j ), n = 1, 2, . . . , N as MMV.
Slices are the two-dimensional sections of a tensor matrix. Figures 1 and 2 show the horizontal, lateral, and frontal slices of a tensor matrix [15]. The above analyzes the SR of the array output covariance vectors at frequency f j . For all narrowband data of every frequency bins, we have MMS, which consist of corresponding array output covariance vectors. Just as shown in Figure 1, the MMS with respect to corresponding frequency bins are stacked up as frontal slices. As a result, B j , j = 1, 2, . . . , J are stacked up in the form of three-dimensional (3-D) tensor matrix B. The new tensor matrix B ∈ C N s ×N×J illustrated in Figure 3 has a property of nonzero elements occupying only K horizontal slices, so that it is called joint K-slice sparse. Therefore, the DOA estimation problem can be posed as the MMS problem of finding a slice-sparse tensor matrix B to represent the 3-D measurement matrix R ∈ C N×N×J under the tensor basic matrix ∈ C N×N s ×J .
For the above joint K-slice sparse model, we formulate our algorithm based on the fundamental OMP algorithm. In order to facilitate the presentation before introducing our algorithm, we define some concepts as follows.
(1) · 2 F and · 2 denote the Frobenius norm for a matrix and the Euclidean norm for a vector; Figure 2 denotes the residual matrices after i th iteration, where R which is related to the frequency f j . Its orthogonal complement is

OMP algorithm
Orthogonal matching pursuit algorithm [11] can process the underdetermined equation solution in (12), i.e., r n (f j ) = j b n,j . Since b n,j has only K nonzero components, the vectorr n (f j ) = j b n,j is a linear combination of K columns from j . To identify the ideal b n,j , we need to determine which columns of j participate in the measurement vectorr n (f j ). The idea behind the algorithm is to pick columns in a greedy fashion. At ith iteration, we choose the column s i of j that is most strongly correlated with the remaining part ofr n (f j ) by wherer (i−1) j is the (i−1)th residual andr (0) j =r n (f j ). Then we subtract off its contribution tor n (f j ) byr Then, we iterate on the residual.
One hopes that, after K iterations, the algorithm will have identified the correct set of columns.

MMS-OMP algorithm
Considering the standard OMP algorithm [11], the sparse solution can be found by building up a small subset of column vectors selected from j (we denote it with support set) to representr n (f j ) sequentially. The selection of a nonzero row of b n,j equals to selecting a column of j .
In the joint K-slice sparse scenario, we use an extensional criterion to select the support set. Consequently, favorable performance and lower computational complexity are achieved. So our method presented in the followings is motivated by and is extensions of method developed for J = 1 and n = 1. In the MMS-OMP algorithm, we first find the lateral slice in the tensor matrix , which is best aligned with the MMV (i.e., r j , and the residual R (1) is obtained by stacking each residual R (1) j from front to back. Next, the lateral slice ψ s 2 in , which is best aligned with R (1) , is found, and a new residual R (2) is formed. The algorithm proceeds by sequentially choosing the lateral slices that best matches the residual matrix. Now, we go into the detail of the MMS-OMP algorithm by looking at the ith iteration.
In the ith iteration, we find the slice matrix most closely aligned with the residual R (i−1) by examining the total matching error where (i) is the jth sub matching error, a s,j is the sth column vector of j . The vector which minimizes the Frobenius norm of the jth sub matching error is selected, i.e., The minimization is achieved by maximizing the second term in the above expression. P a s,j = a s,j a H s,j is employed to select the index of lateral slice as The support set is update by where P i j is the projection matrix of i j . The MMS-OMP algorithm is mainly defined by (19) and (20). The main steps of the MMS-OMP algorithm can be described as follows.

Input:
Tensor matrix R; Tensor basic matrix ; Iterations K. Output: Tensor matrixB, The spatial spectrum P mms omp . Initialization: R (0) = R Initial residue; S 0 = ∅ Initial support set; 0 = ∅ Initial atom set; i = 0 Initial iteration counter. (2) Update the support set and the selected slices (4) If i < K, i = i + 1 and return to step 1), otherwise quit the iteration. We will quit the iteration if halting condition is true. For example, we just need K iteration, since the number of incident signals is known. The algorithm finally obtains the support set S K . On one hand, we can transform it into angle value directly. The kth DOA can be computed by θ k = θ min + (s k − 1)h, where θ min and h are the minimum angle and the step size in the direction grid, respectively. On the other hand, the estimated tensor matrixB can be employed to form the spatial-frequency spectrum, and then the spatial spectrum is obtained via the average of each frequency counterpart, i.e., Then, the DOAs can be found by searching the spectral peak.

Simulation results
In this section, some simulations are carried out to check the performance of MMS-OMP in wideband DOA estimation. The proposed method is compared with the previous methods, including two conventional subspace algorithm-MUSIC [2], CSSM [5] and a kind of SR algorithm-W-CMSR [10], and W-CMSR(SK) [10] represents that the prior spectrum information on W-CMSR is known.
In the simulation I, an intuitionistic spectrum comparison is given between MMS-OMP and previous methods. Assume that there are two far-field uncorrelated wideband signal sources impinging onto an ULA with N = 10 array elements from the directions of 10 • and 30 • . The incident signals are BPSKs with the central frequency 10 MHz and the bandwidth 4 MHz, so that the ratio of bandwidth to central frequency is 40%. The distance between adjacent array elements equals the half-wavelength with respect to the highest signal frequency. The distance between adjacent array elements equals the half-wavelength with respect to the highest signal frequency. We take N s = 360 by searching the angle space from −90 • to 89.5 • with step size 0.5 • . The SNR is set to −10 dB and the number of time domain snapshots is set to 80. We do the temporal-spectral transformation and collect the spectral snapshots by decomposing the sources into 80 frequency bins for MUSIC, CSSM and MMS-OMP. We retain the whole spatial-frequency spectrum of these two methods (MUSIC and MMS-OMP) in Figure 4a,b to give a more composite view of these methods. The 3-D spectrum of MUSIC shows that the peak locations and amplitudes for different frequency points are quite different, due to the fact that it treats each frequency independently. To mitigate this artifact, the MMS-OMP algorithm can incorporate a prior on the continuity of the frequency spectrum of the sources because of the joint frequency processing technique. Then, we average the 3-D spectrum of MUSIC and MMS-OMP in the frequency domain and reprint the corresponding 2-D spectrum in Figure 4c for a more convenient comparison with CSSM and W-CMSR. The results indicate that MUSIC, CSSM, W-CMSR all cannot express angular distribution correctly and pseudo-peaks exist in the spectrum of the three while the MMS-OMP algorithm can gain two correct peaks at 10 • and 30 • .
In the simulation II, we examine the performance of MMS-OMP in the presence of spatial aliasing. The interspacing of the array is enlarged to 1.5 times wavelength associated with the highest signal frequency. 1,200 time domain snapshots are collected and the SNR is 10 dB. Figure 5 shows that MUSIC and CSSM suffer from severe aliasing, and are unable to give any meaningful results. MMS-OMP can suppress the aliasing effects by employing the joint frequency processing technique over the whole frequency range.
In the simulation III, we consider the root mean square error (RMSE) of MMS-OMP and several existing DOA estimation methods, in the statistical sense for signals of various snapshots and SNRs. The curves were obtained by averaging the results of 500 Monte-Carlo simulations. We assume that two BPSK signals with bandwidth of 40% impinge onto the 10-ULA and the direction of the signal is fixed at 10 • and 30 • . Firstly, we compare the RMSE of the four methods in the case of collecting 1,200 snapshots. We can see that the RMSE of the proposed method are smaller than the others in different SNRs indicated in Figure 6. Although W-CMSR(SK) makes full use of prior spectrum information such as code-rate for BPSK, its RMSE is higher than MMS-OMP. These confirm that our method is a more robust algorithm against DOA estimation error. Secondly, we give a comparison of the RMSE http://asp.eurasipjournals.com/content/2013/1/18 by fixing the SNR at 10 dB when the snapshots vary from 320 to 1,280. Figure 7 shows that our method also has a good adaptability to snapshots.
Simulation IV compares the estimated probability of the four algorithms with the change of SNR, which is derived from 500 Monte-Carlo simulations whose condition is the same as simulation III. Our method gives higher estimated probability than other algorithm when the SNR  is deteriorated just as illustrated in the introduction. The proposed method surpasses the W-CMSR algorithm by the absolute superiority in Figure 9, because our method relies less on a priori information of the incident signals than that of the W-CMSR algorithm. Finally, the average computation time of those four methods are compared in different snapshots and signal bandwidths in the simulation V. In order to get a more conveniently adjustable bandwidth, we assume that two 0 dB signals which are formed as mixtures of equal-power sinusoids impinge onto the array from directions of 10 • and 30 • simultaneously. The signal bandwidth is set at 20, 40, 80%, and the time domain snapshots are set at 80, 320 and 1,280. For each bandwidth-snapshot pair, 500 Monte-Carlo trials are carried out to obtain the average computation time of those four methods, and the results are given in Table 1. The results show that MUSIC and CSSM are the two fastest methods. W-CMSR consumes large amount of computation time although its computational efficiency is not affected by the signal bandwidth. And MMS-OMP outperforms W-CMSR significantly in computational efficiency. Although MMS-OMP has larger computation time than the conventional methods, it is more important to note that it performs much better in DOA estimation.

Conclusion
In this article, we are engaged in the DOA estimation of wideband signal sources. We stack all the array covariance matrices, which are derived from passing the observed signal through a bank of narrowband filters, in the form of 3-D tensor matrix to put forward the slice-SR model in the frequency domain. We formulize the corresponding algorithm called MMS-OMP algorithm. This algorithm processes the array covariance data of each frequency bin jointly. The DOA estimation is achieved by jointly finding the sparsest coefficients of the MMS under the basis matrices. Its superiority embodies not only in the capability to resolve coherent sources, but also to incorporate prior information on the frequency spectra of the sources. Furthermore, the MMS-OMP is very efficient in relax the aliasing effects when the distance between adjacent array element is larger than the half-wavelength with respect to the highest frequency. The proposed method and the representative methods are compared. Simulation results show that our method has better noise robustness, higher resolution, and smaller computational complexity.