Reduced-Complexity Direction of Arrival Estimation Using Real-Valued Computation with Arbitrary Array Configurations

. A low-complexity algorithm is presented to dramatically reduce the complexity of the multiple signal classi ﬁ cation (MUSIC) algorithm for direction of arrival (DOA) estimation, in which both tasks of eigenvalue decomposition (EVD) and spectral search are implemented with e ﬃ cient real-valued computations, leading to about 75% complexity reduction as compared to the standard MUSIC. Furthermore, the proposed technique has no dependence on array con ﬁ gurations and is hence suitable for arbitrary array geometries, which shows a signi ﬁ cant implementation advantage over most state-of-the-art unitary estimators including unitary MUSIC (U-MUSIC). Numerical simulations over a wide range of scenarios are conducted to show the performance of the new technique, which demonstrates that with a signi ﬁ cantly reduced computational complexity, the new approach is able to provide a close accuracy to the standard MUSIC.


Introduction
Direction-of-arrival (DOA) estimation is a fundamental array processing problem with many applications, for example, radar, sonar, wireless communication, and smart antenna design [1,2].This topic has been extensively studied during the past four decades, resulting in many efficient and accurate algorithms including multiple signal classification (MUSIC) [3], maximum likelihood (ML) [4], subspace fitting (SF) [5], minimum norm (MN) [6], and their derivations.Among those methods, the MUSIC algorithm is able to provide super resolution DOA estimation for closely spaced sources and has an easy implementation with arbitrary array configuration [7].However, such advantages are archived at the expense of huge complexity which is mainly caused by an involved eigenvalue decomposition (EVD) step and a tremendous spectral search step, both are generally implemented with complex-valued computations.Therefore, the standard MUSIC is prohibitively expensive in scenarios where real-time processing is required [8].
It is well known that the problem of DOA estimation can be simplified with some specified array structures, based on which a variety of low-complexity methods has been proposed for high computational efficiency.For example, when the array geometry satisfies the rotational invariant property, DOA can be immediately computed without spectral search by using the estimation of signal parameters via rotational invariance technique (ESPRIT) [9].Exploiting the special geometry of a uniform linear array (ULA), one can also find source DOA by polynomial rooting with search-free computations using the well-known root-MUSIC algorithm [10].Besides, numerous useful modifications of ESPRIT and root-MUSIC have been proposed with circular and other kinds of arrays [11,12].As the array structures used in these techniques are still rather specific, many promising attempts such as array interpolation (AI) [13], manifold separation technique (MST) [14], and Fourier-domain (FD) root-MUSIC [15] have been proposed recently to extend ESPRIT and root-MUSIC to more general classes of array configurations.Unfortunately, those approaches usually use a polynomial with a sufficiently high order to warrant that the truncation errors are small, and hence, the complexity for finding the roots of this polynomial may be higher than expected [16].
Considering one multiplication between two complex variables generally require four times that between two real ones, algorithms with real-valued computation can reduce about 75% complexity as compared to their complex versions.Following this idea, Huarng and Yeh [17] introduced an outstanding unitary transformation for low-complexity DOA estimation with real-valued computation.Whereafter, Linebarger et al. [18] proposed a forward/backward (FB) averaging method which obtains real-valued arithmetic for fast DOA estimation as well.These techniques are firstly exploited by the unitary MUSIC (U-MUSIC) algorithm [19] and then extended to many other derivations including unitary root MUSIC (U-root-MUSIC) [20], unitary ESPRIT (U-ESPRIT) [21], unitary method of direction-of-arrival estimation (U-MODE) [22], and unitary matrix pencil (U-MP) [23].Based on the special structure of a centrosymmetrical array (CSA), those algorithms transform the complex array covariance matrix into a real one.It has been proven that this real matrix is symmetrical, and hence, both matrix decomposition and spectral search can be implemented with real-valued computation.Besides, it has been found that real-valued algorithms also show improved accuracy as compared to complex-valued approaches.Despite their increased estimation accuracies with reduced costs, almost all of the state-of-the-art real-valued methods are only suitable for CSAs [24], which severely limits their applications [25,26].
As a further development of our previous work in [27], the objective of the present paper is to propose a new real-valued subspace-based algorithm for fast DOA estimation.However, unlike the abovementioned methods that attempted to overcome the high-complexity problem using special array structures or array mapping techniques [28], we show that exploiting only the real part of the array covariance matrix can lead to a real-valued version of the MUSIC algorithm, which has implementation advantages over most state-ofthe-art real-valued techniques since it can be used with arbitrary arrays.Furthermore, we prove that the real part of the array covariance matrix can be equivalently reformulated as an entire array covariance matrix received by a virtual array with an extended array manifold which is real.Because the virtual array covariance matrix and the virtual array manifold are real, both tasks of matrix decomposition and spectral search can be implemented with efficient real-valued computations instead of complex-valued ones.Thanks to these merits, the new algorithm is able to reduce about 75% computational complexity as compared to the standard MUSIC.Finally, we provide in-depth insights into the new approach including the EVD of the real part of the array covariance matrix as well as performance simulation over a wide range of application scenarios, which demonstrate that with significantly reduced complexity, the proposed approach can still provide acceptable accuracy close to the standard MUSIC.
1.1.Mathematical Notations.Throughout the paper, matrices and vectors are denoted by upper and lower boldface letters, respectively.Complex and real vectors and matrices are denoted by single-bar-and double-bar upper boldface letters, respectively, and detailed mathematical notations are defined in Table 1.

