Abstract
We propose a low complexity subspace-based direction-of-arrival (DOA) estimation algorithm employing a direct signal space construction method (DSPCM) by subsampling the autocorrelation matrix of a uniform linear array (ULA). Three major contributions of this paper are as follows. First of all, we introduce the method of autocorrelation matrix subsampling which enables us to employ a low complexity algorithm based on a ULA without computationally complex eigenvalue decomposition or singular-value decomposition. Secondly, we introduce a signal vector separation method to improve the distinguishability among signal vectors, which can greatly improve the performance, particularly, in low signal-to-noise ratio (SNR) regime. Thirdly, we provide a root finding (RF) method in addition to a spectral search (SS) method as the angle finding scheme. Through simulations, we illustrate that the performance of the proposed scheme is reasonably close to computationally much more expensive MUSIC- (MUltiple SIgnal Classification-) based algorithms. Finally, we illustrate that the computational complexity of the proposed scheme is reduced, in comparison with those of MUSIC-based schemes, by a factor of , where is the number of sources and is the number of antenna elements.
1. Introduction
Subspace-based spectral estimation and direction-of-arrival (DOA) estimation schemes have been widely studied during the last several decades [1–3]. Right after the introduction of Pisarenko and MUSIC (MUltiple SIgnal Classification) algorithms, various alternative schemes such as min-Norm method [4], ESPRIT (Estimation of Signal Parameter via Rotational Invariance Technique) [5], and root-MUSIC method [6] were developed during the first decade. In more recent years, such early ideas and concepts have been extended to deal with various issues such as nonlinear or nonuniform array shapes [7–9], coherent sources [10], and 2-dimensional (2D) angle estimation [11]. Among such issues drawing continued research interest is the computational complexity, which is particularly important for arrays with large number of antenna elements.
One of the main reasons for large computational complexity with existing subspace-based algorithms is the computational burden to construct the subspace, which usually involves the process of eigendecomposition (ED) or singular-value decomposition (SVD). To reduce the computational burden that arises due to ED or SVD, various algorithms that do not require ED or SVD have been proposed. Among the first such examples was the method proposed in [12] which computes the projection matrix by using a submatrix of the correlation matrix. Later, in [13], an improved algorithm called orthogonal propagator method has been proposed that achieves better performance particularly at medium and high signal-to-noise ratio (SNR). Recently, a 2-dimensional (2D) DOA estimation algorithm based on cross-correlation matrix was proposed in [14] that similarly computes the required projection operators without using ED or SVD. However, even though these methods do not require ED or SVD process, they are still computationally burdensome particularly for an array with large antenna elements.
More recently, Xi and Liping proposed, in [15], a very impressive 2D DOA estimation algorithm with L-shaped array. Most of all, it is truly computationally efficient reducing the computational complexity by a factor of , additionally from the algorithms in [12–14]. The authors also introduced a method to exploit the conjugate symmetry to enlarge the effective array aperture, which resulted in reasonably impressive performance. However, despite the fact that the performance of the algorithm [15] was reported to outperform that of [14], its performance is still inferior to computationally more expensive MUSIC scheme particularly at low signal-to-noise ratio (SNR). Moreover, based on the cross-correlation matrix between data collected by two component ULAs of L-shaped array, it cannot be used for a single ULA.
In this paper, we propose a low complexity subspace-based DOA estimation algorithm employing direct signal space construction method (DSPCM) similar to that in [15] but by subsampling the autocorrelation matrix of a ULA. The use of autocorrelation matrix in signal space construction enables us to use the low complexity method with a single ULA and to exploit the Toeplitz structure of the autocorrelation matrix in improving the accuracy of estimating signal vectors. We further propose a method of signal vector separation to improve the distinguishability of the signal vectors and hence to improve the performance particularly in the low SNR regime. Finally, we derive root finding (RF) method in addition to the spectral search (SS) method as the angle finding scheme.
The rest of this paper is organized as follows. In Section 2, the system model is described together with the definition of the notations used throughout this paper. In Section 3, the theoretical background of the algorithm is provided. In particular, a sufficient set of conditions are provided for the applicability of the algorithm. Next, the issues of practical algorithm implementation are discussed in Section 4. In this section, we provide two realizations, namely, SS and RF schemes of the proposed DSPCM. Then, the performance of the proposed scheme is compared with existing schemes, particularly with the MUSIC-based schemes and with that proposed in [15], and the computational complexity is analyzed in Section 5. Finally, we draw conclusions and discuss the directions of future work in Section 6.
Notation. The set of complex numbers will be denoted by and the set of all complex matrices will be denoted by . The complex conjugate, the transpose, and the conjugate transpose of matrix are denoted, respectively, by , and Boldface letters are used for matrices and vectors. We will denote by the matrix whose , and th columns are , and , respectively. identity matrix and zero matrix are denoted, respectively, by and We will denote by the Dirac delta function. For a real number , denotes the largest integer that does not exceed
2. System Model
In this paper, we consider a uniform linear array (ULA) and propose a low complexity subspace-based direction-of-arrival (DOA) estimation algorithm based on correlation matrix subsampling. We assume that antenna elements are located along the -axis at coordinates as shown in Figure 1. We assume that narrowband far-field plane waves of wavelength impinge on the array making angles with the -axis. We let denote the signals corresponding to these waves arriving at the origin, namely, at the leftmost antenna element in Figure 1. We denote by the data collected by the antenna elements located at , respectively. For notational convenience, we let
Here, we assume that represent the complex baseband representation of the relevant signals at time For compact mathematical description, we let and Then, the system can be modeled bywhere represent the additive Gaussian noise vector and is the dimensional array response matrix whose element is given by . We assume that the signal wavelength is known a priori and that the number of impinging signals is predetermined by appropriate methods such as that in [16].
We assume that and are statistically independent complex random vector processes with zero mean satisfying where is the one-sided noise power spectral density and is a positive definite matrix. Throughout this paper, we assume that sources ’s are statistically uncorrelated so that is a diagonal matrix with positive diagonal elements
3. Direct Signal Space Construction by Autocorrelation Matrix Subsampling
In [15], a subspace-based DOA estimation algorithm using the cross-correlation matrix between two data vectors collected by the component ULAs of an L-shaped array has been proposed that does not require computationally expensive ED or SVD in constructing the signal or noise subspace. In this paper, we show that similar signal space construction is possible by subsampling the autocorrelation matrix of a ULA with significantly improved performance.
First, we note that the autocorrelation matrix of is given by In particular, the autocorrelation matrix is a Hermitian Toeplitz matrix with its element given by where Consequently, does not involve noise component except when In particular, we note that for any with values. Noting this fact, we propose methods to construct the signal space by subsampling using only the nondiagonal elements of .
To be more precise, to construct dimensional signal space, we define -dimensional column vectors byHere, we assume that the integers , and are chosen to satisfy the following four conditions: (1) is larger than (2) is larger than (3).(4)None of the indices is zero.
For illustration of the choice of the indices, let us consider the case of , which is, in fact, of our primary interest. In this case, we might define where denotes the integer , namely, the largest integer not exceeding
Next, we define the signal space as the dimensional subspace of spanned by To understand the structure of the signal space , we first note that where is matrix given by and . We recall that the matrix is the diagonal matrix with diagonal elements introduced in Section 2. Consequently, the following relation holds: where is matrix given by Here, we note that the matrix contains the information about the directions of arrivals, namely, Moreover, if is invertible, we have since has already been assumed to be invertible. This means that the columns of containing the information about the directions of arrivals can be obtained as a linear combination of , and . Consequently, the subspace generated by , and can play the role of signal space if is invertible.
It is easy to note that the Vandermonde matrix has rank if and only if are all distinct. However, the rank of the matrix , which is a generalized Vandermonde matrix, is not easily determined in the most general case. However, if we specialize in the case of , which is of particular practical importance, the determinant of is given by which is nonzero if and only if are all distinct. For this reason, we further assume, in what follows, that and are so chosen that the following three properties hold: (1),(2) are all distinct,(3) are all distinct.
Here, the constant will be referred to as the signal vector separation index. We note that the overlap among signal vectors can be reduced by increasing the value of . Even though the signal space dimension is decreased, the increase in the value of can improve the distinguishability among the signal vectors and result in improved performance. We will call this technique signal vector separation.
Now, let us consider dimensional steering vector defined by Then, just as in ordinary subspace-based schemes such as the MUSIC scheme [17], the following four statements become equivalent: (1)The vector belongs to (2)The orthogonal projection of onto equals itself.(3)The augmented matrix , which is a Vandermonde matrix, has rank (4)The complex number equals one of the complex numbers
For the realization of reduced complexity algorithms, we particularly note the equivalence between second and fourth statements. To effectively use this equivalence, we let denote the orthonormal vectors obtained by applying the Gram-Schmidt orthonormalization process on the vectors . Before proceeding, we note that the Gram-Schmidt orthonormalization can be carried out efficiently via matrix QR-decomposition without requiring computationally expensive ED or SVD process. Since the orthogonal projection of is given by it follows that But since with equality if and only if , the following theorem holds.
Theorem 1. The complex number equals one of if and only if takes its maximum possible value .
4. DOA Estimation Algorithm Implementation
4.1. Signal Space Construction and Correlation Estimation
In this paper, we estimate the directions of arrivals using Theorem 1. We note that the direction information is stored in the signal space basis vectors or, equivalently, in the orthogonal basis vectors We recall that the vectors ’s are constructed using nondiagonal elements of the correlation matrix . Consequently, we only need to estimate the values of ’s for nonzero values. Before further proceeding, we note that Consequently, we only need to estimate the values of ’s for positive values.
To discuss how we obtain , we first note for As in most subspace-based schemes, the estimate of can be obtained by Based on this estimate of , we can obtain the estimates of ’s with varying degree of accuracy and computational complexity. The first and simplest way is to use one of as the estimate of . Instead of this estimate that poorly exploits the Toeplitz structure, we may well employ an averaging method to obtain the improved estimates defined by for Although these second estimates greatly improve the estimation accuracy over the first estimates ’s, they require more computational burden. Consequently, we also consider a partial averaging method to improve the accuracy with limited increase in computational complexity. To be more specific, we define the third estimates by if , and by if , with a positive integer suitably chosen depending on the situation at hand.
4.2. Choice of Parameters
Once we obtain the estimates of ’s, we can readily construct the orthonormal basis vectors of the signal space . However, we need to choose the right values for the parameters , , , , and to achieve optimal performance in a computationally efficient way. The DOA estimation in this paper will be carried out by equating with one of using Theorem 1. However, the exponential function exhibits periodicity with period Hence, we need to choose the values of and to satisfy to be able to identify uniquely the values of without any prior information about their values. To increase the array aperture size and hence to improve the resolution we need to choose the value of as large as possible. Consequently, we usually put and for best estimation accuracy.
Next, we discuss the choice of the parameters , , and for given value of We can achieve maximum utilization of the elements of if we put , , and . This is a reasonable choice for most situations since we can generally achieve better estimation accuracy with increased effective observation space dimension . However, we note that there is significant overlap between and , which can make the signal vectors ’s indistinguishable when the signal-to-noise ratio (SNR) is low. One possible way to reduce the overlap among the signal vectors ’s is to choose the signal vector separation index large. However, if we choose a large value for , it will decrease the effective observation space dimension , which can then reduce the resolution. Consequently, it is very important to choose an appropriate value for . We also note that the estimates of ’s are more accurate for small values and less accurate for large values if we adopt or for the estimation of Consequently, we may avoid using ’s with large values, which is particularly useful when the number of antenna elements is large. If we avoid using and , the maximum possible effective observation space dimension is given by .
4.3. DOA Estimation Algorithms: Spectral Search and Root Finding
Once we have chosen the parameters , , , , and appropriately as discussed in the previous subsection, we can apply Theorem 1 to determine the directions of arrivals. As usual, the directions of arrivals can be obtained by spectral search directly applying Theorem 1. This method is suitable when the values of ’s are already estimated and only additional tracking or finer estimation is necessary. However, the spectral search method may prove to be computationally expensive for initial estimation. To avoid the burden of spectral search, we can also adopt root finding method with slight modification on the usual method.
To discuss how the root finding scheme can be derived, let us note that the spectral search method searches for that satisfy From (24), we note that we can alternatively search for unit-norm zeros of the following equations:obtaining the root finding realization of the proposed algorithm. We note that the matrix is Hermitian and hence is also a zero of (25) if is a zero of (25).
5. Performance and Computational Complexity
5.1. Performance Evaluation
In this section, we compare the performance of the proposed algorithm with existing schemes. As the measure of performance, we use the root mean squared error (RMSE) defined by where denotes the estimate of . In what follows, the angles and RMSE will be measured in degrees rather than radians.
For performance evaluation, we first consider the situation in which equal strength uncorrelated narrowband signals arrive on a ULA with antenna elements separated by . We assume that the incident angles are , , , and . We assume that and that the second estimate is used to obtain unless stated otherwise. While we have developed the proposed algorithm valid for any integer value greater than or equal to , we will assume that is chosen to be in what follows to discuss various issues in a space-effective manner. We will also assume that is chosen to be in what follows. We first study the impact of the choice of the signal vector separation index on the performance. To study the trade-off of choosing different values, we consider , , and for the value of (We note that we cannot avoid using if we choose an odd number for ) For these values, the maximum possible effective observation space dimension is given by , , and , respectively. We have chosen these values in our simulations, the results of which are illustrated in Figure 2. For the simulations, we considered both spectrum search (SS) and root finding (RF) methods in implementing the proposed algorithm. In Figure 2, we clearly see that the different choice of value can result in significant difference in the performance. In particular, we note that the performance is significantly improved by increasing from to , although the performance is degraded if we further increase to Consequently, Figure 2 confirms the importance of reducing the overlap between signal space vectors by appropriately separating the elements. To see more intuitively how the choice of value affects the performance, we include two figures, namely, Figures 3 and 4, which depict the value defined by with . Figure 3 depicts the results with at three different signal-to-noise ratios, while Figure 4 represents those with To draw each of these figures, we randomly generated signals according to the specified average signal-to-noise ratio and the results are plotted as a function of the incident angle We clearly see that the peak values are much more easily recognizable with than with at the same signal-to-noise ratio. Before proceeding, we also note that the root finding method provides slightly lower RMSE values (except when ) in comparison with spectral search method, as in the case of MUSIC schemes [6, 18].
In obtaining the results given in Figure 2, we used the second estimate to obtain . However, the computation of can be burdensome in some situations. For this reason, we propose using in such situations with appropriately chosen values. To see how well we can approach, with reduced complexity, we consider the above system setup with The results of performance evaluation are given in Figure 5. We note that the second estimate reduces to the first and third estimates and , respectively, at the extreme cases of and In Figure 5, we observe that the performance is enhanced with increased values. We also note that the performance enhancement is particularly conspicuous when is increased from a small value. For example, by increasing from to or , we can achieve about dB or dB gain, respectively, at RMSE value of . However, we can achieve only about dB gain by increasing value from to Consequently, it is possible to reduce the computational complexity by employing the third estimate with reasonably small value without significant performance sacrifice.
Now, we compare the proposed algorithm with existing schemes, particularly with MUSIC-based schemes. The result is given in Figure 6, in which the proposed scheme is referred to as DSPCM standing for direct signal space construction method. In Figure 6, SS DSPCM and RF DSPCM represent the same curves labeled, respectively, as SS and RF in Figure 2. We observe that DSPCM holds good performance against the MUSIC schemes. In Figure 6, MUSIC refers to ordinary MUSIC scheme employing (19) to estimate the covariance matrix, while TMUSIC stands for the MUSIC scheme employing the second estimate to estimate the covariance matrix with Toeplitz property. By comparing SS MUSIC and SS TMUSIC, we observe that the enhanced accuracy with Toeplitz property provides better performance at low SNR. We note that the proposed scheme with is based on submatrix of the autocorrelation matrix, which is essentially a cross-correlation matrix between two ULAs with antenna elements with antenna spacing For this reason, we considered, to compare with the proposed scheme, two antenna configurations for the MUSIC algorithm, one consisting of antenna elements separated by distance and the other consisting of antenna elements separated by distance We note that RF provides better performance (compared with SS) and array configuration leads to better performance (compared with configuration). We have also evaluated the performance of the method proposed in [15]. Other subspace-based algorithms not requiring ED or SVD are compared in [14, 15]. According to [14], the algorithms proposed in [13] and [19] are outperformed by that in [14]. Also, the algorithm proposed in [15] is reported, in [15], to outperform that in [14]. Consequently, we can comfortably conclude that the algorithm proposed in this paper outperforms the algorithms in [13–15, 19], since it outperforms that in [15] by a significant margin as illustrated in Figure 6. Finally, we included the RMSE value obtained by Cramer-Rao bound [20]. The performances of the proposed algorithm and TMUSIC scheme are very close to the CRB for SNR greater than −8 dB. The performance of TMUSIC with even falls slightly below, at certain SNR values, with some small amount of bias induced in the simulations.
Figures 2–6 were obtained assuming that the incident angles are , , , and . To examine how the proposed scheme performs in more severe situations, we next consider the case in which the incident angles are , , , and . The results of performance evaluation are given in Figure 7. In this case, we clearly see the advantage of having more antenna elements in achieving better resolution. We also note that RF method which allows for zero slightly outside the unit circle is much better suited than the SS method when two or more incident angles are very closely spaced. Although the performance gap between RF DSPCM and RF MUSIC with becomes wider in Figure 7 than in Figure 6, we believe that the performance of RF DSPCM is still very competitive considering the relative computational simplicity. Unlike in Figure 6, the algorithm proposed in [15] has not been included since its performance cannot be included in the framework. For example, the RMSE values do not fall below even at SNR of dB. In Figure 7, we note two salient differences in comparison with Figure 6. First, the performance of the proposed scheme as well as that of TMUSIC is not very close to that of CRB. Secondly, the original MUSIC, while inferior to TMUSIC at low SNR, outperforms TMUSIC at very high SNR. From this, we note that the exploitation of Toeplitz property is particularly useful at low SNR.
To further study the resolvability of more closely spaced incident waves, we depicted values with three different scenarios in Figure 8. In each of the three scenarios, it is assumed that incident waves arrive in a ULA with antenna elements and that the value is chosen for the value of In the first scenario, the incoming waves are assumed to arrive at incident angles at and . We note that the peak values of are easily identifiable even at very low SNR of dB. In the next two scenarios, the incident waves are more closely spaced with angle separations and , respectively. We note that the peaks are again separable even though high signal strength is required in these situations.
Finally, we study the number of antenna elements required to resolve a given number of incident waves. We recall that the dimension of the reconstructed observation vector is given by . Consequently, to resolve distinct incident waves, we should have at least antenna elements satisfying In principle, we may choose values to satisfy . However, according to our evaluations, must be chosen to be no less than to resolve waves at reasonable SNR values. In Figure 9, we illustrate the simulation results with various and values. In each of the results in Figure 9, the value 2 is chosen for the value of and the incident angles of the waves are separated by For example, for , the incident angles are assumed to be , and . In each of the situations, is chosen to be such that except when From Figure 9, we see that the choice generally provides reasonable estimation performance except for The reason for relatively large required value for is related to the antenna aperture size. For example, even for , the choice that yields cannot provide enough antenna aperture size and hence cannot yield excellent performance. Consequently, even for , it is generally required to choose greater than
5.2. Computational Complexity
Three most computationally expensive steps in subspace-based DOA estimation schemes are the estimation of correlation matrix, the acquisition of signal or noise subspace, and the step of finding the DOA angles. In this subsection, we consider only the complexity of the second step, namely, the computational complexity of the signal or noise subspace acquisition step. There are several reasons for this. First of all, the first and third steps are required equally for other subspace-based schemes. Secondly, the complexity of these steps depends heavily on the situations at hand. For example, the first step of estimating the correlation matrix can be computationally burdensome for a single-shot or initial DOA estimation. However, adaptive differential correlation matrix update can be achieved with negligible or at least with reduced complexity in continuous target tracking applications. Similar argument holds true for the case of the third step. For example, the computational complexity of this step can greatly be reduced by adopting root finding method instead of spectral search. Moreover, in the case of continuous target tracking, extensive search is not strictly necessary since the target is likely to be in the neighborhood of the previous location. Consequently, we will focus solely on the computational complexity of the subspace acquisition step, which depends only on the algorithm we adopt.
First of all, we note that the proposed algorithm requires only about real floating point operations (rflops). As long as the number of targets and the effective observation space dimension are the same, the algorithm proposed here and the one in [15] require essentially the same computational complexity because similar signal subspace construction methods are used with only slight performance enhancing modifications. We next consider the complexity of subspace acquisition in MUSIC scheme. We recall that eigendecomposition of the correlation matrix is required in MUSIC scheme. The eigendecomposition of the correlation matrix, which is Hermitian, is performed in three steps, in general, to reduce the computational complexity. In the first step, the given Hermitian correlation matrix is converted into a Hermitian tridiagonal matrix. Next, the eigenvalues and eigenvectors of the tridiagonal matrix are obtained through algorithms such as divide-and-conquer or tridiagonal QR iteration method as discussed in [21, 22]. Then, the eigenvectors of the original Hermitian matrix are obtained by multiplying a matrix stored during the first step by the eigenvectors of the tridiagonal matrix. Since the computational complexity of the first step is more significant for our applications, we will consider only the complexity of the tridiagonal reduction process. As discussed in Sections 4.4.7 and 5.3.1 of [21], the first step of tridiagonal reduction is usually carried out with Hessenberg reduction process, which requires about arithmetics if both eigenvalues and eigenvectors of dimensional symmetric matrix are desired to be found. The same reduction process can be used for a complex Hermitian matrix but with complex arithmetics, approximately half of which are complex multiplications. We recall that a complex multiplication requires six real arithmetics while a complex addition corresponds to two real arithmetics. Consequently, we can conclude that the eigendecomposition of dimensional correlation matrix requires at least about rflops, which is significantly larger than rflops.
For numerical comparison of the computational complexity, let us consider the simulation setups used to produce the results in Figures 6 and 7. We recall that the number of sources was and that the effective observation space dimension was for the proposed scheme. Consequently, the proposed DSPCM requires 1,664 rflops. In contrast, the MUSIC schemes require at least about 10,667 and 85,333 rflops, respectively, for and element antenna array setups. Consequently, the MUSIC schemes in Figures 6 and 7 require at least about to times the computational complexity in constructing the signal or noise subspace. In fact, we note that the complexity difference will grow as a function of as the number of array antenna elements increases, which means more than a thousand times complexity reduction can be achieved when
6. Conclusion
In this paper, we proposed a low complexity subspace-based DOA estimation algorithm for a uniform linear array. In the proposed scheme, the signal space is constructed directly, without using computationally expensive ED or SVD process, by subsampling the autocorrelation matrix of a ULA. To avoid possible lack of distinguishability among the constructed signal space basis vectors, particularly at low SNRs, we proposed a signal vector separation method by slightly trading off the signal space dimension. Through simulations, we showed that the signal vector separation method can indeed improve the performance significantly. In addition to spectral search method, we provided root finding scheme as the angle finding methods. We showed that the performance of the proposed scheme is significantly better than the one in [15] almost matching that of the MUSIC schemes under various scenarios. Finally, we illustrated that the computational complexity of the proposed scheme is reduced by a factor of in comparison with MUSIC schemes, where is the number of sources and is the number of antenna arrays.
While the proposed algorithm exhibits a number of useful features, there are several issues to be further resolved for its practical use. First, we recall that the computational complexity of subspace-based DOA estimation schemes stems from three computationally burdensome stages, namely, the acquisition of correlation matrix, the construction of subspace, and the final angle finding stage. For this reason, we are planning to develop methods to further reduce the computational burden of correlation matrix acquisition and the final angle finding stage of spectral search or root finding. Next, we will study methods to avoid or reduce the detrimental effect of mutual coupling, since we are dealing with antenna elements spaced more densely than usual and hence the mutual coupling effect can become very critical.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2010-0025062 and NRF-2015R1D1A1A01060234).