Computationally Efﬁcient Direction-of-Arrival Estimation Based on Partial A Priori Knowledge of Signal Sources

A computationally e ﬃ cient method is proposed for estimating the directions-of-arrival (DOAs) of signals impinging on a uniform linear array (ULA), based on partial a priori knowledge of signal sources. Unlike the classical MUSIC algorithm, the proposed method merely needs the forward recursion of the multistage Wiener ﬁlter (MSWF) to ﬁnd the noise subspace and does not involve an estimate of the array covariance matrix as well as its eigendecomposition. Thereby, the proposed method is computationally e ﬃ cient. Numerical results are given to illustrate the performance of the proposed method.


INTRODUCTION
It is desired to estimate the directions-of-arrival (DOAs) of incident signals from noisy data in many areas such as communication, radar, sonar, and geophysical seismology [1]. The classical subspace-based methods, for example, the MUSIC-type [2] algorithms that rely on the decomposition of the observation space into signal subspace and noise subspace, can provide high-resolution DOA estimates with good estimation accuracy. Normally, the classical subspace-based methods are developed without considering any knowledge of the incident signals, except for some general statistical properties like the second-order ergodicity. Nevertheless, the subspace-based methods typically involve the eigendecomposition of the array covariance matrix. As a result, these methods are rather computationally intensive, especially for large arrays.
To attain better DOA estimation accuracy and, perhaps, reduce the computational complexity, a number of algorithms that assume some a priori knowledge, such as the waveforms, of the incident signals have been developed in [3][4][5][6][7][8][9]. The assumption is reasonable in friendly communications, such as wireless communications and GPS, where certain a priori knowledge of the incident signals is available to the receiver. The a priori information may or may not be explicit. For example, in a packet radio communication system or a mobile communication system, a known preamble may be added to the message for training purposes. In a digital communication system, the modulation format of the transmitted symbol stream is known to the receiver, although the actual transmitted symbol stream is unknown [10]. With the assumption that the waveforms of the incident signals are known, several computationally efficient maximum likelihood (ML) estimators, for example, the methods named DEML [3], CDEML [4], and WDEML [5] were presented for DOA estimation. Using the known waveforms of the signals, these methods decouple the multidimensional nonlinear optimization of the exact ML estimator to a set of one-dimensional (1D) optimization and, thereby, are relatively computationally simple. To reduce the computational complexity, several algorithms for DOA estimation have been developed by exploiting the partial a priori knowledge of signal sources such as the special features of cyclostationary signals [6] and constant modulus (CM) signals [7]. The authors in [6] utilized the cyclic correlation matrix to calculate the noise subspace through a linear operation. Since this method can avoid the eigendecomposition of the covariance matrix, it is computationally efficient. With the CM assumption [7], it is possible to find the estimate of the array response matrix, and then use a scheme similar to the ESPRIT method to directly achieve the DOA estimation, therefore avoiding the 1D search and reducing the computational complexity. Nevertheless, these methods are suitable only for signals with the appropriate special temporal properties. Recently, the reduced-order correlation kernel estimation technique (ROCKET) [11] and ROCK MUSIC algorithms [8,9] were applied to high-resolution spectral 2 EURASIP Journal on Applied Signal Processing estimation. Exploiting the received signal of the first array element to initialize the multistage Wiener filter (MSWF) [12], the ROCKET algorithm only needs the forward recursion of the MSWF to find a subspace of interest and use that subspace to calculate a reduced-rank data matrix and a reducedrank weight vector for a reduced-rank autoregressive (AR) spectrum estimator. Given the direction or spatial frequency of one signal, the ROCK MUSIC method can find a nonunitary basis for the signal subspace by using the forward and backward recursions of the MSWF. The ROCKET and ROCK MUSIC algorithms do not resort to the eigendecomposition of the array covariance matrix, giving them a computational advantage. Nevertheless, the ROCK MUSIC algorithm still needs the forward and backward recursions of the MSWF, which increases the complexity of the algorithm since the backward recursion coefficients completely change with each new stage that is added. To find the reduced-rank data matrix and the reduced-rank AR weight vector, the ROCKET method still involves complex matrix-matrix products, implying that additional computational cost is incurred.
In this paper, we propose a computationally efficient method for DOA estimation, based on partial a priori knowledge of signal sources. Using the orthogonal property of the matched filters of the MSWF, we show that the signal subspace and the noise subspace can be spanned by the matched filters. The estimated noise subspace is then exploited to super-resolve the incident signals instead of using the eigendecomposition-based MUSIC method, thus reducing the computational complexity of calculating the noise subspace. To cure coherent signals, we apply the spatial smoothing technique merely to the array data matrix and the training data vector, and therefore avoid the estimate of the array covariance matrix. Unlike the ROCKET and ROCK MUSIC techniques, the proposed method merely needs the forward recursion of the MSWF to obtain the noise subspace and does not require any complex matrix-matrix products, thereby further reducing the computational complexity of the algorithm. Compared to the classical MUSIC estimator and the fast subspace decomposition (FSD) method [13], the proposed method does not involve the estimate of the array covariance matrix or any eigendecomposition. Thus, the novel method is computationally attractive and can be used in the case of small samples where the array covariance matrix cannot be estimated efficiently. While operationally similar to the classical MUSIC estimator, the proposed method finds the noise subspace in a more computationally efficient way, which is the distinguishing feature of the new method.
This paper is organized as follows. Section 2 gives the data model and reviews the MSWF. Section 3 presents the new method for DOA estimation. In Section 4, numerical results are given. Finally, conclusions are drawn in Section 5.