Signal Model and Standard MUSIC
2.1.Signal Model and Basic Assumptions.Without loss of generality, let us consider an arbitrary linear array (as to be shown that the proposed method has no dependence on array geometries, and it can be easily extended to arbitrary plane array for two-dimensional DOA estimate) with M omnidirectional sensors labeled by 1, 2, … , M, where the coordinate of the mth sensor is denoted by x m , m = 1, 2, … , M. Suppose that there are L uncorrelated narrow-band plane waves with unknown DOAs Θ ≜ θ 1 , θ 2 , … , θ L impinging on the array from far field, the array received data can be written as where N is the number of snapshots, is the additive white Gaussian noise (AWGN), and is the steering vector matrix (also known as array manifold), and each column of A Θ is defined as a so-called steering vector, which can be expressed as a θ ≜ e j⋅2π/λ⋅x 1 ⋅sin θ , e j⋅2π/λ⋅x 2 ⋅sin θ , ⋯, e j⋅2π/λ⋅x M ⋅sin θ T ∈ ℂ M×1 , 3 where j = −1 and λ are the center wavelength.The goal of the DOA estimation is to find the DOA set Θ, based on N snapshots of array received data x t N t=1 .Throughout the paper, it is assumed that the number of sources is known [29] and is smaller than half that of sensors such that L ≤ M − 1 /2 .It is also assumed that positions and displacements of array elements do not need to be uniform or follow regular patterns, provided that they satisfy the rank M − 1 ambiguity restriction [30].The AWGN n t is assumed to satisfy the following statistical performances [31][32][33].
and the incident signal s t is uncorrelated and satisfies [34,35] where R s is the L × L source covariance matrix and is a real diagonal matrix.
and further obtained the EVD of R as follows: where span S and span G are the signal and noise subspaces, respectively.In practical situations, the theoretical R is estimated by Hence, the EVD of the array output covariance matrix is in fact given by Based on some facts span S ⊥span G and span S = span A , the standard MUSIC suggests to search the L peaks of the following function over −π/2, π/2 to find source DOAs. 10 1 The primary advantage of the MUSIC algorithm is its easy implementation with arbitrary array configurations, which is archived at the expense of huge complexity caused by an involved tremendous spectral search step with expensive complex-valued computation.

The Proposed Algorithm
We begin our algorithm by considering the M × M real part of the array output covariance.Using (6), we have where B Θ is defined as Therefore, we have where A Θ c and R s ~are given by Using (5) 16) into (11), we obtain where x t ~is given by Comparing ( 18) and ( 19) with ( 1) and ( 6), the real part of the array covariance matrix can be regarded as an entire array covariance matrix of a virtual array with a manifold A Θ .The output sensor vector of this virtual array is x t ~and the incident signal on this virtual array is s t ~with its covariance matrix given by ℝ s ~.Therefore, (19) provides a virtual signal model for DOA estimation, which is similar to the conventional one given by (1).Because DOA information is also contained in A Θ , one can exploit the idea of eigenvalue analysis to perform the EVD of the virtual array covariance matrix R and construct a MUSIC-like spectral for DOA estimate.Noting that R T = R, which implies that R is a symmetrical real matrix.Therefore, the EVD computation requires only real-valued flops [36].Using (6), it is shown in the appendix that the eigenvalues of R are given by and ℍ Θ is defined as According to (21), R has 2L − K significant and M − 2L + K smallest eigenvalues.Therefore, (20) can be further written as where ℚ s contains the 2L − K significant eigenvalues and S and G are the signal and noise matrices such that

25
Since A Θ is the virtual manifold matrix associated with R, by using the orthogonality between the signal and the noise subspaces, we obtain Because A Θ = Re A Θ Im A Θ is in fact a combined manifold matrix covering both Re A Θ and Im A Θ , we have where a re θ ≜ Re a θ ∈ ℝ M×1 and a im θ ≜ Im a θ ∈ ℝ M×1 are two real vectors, given by Based on ( 26) and ( 27), we can define a MUSIC-like spectral for DOA estimation as follows Since a re θ , a im θ , and G are all real, we must have which means that the minima of f Proposed θ over half of the angular field of view, that is, −π/2, 0 or 0, π/2 , is either true 4 International Journal of Antennas and Propagation DOAs or their images.Thus, we can search over half of the angular field of view 0, π/2 to obtain T candidate angles and select the L true DOAs from Θ + and its mirror The conventional beamformer (CBF) [37] is exploited to select the L true DOAs among Θ + and Θ − .As the CBF takes its maxima only at the true DOAs, the CBF spectral amplitudes responding to the true DOAs must be much larger than those associated with their symmetrical mirrors.On the other hand, as number of the true DOAs is known in advance, the L true DOAs can be easily selected among the candidate angles by where f CBF θ ≜ 10 log 10 a H θ Ra θ dB 35 Although using CBF to exclude the symmetrical mirror DOAs means that there is an additional step involved in the proposed estimator, the complexity of this further step is substantially lower since we only need to compute the product a H θ Ra θ for at most 2L spectral points.Detailed for implementing the proposed algorithm are summarized in Algorithm 1.
We end up this section by a comparison of primary realvalued computational flops required by MUSIC [3], U-MUSIC [19], and the proposed method, which is shown in Table 2.In the table, J stands for the total number of search points over −π/2, π/2 .Because EVD and spectral search dominate the total complexity, the common flops required by the three algorithms for the estimation of R and the slight flops required by the proposed method for solving the ambiguity problem are omitted.
Using the fast subspace decomposition (FSD) technique [38], the EVD computation on a M × M real matrix for a signal subspace of dimension L costs M 2 L + 2 real-valued flops.As the proposed method has 2L − K equivalent signals, this term for the proposed method is given by M 2 2L − K + 2 .Since the standard MUSIC needs to compute the EVD of the complex matrix R, this term for MUSIC is given by 4 × M2 L + 2 .Because a θ ∈ C M×1 and G ∈ ℂ M× M−L , computing the MUSIC spectral for J sample points over −π/2, π/2 costs 4 × J M + 1 M − L real-valued flops [39].Since the proposed method involves a compressed search over 0, π/2 or −π/2, 0 , it only needs to compute J/2 sample points.Taking into account that are θ ∈ ℝ M×1 , aim θ ∈ ℝ M×1 and G ∈ ℝ M× M−2L+K , the complexity of spectral search for the proposed method is significantly reduced to J M + 1 M − 2L + K .
It can be concluded from Table 2 that our method has the lowest complexity among the three algorithms.Moreover, because we generally have J ≫ M [8,15], about 3/4 = 75% complexities are reduced by the new method as compared to the standard MUSIC.