Data model
Consider a uniform linear array (ULA) composed of M isotropic sensors. Impinging upon the ULA are P narrowband signals from distinct directions θ 1 , θ 2 , . . . , θ P . The M ×1 vector received by the array at the kth snapshot can be expressed as where s i (k) is the scalar complex waveform referred to as the ith signal, n(k) ∈ C M×1 is the additive noise vector, N and P denote the number of snapshots and the number of signals, respectively, a(θ i ) is the steering vector of the array toward direction θ i and takes the following form: where ϕ i = (2πd/λ) sin θ i in which θ i ∈ (−π/2, π/2), d and λ are the interelement spacing and the wavelength, respectively, and the superscript (·) T denotes the transpose operator. Assume that the first signal is the desired signal whose waveform or training data is known. In matrix form, (1) becomes where are the M × P steering matrix and the P × 1 complex signal vector, respectively. Throughout this paper, we assume that M > P. Furthermore, the background noise uncorrelated with the signals is modeled as a stationary, spatiallytemporally white, zero-mean, complex Gaussian random process.

Multistage Wiener filter
It is well known that the Wiener filter (WF) w w f ∈ C M×1 can be used to estimate the desired signal d(k) ∈ C from the array data x(k) in the minimum mean square error (MMSE) sense. Thereby, we have the following design criterion: where d(k) = w H x(k) represents the estimate of the desired signal d(k), and w ∈ C M×1 is the linear filter. Solving (5) leads to the Wiener-Hopf equation where The classical Wiener filter, that is, w w f = R −1 x r xd , is computationally intensive for large M since the inverse of the array covariance matrix R x is involved. The MSWF developed by Goldstein et al. [12] is to find an approximate solution to the Wiener-Hopf equation, which does not need the inverse of the array covariance matrix. The MSWF of rank D based on the datalevel lattice structure [14] is shown in Algorithm 1. Figure 1 illustrates the lattice structure of the MSWF. The reference signal d 0 (k) is the training data of the desired signal, which is available in friendly communications. In this paper, let d 0 (k) = s 1 (k). The observation data x i−1 (k) at the ith stage are partitioned into an interesting signal d i (k) and its orthogonal component The array data matrix is partitioned stage-by-stage in the same manner. As a result, we can readily achieve the prefiltering matrix

COMPUTATIONALLY EFFICIENT ALGORITHM FOR DOA ESTIMATION
Since all the matched filters h 1 , h 2 , . . . , h M are mutually orthogonal for the special choice of the blocking matrix B i = I − h i h H i , the matched filters after the Pth stage of the MSWF are orthogonal to the signal subspace, that is, h i ⊥ col{A(θ)} for i = P + 1, P + 2, . . . , M. Therefore, the last M − P matched filters span the orthogonal complement of the signal subspace, namely the noise subspace: Equation (8) indicates that the noise subspace can be readily obtained by performing the forward recursion of the MSWF, and thus the MUSIC estimator based on the noise subspace can be exploited to produce peaks at the DOA locations. For coherent signals, however, the noise subspace estimated by this method is no longer correct. That is to say, the last M − P matched filters do not span a noise subspace for the case where the signals are coherent. As a result, we must resort to the smoothing techniques to decorrelate the coherent signals. Since the array covariance matrix is not involved in computing the basis vectors for the noise subspace, we perform the spatial smoothing method [16] merely to the array data matrix.
For the spatial smoothing technique, an array consisting of M sensors is subdivided into L subarrays. Thereby The selection matrix J l is used to select part of the M × N array data matrix X 0 = [x 0 (0), x 0 (1), . . . , x 0 (N −1)], which corresponds to the lth subarray. Hence, the spatially smoothed (i) Initialization. d 0 (k) and x 0 (k) = x(k).
(ii) Forward recursion. For i = 1, 2, . . . , D, M L × LN data matrixX 0 is constructed as Similarly to the spatially smoothed data matrixX 0 , the "spatially smoothed" training data vector should have the form ] T ∈ C N×1 and ";" denotes vertical concatenation. Accordingly, the ith spatially smoothed matched filter of the MSWF is computed as Thus, the computationally efficient algorithm for DOA estimation can be summarized as shown in Algorithm 2.
Remark 1. Notice that the lattice structure of the MSWF avoids the formation of blocking matrices, and all the operations of the MSWF only involve complex vector-vector products. Consequently, the proposed method merely requires O(MN) flops to calculate each basis vector h i and thereby Step 1. Apply the spatial smoothing technique to the M × N data matrix X 0 and obtain the spatially smoothed M L × LN data matrixX 0 .
Step 3. Perform the following M L recursions.
Obtain the estimated noise subspace  Remark 2. It should be noted that the proposed method can determine the directions of the desired signal with the knowledge of training data and the interferences without the knowledge of training data. That is to say, the proposed method only needs partial a priori knowledge of signal sources, such as the training data of the desired signal, to estimate the DOAs of all the incident signals.