Simulation Results
Numerical simulations with 500 Monte Carlo trials are conducted to demonstrate the performance of the proposed method on a linear array composed of M = 10 sensors, where L = 2 sources are considered throughout the simulations.For the root mean square error (RMSE) comparison, the Cramer-Rao lower bound (CRLB) given in [40] is also applied for reference and the RMSE is defined as where θl i is the ith estimated value of the lth signal incident angle θ l .For L sources, the signal-to-noise ratio (SNR) is defined as where P avg = 1/L∑ L l=1 P l denotes the average power of all sources, and P l = E s 2 l t is the power of thelth, l = 1, 2, … , L source.
In the first simulation, we use a ULA to estimate the DOAs of two sources located at θ 1 = 30 °and θ 2 = 40 °.Five algorithms including the standard MUSIC [3], ESPRIT [9], U-MUSIC [19], RV-MUSIC [39], and the proposed method are applied in the simulation.
Figure 1 compares the RMSEs of different algorithms against the SNR, where the number of snapshots is set as N = 200.It is seen from Figure 1 that U-MUSIC is the most accurate one among the four methods as expected, although it can be only used with CSAs.It is also seen from the figure that with a significantly reduced complexity, our method can still provide a close accuracy to the standard MUSIC and the RV-MUSIC algorithms, especially with SNR > 0 dB.In addition, the proposed has a much better accuracy than ESPRIT.
To see more clearly the performance of the new method with a ULA, Figure 2 plots the RMSEs of different algorithms as functions of the number of snapshots, where the SNR is International Journal of Antennas and Propagation fixed as SNR = 5 dB.It can be concluded again from Figure 2 that the proposed method can provide acceptable estimation accuracy which is very close to MUSIC and RV-MUSIC and is much better than ESPRIT.
In the second simulation, we examine fast DOA estimation with nonuniform linear arrays (NULAs), where two closely spaced sources located at θ 1 = 25 °and θ 2 = 28 °are considered.We use a minimum redundancy linear array (MRLA) [41] to assess the performance of the proposed method and compare it to that of the standard MUSIC.Note that neither the shift invariant geometry or the centrosymmetrical structure exists in a MRLA, and hence ESPRIT [9] and most state-of-the-art real-valued methods [19][20][21][22][23] including U-MUSIC cannot be used in such an array directly without array mapping techniques [28].
In Figure 3, we set again the number of snapshots as N = 200 and compare the RMSEs of different algorithms against the SNR.In Figure 4, we fix SNR = 5 dB and plot RMSEs as functions of the number of snapshots, where the amounts of snapshots vary from a wide range over N = 20 to N = 5120.
It can be seen clearly from Figure 3 and Figure 4 that the three methods perform closely to each other.More specifically, the standard MUSIC slightly outperforms RV-MUSIC and the proposed technique, especially with low SNRs and small numbers of snapshots.However, as the proposed method is a real-valued DOA estimator which reduces about 75% complexities, it makes an efficient trade-off between complexity and accuracy as compared to the standard MUSIC.
In the third simulation, we examine the performance of the proposed method with different numbers of signals in Figure 5 and Figure 6, where the simulation parameters are given in the captions.As can be seen clearly from the two figures that our method provides very close RMSEs to

Algorithms
Real-valued computational flops   It is also seen from the two figures that the performances of both the two method decreases slightly when L increases.This is because when L increases, the dimensions of noise subspaces decreases oppositely.
In the last simulation, we use a ULA to plot the simulation times of DOA estimation by MUSIC [3], U-MUSIC [19], RV-MUSIC [39] and the proposed method as functions of the number of sensors in Figure 7, where we set SNR = 10 dB, N = 200, and two sources located at θ 1 = 10 °and θ 2 = 20 °are considered.The simulated results are given by a PC with Intel(R) Core(TM) Duo T5870 2.0 GHz CPU and 1 GB RAM by running the Matlab codes in the same environment, where a common fine search gird Δθ = 0 0053 °is used the three methods.It can be concluded from Figure 7 that although the proposed method is slightly efficient than the RV-MUSIC algorithm, it is the most efficient one among the four techniques, and it costs a simulation time being about 4 times lower than that of MUSIC.

Conclusions
We have proposed an efficient DOA estimation algorithm, in which both tasks of eigenvalue decomposition and spectral search are implemented with low-complexity real-valued computations.In addition, the proposed method can be used

7
International Journal of Antennas and Propagation with arbitrary array configurations, which shows implementation advantage over most existing unitary methods.Complexity analysis and simulation demonstrate that with about 75% reduced complexity, the new method is able to provide a good DOA estimation accuracy which is very close to the standard MUSIC.

Appendix Proof of Equation (21)
We Begin the Proof by Using (6) to Rewrite R as Since there are K pairs of symmetrical sources in Θ, we can divide Θ into two subsets as follows: where Θ 1 contains all the symmetrical sources and we can write Note for any two DOAs Since R s ~is a diagonal matrix with full rank (note that n t is assumed to be uncorrelated), we can write R s ~= VV T , where V is also a diagonal matrix with full rank.Thus, we further obtain [28] rank As β 1 , β 2 , … , β 2L−K are the 2L − K eigenvalues of matrix D Θ , we have Using (A.6), we have which implies that the eigenvalues of R are given by and the proof is completed.

Figure 6 :
Figure 6: RMSE versus the number of snapshots respect to different numbers of signals, ULA, SNR = 10 dB and M = 14.
θT , 32 by a further step to solve the ambiguity problem.Because for one true DOA θ k , f Proposed θ shows two symmetrical peaks at angles θ k and −θ k simultaneously.Therefore, for L true DOAs, f Proposed θ gives at most L peaks over half of the angular field of view.This happens if and only if any two of the L true DOAs are not symmetrical.Oppositely, f Proposed θ generates at least L/2 peaks over half of the angular field of view, and this happens if and only if all the L true DOAs are paired symmetrical themselves.In other words, the number of candidate angels T satisfies