Uncorrelated signals
Assume that there are two uncorrelated signals with equal power impinging upon the ULA composed of 10 sensors from directions {0 • , 5 • }, and that signal 1 is the desired signal whose waveform is known a priori. We also assume that the number of signals is known. The background noise is a sta- tionary, spatially-temporally white, complex Gaussian random process with zero-mean and the variance σ 2 n . The spatial spectra of the proposed method and the classical MUSIC algorithm are shown in Figure 2, where N = 64 and signal-to-noise ratio (SNR) is 10 dB. SNR is defined as 10 log(σ 2 s /σ 2 n ), where σ 2 s is the power of each signal in a single sensor. From Figure 2, it can be observed that the proposed method works very much like the classical MU-SIC algorithm. To evaluate the estimation performance of the proposed method, we exploit the root-MUSIC algorithm to yield the DOAs of the incident signals and make 500 Monte Carlo runs to compute the root-mean-squared errors (RM-SEs) of estimated DOAs. The RMSEs of estimated DOAs versus SNR are shown in Figure 3, where N = 64. The Cramér-Rao bounds (CRBs) [17] are also plotted for comparison. As shown in Figure 3, when SNR is lower than 6 dB the proposed estimator surpasses the classical MUSIC algorithm, especially in the estimation of the first signal since its waveform is known and used to calculate the basis vectors for the noise subspace. As SNR increases, the proposed method provides the same estimation accuracy as the classical MUSIC algorithm. The RMSEs of the two signals for the two methods approach to the corresponding CRBs when SNR becomes high. The RMSEs of the estimated DOAs for the two methods versus the number of snapshots are demonstrated in Figure 4, where SNR = 5 dB. It can be observed from Figure 4 that the estimation accuracy of the proposed method is higher than that of the classical MUSIC estimator when the number of snapshots is less than 64. As the samples become large, the proposed method yields the same estimation accuracy as the classical MUSIC method.
Lei Huang et al.

Coherent signals
Consider the case where there are two signals impinging upon the ULA consisting of 12 sensors from the same signal source whose waveform is known a priori. The first is a direct-path signal and the other refers to the scaled and delayed replicas of the first signal that represent the multipaths or the "smart" jammers. The propagation constants are {1, −0.8 + j0.6}. We assume that the true DOAs are {0 • , 5 • } and the number of signals is known. The background noise is identical to that in the case of uncorrelated signals. To decorrelate the incident coherent signals, the spatial smoothing technique is also applied to the classical MUSIC algorithm.
The spatial spectra of the proposed method and the classical MUSIC algorithm are shown in Figure 5, where N = 64, SNR = 10 dB, and the number of sensors of the subarray is 9, namely M L = 9. Figure 5 indicates that the proposed estimator works very much like the classical MUSIC estimator in the case of coherent signals. The following results are based on 500 Monte Carlo trials. The RMSEs of estimated DOAs versus SNR are shown in Figure 6, where N = 64. For comparison, the CRBs [18] for coherent signals are given as well. From Figure 6, it can be observed that the proposed method clearly outperforms the classical MUSIC algorithm when SNR ≤ 6 dB, and provides the same estimation accuracy as the latter when SNR > 6 dB. The RMSEs of estimated DOAs for the two methods versus the number of snapshots are plotted in Figure 7, where SNR = 5 dB. It is shown in Figure 7 that the proposed method surpasses the classical MUSIC estimator when the number of snapshots is less than 96 and provides the same estimation accuracy as the latter when the samples become large. Since the waveform of the desired signal is known and exploited to compute the new basis vectors for the signal subspace and the noise subspace, the new signal subspace is capable of capturing the signal information while excluding a large portion of the noise. On the contrary, its orthogonal complement can eliminate the signals more accurately from the noisy data and, thereby, is a 6 EURASIP Journal on Applied Signal Processing   cleaner noise subspace that leads to the enhanced estimation performance.

CONCLUSION
We have presented a computationally efficient method for DOA estimation in this paper. The proposed method only needs the forward recursion of the MSWF and does not resort to the eigendecomposition of the array covariance ma-trix, thereby requiring lower computational cost than the classical MUSIC algorithm especially in the case of a large array. Numerical results indicate that the proposed method surpasses the classical MUSIC estimator for the case of small samples and/or low SNR and provide the same estimation performance as the latter when the samples become large and/or SNR increases.