2-D unitary ESPRIT-like direction-of-arrival (DOA) estimation for coherent signals with a uniform rectangular array.

A unitary transformation-based algorithm is proposed for two-dimensional (2-D) direction-of-arrival (DOA) estimation of coherent signals. The problem is solved by reorganizing the covariance matrix into a block Hankel one for decorrelation first and then reconstructing a new matrix to facilitate the unitary transformation. By multiplying unitary matrices, eigenvalue decomposition and singular value decomposition are both transformed into real-valued, so that the computational complexity can be reduced significantly. In addition, a fast and computationally attractive realization of the 2-D unitary transformation is given by making a Kronecker product of the 1-D matrices. Compared with the existing 2-D algorithms, our scheme is more efficient in computation and less restrictive on the array geometry. The processing of the received data matrix before unitary transformation combines the estimation of signal parameters via rotational invariance techniques (ESPRIT)-Like method and the forward-backward averaging, which can decorrelate the impinging signalsmore thoroughly. Simulation results and computational order analysis are presented to verify the validity and effectiveness of the proposed algorithm.

Introduction

Two-dimensional (2-D) direction-of-arrival (DOA) estimation of coherent signals has received much attention in many applications, such as radar, wireless communication and sonar in the multipath environment [1-5]. There are several high resolution techniques proposed to solve the rank deficiency of spatial covariance matrix caused by the presence of coherent signals. The conventional solution to this problem is the spatial smoothing method [6,7], which partitions the original array into a series of overlapping subarrays. Although it is efficient to decorrelate the incoming signals, peak searching of the spectrum in a 2-D space is required, which costs a large amount of computations. In order to reduce the computational complexity, an efficient method is performed by Hua [8]. This method, called the matrix enhancement and matrix pencil (MEMP) algorithm, exploits the structure inherent in an enhanced matrix from the original data. It estimates the azimuth and elevation separately in each dimension and combines them using a pairing method. However, the pairing result is not always correct when there are repeated parameters. Fortunately, a modified MEMP (MMEMP) method [9] is proposed to successfully solve the pairing problem.

In order to decorrelate the coherent signals thoroughly, recently, Han et al. [10] proposes an estimation of signal parameters via rotational invariance techniques (ESPRIT)-like algorithm for coherent DOA estimation. By reconstructing a Toeplitz matrix from the covariance matrix, this approach can decorrelate the impinging waves thoroughly. In [11], Chen extends it to the 2-D situation, namely the 2-D ESPRIT-like method, in conjunction with the MMEMP method, which outperforms the spatial smoothing method significantly in terms of the estimation accuracy. Although there is no peak searching existing in this algorithm, the computational burden is still heavy, due to the complex eigenvalue decomposition (EVD) and singular value decomposition (SVD) involved.

In this paper, we present a 2-D unitary ESPRIT-like (2-D UESPRIT-like) algorithm to reduce the computation complexity. Based on the block Hankel matrix obtained from [11], we preprocess it through a forward-backward average-like method convenient for unitary transformation. It can therefore transform the complex computations into real-valued ones and provide significant computational savings. The following DOA extractions are achieved simply by the one-dimensional (1-D) unitary ESPRIT [12], avoiding the computations of 2-D matrices. Simulation results will show that the real computations required for our new algorithm are much less than that of the 2-D ESPRIT-like method. It becomes especially obvious when the dimensionality of the Hankel matrix tends to be large. We also show that the variance of the estimates from our proposed method is close to the Cramer-Rao bound, and the resolution ability is superior to the others for the forward-backward average processing.

Background

Signal Model for URA

Consider K narrowband, far-field and coherent radiating sources with wavelength λ impinging on a URA of N × M identical and omnidirectional sensors with interelement spacing, d = d = λ/2. Using analytic signal representation, the received signal at the (n, m)th sensor can be expressed by: [7] where s(t) is the complex envelope of the kth wavefront, (β, γ) = (ejπ sincos, ejπ sinsin) and ϕ and θ are the elevation and azimuth angles of the kth source, respectively, and n(t) is the additive spatially white noise with variance . Figure 1 shows the sensor-source geometry configuration in the 2-D scenario. For simplicity, we define:

Figure 1.

Sensor-source geometry configuration for 2-D direction-of-arrival (DOA) estimation.

Therefore, (β, γ) can be expressed as (ej,ej). Rewriting Equation (1) in vector notation, we get: where (t) is the NM × 1 the data vector: (t) is the K × 1 signal vector: is the NM × K steering vector matrix: with as the NM × 1 steering vector of the kth source: and (t) is the NM × 1 noise vector: Here. we denote [·]T as the transpose and ⊗ as the Kronecker product.

From Equation (3), it follows that the covariance matrix of the received signal is given by: where E [·] denotes the expectation and (·)H represents the complex conjugate transpose.

According to the derivation in [11], the element of , for example, the cross correlation of the signals received at the (n,m)th and (p,q)th sensor for n,p = 0,…, N − 1 and m, q = 0,…, M − 1 is given by: where with (·)* being the conjugate and .

Real-Valued Processing for 1-D ULA

As the real-valued processing with a uniform linear array (ULA) provides the important preliminary knowledge to our new algorithm, we give a quick review of the definition of the unitary matrix and the real-valued processing based on it, which have been widely used in certain kinds of unitary transformation algorithms ([12,13]etc). Suppose there are only N sensors located on the x axis; left in Figure 1. N is odd, and the center of the ULA is the reference. The DOA of the incoming signal is denoted by (θ,ϕ = π/2). Then, the N × 1 steering vector of ULA can be written as: which is conjugate centro-symmetric. Such a property can be expressed mathematically as: where Π is the N × N exchange matrix with ones on its antidiagonal and zeros elsewhere: Define as the 1-D unitary matrix. By multiplying it, the elements in )(θ) can be transformed into real quantities, such that: If N is even, a similar result can be obtained using:

Proposed Algorithm

Signal Decorrelation

The proposed method is developed in the 2-D scenario, which has been introduced in Section 2.1. In order to resolve the rank deficiency problem caused by signal coherency, we first construct the following Hankel matrix from Equation (5) [11]: Then arranging a series of the above Hankel matrices into a block Hankel one, we have: The analytic expression of (n, m) is given by [11]: where = [1,…,] with , (n,m) = diag{d1,,,…,d} and with . It has been proven that (n, m) is full-rank, provided the values of P and Q are selected properly. In this case, the rank of (n, m) equals the number of incoming signals. Therefore, lots of high resolution methods for 2-D DOA estimation of the uncorrelated or partly correlated signals can be used.

It should be noted that the following derivation will be performed under the assumption that there is no noise existing in the received data, which can be seen from Equation (11). Further study on the complex situation with spatially white noise will be carried out in Section 5 through several simulations.

Real- Valued Processing

Although we can apply the eigenstructure techniques to estimate 2-D DOA based on the full-rank (n,m), the computational burden is much heavier, because of the complex computations involved in it. In this note, we develop a 2-D unitary transformation method to reduce the complex computations to real ones.

If we premultiply and postmultiply (n, m) with unitary matrices directly, (n, m) cannot be transformed into real-valued, because the matrix (n, m) is complex. Therefore, we need to construct a new matrix associated with (n, m) before unitary transformation to guarantee this property. Let us define: Then we have: With the property of the Kronecker product α1α2(1⊗2) = α11⊗α22, where α1 and α2 represent any constant and 1 and 2 represent matrices, can thus be converted into: Furthermore, we have where = diag{ej [(−1)1+(−1)1]/2,…, ej [(−1)+(−1)]/2}, with , (u) = [e−j(−1)/2,…,ej(−1)/2]T and (υ) = [e−j(−1)/2,…,ej(−1)/2]T. (u) and (υ) are both conjugate centro-symmetric, the same as the steering vector given by Equation (6).

Substituting Equations (14) and (15) into Equation (13) can be rewritten as: where + * is real, with = (n,m)T*H(n,m)H. Compared with Equation (4)′ in Equation (16) can be viewed as a new array response matrix and + * as the equivalent covariance matrix of the incoming signals. Notice that the achievement of is consistent with the forward-backward average processing [14,15]. Compared with (n,m), decorrelates the coherent signals more thoroughly and, therefore, can provide higher estimation accuracy.

Since (u) and (υ) in ′ have the similar form as Equation (6), we can make use of the 1-D real-valued processing-like Equation (8) to perform the 2-D unitary transformation. Define the 2-D unitary matrix as: where and use the definition of Equation (7) or (9). Premultiplying and postmultiplying by , yields: According to the useful property: φ becomes: with . As the vector in each column of , for example , satisfies Equation (8), the matrix is real. Combined with the real-valued + *, we can easily deduced that ϕ is real.

Extracting βk and γk

To avoid peak searching, to retain the 2-D DOA estimation real-valued and reduce the computational complexity, we develop a simple implement of the 2-D unitary ESPRIT based on the 1-D solution [12].

Let eigenvectors associated with the K largest eigenvalues of be denoted by ∈ C×. Since and ′ both span the signal subspace of , as Equation (16) has shown, there exists a unique, nonsingular matrix, ∈ C×, such that = ′. Define 1 = 1 and 2 = 2 as the first and last, (P − 1)Q, rows of matrix, , respectively, with 1 = [(−1), 0(−1)×] and 2 = [0(−1)Q×Q, (−1)]. Then, replacing by ′, we have: Observing the structure of ′, we can find that: with Ω = diag{eju1,…, ej}. Combined with Equations (19) and (20), we can write the relationship between 1 and 2 as: It is clear that Ψ is the total least squares (TLS) solution of Equation (21). The eigenvalues of Ψ are equal to the diagonal elements of Ω. Therefore, u corresponding to β can be extracted by an EVD of Ψ. In [12], the author has proved that the TLS problem can be solved through a SVD of [1, 2]. Here, we define: By reconstructing the matrix, [1, 2], and performing the unitary transformation, the TLS problem can be solved by computing the SVD of the real matrix: with 2 defined by Equation (9). Note the eigenvectors are associated with the K largest eigenvalues of matrix, ϕ, as . From Equation (17), we have: Substituting Equations (22) and (24) into Equation (23) and rearranging it, yields: where and Since the computations of 1, 2 require a large amount of multiplies, we intend to simplify their expressions with lower dimensions. Let which are selection matrices used in the 1-D case [12]. It follows that and . Using the property Equation (18)1 can be simply determined by: where . As [12] stated, the matrix is real. Similarly, we have: with . Therefore, the computations of 1 in Equation (25) can be greatly simplified by using Equations (26) and (27).

Denote the right singular vector matrix of 1 by and partition it into submatrices: Finally, we get the real-valued TLS solution corresponding to Ψ as: To extract the parameter υ, we define different P(Q − 1) rows of as 3 = 3 and 4 = 4s, in which and with: Similar to the derivation of Equation (21), we have: where Ω = diag{ej1, …,ej} and are the nonsigular matrix in Equation (19). Ψ is the TLS solution, which can be obtained by the SVD of real matrix: where with: Partitioning the right singular vector matrix of the real P(Q − 1) × 2K matrix 2 into a block one, we obtain: Then, it follows that the real-valued TLS solution associated with Ω can be computed by:

Employing the existing automatic pairing method [16,17], the estimates, υ and υ, can be achieved by computing the EVD of the “complexified” matrix ϒ + jϒ. Denote the real and the imaginary parts of the K eigenvalues as and . The estimation of u and υ will be: From Equation (2), we can get the final result:

Summary of the Algorithm

The steps of the proposed method are described as follows:

Obtain = [x(t1),…, x(t)] with x(t) as the snapshot at time, t. Then, compute the covariance matrix approximately by = with L as the snapshot number.

Construct the block Hankel matrix (n, m) by Equation (10) to decorrelate the coherency of signals and obtain through Equations (12) and (13) to facilitate the following unitary transformation. Then, compute the real-valued matrix .

Compute as the K dominant eigenvectors of ϕ and calculate 1, 2 through Equations (25) and (29). Conduct the SVD of 1, 2 to obtain the right singular vector matrices, , , and get ϒu, ϒ from Equations (28) and (30), respectively.

Conduct EVD of the complex-valued matrix, ϒu + jϒ. Extract u and υ from the eigenvalues by Equation (31). At last, estimate θ and ϕk using Equation (32).

Computational Order Analysis

In the following, we will first derive an estimate of the order of real multiplications involved in each step. Then, we will compare the computational order of our new method, namely the 2-D UESPRIT-like method, against that of the 2-D ESPRIT-like method.

Computational Order of Step 1

Here, we calculate by H/L with ∈ C×. The direct computation requires real multiplications.

Computational Order of Step 2

From Equation (13), we can see that to obtain , only the computation of = (n, m)H (n, m) is needed. It is because can be achieved by rearranging the elements of y simply, without any multiplication. According to the computational analysis in [8], the minimum number of real multiplications required to compute is: provided P ≫ 1,Q Gt; 1,N − P Gt; 1,M − Q Gt; 1.

In this step, we also need to calculate ϕ. Due to the special structure of unitary matrices, U and U, the computation of = ⊗ only contains that of the product of j and U. Therefore, the order of computing is 2Q. For the same reason, the real multiplications involved in is:

Computational Order of Step 3

As φ is a real-valued matrix, the number of multiplications required (based on the symmetric QR algorithm [18]) for its EVD is: when PQ ≫ 1.

The following real multiplications involved in the ϒ and ϒ achievement is listed in Table 1. Denote the total number involved in it as C1. The computational order of SVD is obtained by the Chan SVD [18].

Table 1.

Real multiplications involved in the computations of ϒ and ϒ.

	Real multiplications	conditions
K1′,K2′	2 [(P − 1)2P + P2(P − 1)]	/
K1,K2	2PQ(P − 1)	/
K1Esϕ, K2Esϕ	2PQ2(P − 1)L	/
SVD of T1	(2L)2(P − 1)Q+ 17(2L)3/3	PQ > (P −1)Q
ϒu	L3	/
K3′,K4′	2 [(Q − 1)2Q + Q2(Q − 1)]	/
K3,K4	2PQ(Q − 1)	/
K3Esϕ, K4Esϕ	2P2Q(Q − 1)L	/
SVD of T2	(2L)2P(Q − 1)+ 17(2L)3/3	PQ > P(Q − 1)
ϒυ	L3	/

Computational Order of Step 4

In this step, only the EVD of complex-valued matrix, ϒ + jϒ, is considered. It requires 4 × 5L3 real multiplications when L≫1.

Comparison to the 2-D ESPRIT-Like Method

According to the above analysis, the order of computations needed by the 2-D UESPRIT-like method is: as long asP > L ≫ 1,Q > L ≫ 1,N ≫ 1,M ≫ 1. It can be seen that the computation of and the EVD of ϕ occupy the major part.

In contrast, we present the computations needed in the 2-D ESPRIT-like method [11], which uses the same decorrelation processing, but a different MMEMP method behind it. As stated in [8], the multiplications required in each step are listed in Table 2. C2 = 3PQL2 + 7L3, the number of complex multiplications used in MMEMP, is obtained under the assumption that there are no repeated β. From Table 2, the real computations of the 2-D ESPRIT-like method can be written as:

Table 2.

Real multiplications required for the 2-D UESPRIT-like and 2-D ESPRIT-like method.

	2D UESPRIT-Like Method	2D ESPRIT-Like Method
Rx	2M N(M N +1)L	2M N(M N +1)L
Ry	4 × [2PQ(N − P/2)(M − Q/2)]	4 × [2PQ(N − P/2)(M − Q/2)]
φ	2(PQ)2 + 4PQ	/
EVD	of φ: 5(PQ)3	of Ry: 20(PQ)3
The rest of the operations	2D Unitary ESPRIT:C1 + 20L3	MMEMP: C2

According to Equations (34) and (35), the computational saving of our method is about 15(PQ)3, caused by the unitary transformation. In [8], the author recommended the scope of P, Q given below: in which the estimation accuracy and computations become both higher as the values increase. For comparison, we estimate the most multiplications cost in computing (n, m) by choosing P = (N + 1)/2 and Q = (M + 1)/2. In this case, the first part of Equations (34) and (35) becomes: In the condition of P > L ≫ 1 and Q > L ≫ 1, Equation (33) is far more than Equation (36), which indicates that the computational saving of our algorithm is considerable.

Figure 2 is plotted to compare the proposed method with the 2-D ESPRIT-like method in the aspect of real computations required as a function of P and Q. The multiplications cost in each method are computed by the sum of the corresponding column in Table 2. The snapshot number is L = 3, and the size of the URA is N × M = 30 × 20. P and Q change in the scope of [L + 1, (N + 1)/2] and [L + 1, (M + 1)/2]. The figure shows that as P and Q increase, the complexities of estimating the 2-D DOA with two different algorithms increase, as well. The computations needed for our new method is much less than that of the 2-D ESPRIT-like method when P and Q go towards the upper bound. Therefore, we conclude that the proposed scheme can obtain dramatic computational savings in estimating the elevation and azimuth.

Figure 2.

Real computations needed to estimate the 2-D DOA as a function of P and Q.

Simulation Results

In this section, we present simulation results that compare the proposed method with several other 2-D DOA estimations in the presence of a zero mean Gaussian white noise. Except the developed scheme, DOA extractions and pairings in other methods are all performed using MMEMP [9].

Suppose K narrowband equipowered signals are incident on a 11 × 9 URA (i.e., N = 11, M = 9) with interelement spacing d = d = λ/2. The correlation factor between the first signal and the others is denoted as γ. When γ equals 1, they are coherent and completely uncorrelated when γ = 0. We generate correlated signals via: where s(t) is the amplitude of the kth signal at time, t, SNR denotes the signal to noise ratio and r(t) is the output of the kth random number generator at time, t. In the range of Equation (36), we select P = 5 and Q = 4.

First, we consider three coherent signals from (ϕ, θ) = (10°, 0°), (5°,−5°), (−5°, 7°) with SNR varying from −20 dB to 10 dB. 1,000 Monte Carlo trials are run for each SNR. Five hundred snapshots are taken to form the estimate of the covariance matrix of the array outputs. Figure 3 shows the probability of identifying all the directions correctly versus SNR. The result illustrates that the performance of both the proposed algorithm and 2-D ESPRIT-like method are better than that of 2-D spatial smoothing (2-D SS). This is because the first two algorithms can eliminate the coherency of signals completely by reconstructing the equivalent covariance matrix (n, m) by Equation 10, while the 2-D SS method can only provide the alleviation of the rank deficiency to some extent. The graph also shows that the new method has some improvement over the 2-D ESPRIT-like method, due to the fact that the achievement of in our method is similar to the forward-backward spatial smoothing process [14], which can further enhance the ability of decorrelation.

Figure 3.

Probability of the DOA identification versus input signal to noise ratio (SNR).

Second, we consider the same scenario as in the first one. Define the root mean square error (RMSE) of the DOA estimates as: where and are the estimates of θ and ϕ for the nth trial,respectively. K is the number of signals. The comparison of Cramer-Rao lower bound (CRB), computed according to formulas provided in [19] and the RMSE of DOA estimates of 2-D UESPRIT-like method, 2-D ESPRIT-like method and 2-D SS are plotted in Figure 4. The SNR changes from −10 dB to 20 dB. Simulation result shows that the estimation errors of all the methods decrease obviously as the SNR increase. Moreover, our proposed method is observed to have a superior performance over the others and to be close to the CRB most. When the SNR is lower than 0 dB, the 2-D SS method fails to distinguish the two closely located sources, while our algorithm can still accomplish it very well. The same phenomenon appears for the 2-D ESPRIT-like method when the SNR is lower than −7 dB.

Figure 4.

Root mean square error (RMSE) of the DOA estimates versus input SNR.

The third simulation considers the same signals as the above two experiments with SNR = 20 dB. The snapshot number, L, changes from 10 to 300. The RMSE of the 2-D estimates as a function of the snapshot number and the CRB are plotted in Figure 5. As expected, the 2-D UESPRIT-like method is shown to have dramatically high accuracy over the other two and can achieve the closest performance to CRB.

Figure 5.

RMSE of the DOA estimates versus snapshot number.

In the last simulation, we investigate the ability of four algorithms to decorrelate coherent signals. Assume there are two narrowband signals with (ϕ, θ) = (10°, 0°), (5°, −5°) and SNR = 5dB. Their correlation factor varies in the range (0,1). For each value of γ, 1,000 independent estimates are carried out to examine the RMSE of DOA estimates. First, we compare the 2-D ESPRIT method with the other three decorrelation algorithms. As Figure 6 has shown, when the signals are uncorrelated or partly correlated with low correlation factor, the conventional 2-D ESPRIT is the best choice for DOA estimation. The reason is the use of decorrelation in the other methods results in a small effective aperture, which can reduce the resolution significantly. However, as γ increases, the performance of the 2-D ESPRIT degrades slowly. Until γ = 0.9, it totally fails to identify DOA, because of the serious rank loss of . Besides, Figure 6 also demonstrates that the 2-D SS method provides a much better performance than the 2-D ESPRIT-like method, as well as our proposed method, when r ≤ 0.32. This is due to the fact that in such a low correlation, it is sufficient for 2-D SS to decorrelate signals, and it can restrain noise to some extent. When the signals become nearly coherent, that is γ ≥ 0.95, the superiority of our proposed scheme and the 2-D ESPRIT-like method appears to be remarkable for the excellent decorrelation ability, while in terms of the veracity of DOA estimates, it is obvious that the new method outperforms 2-D ESPRIT-like method all the time.

Figure 6.

RMSE of the DOA estimates versus correlation factor.

Conclusions

An application of unitary transformation to 2-D DOA estimation of coherent sources has been proposed in this paper. The decorrelation is performed based on the existing 2-D ESPRIT-like method. While the computational load is significantly reduced by transforming the complex matrices into real-valued ones and conducting the EVD with 1-D matrices, the 2-D UESPRIT-like method can also provide better performance in DOA estimation by preprocessing the block Hankel matrix using the forward-backward averaging-like method. Computational analysis and simulation results have shown the significant reduction of computations and the dramatic low RMSE in DOA estimations. A less restrictive requirement for the array geometry is also provided to generalize the application of this method.

It is worth mentioning that our proposed algorithm is fit for dealing with the estimation of highly correlated signals. In the situation where all the signals are uncorrelated or partly correlated, the method given in this paper will suffer degradation to some extent.

It is also interesting to notice that our new algorithm is still practically useful in the presence of mutual coupling, though such an effect was not considered in this paper. For 2-D DOA estimation in URA, if the mutual coupling matrix (MCM) is known, the coupling effect can be easily eliminated by premultiplying the inverse MCM with the original data. If the MCM is unknown, we can also use the method provided in [20] to completely alleviate the mutual coupling by setting the sensors on the boundary of the URA as auxiliary sensors, provided that the MCM satisfies the block banded symmetric Toeplitz assumption. The output of the middle URA is, therefore, an ideal model without a coupling effect. Moreover, the 1-D version of our proposed method, namely 1-D UESPRIT [21], can be extended to the azimuth estimation with the Uniform Circular Array (UCA) in the presence of mutual coupling. In this case, the elevation is assumed to be known. The mutual coupling can be directly compensated for if the coupling effect is not so serious and the MCM is known [22]. Then, the array manifold of the UCA is projected by a coupling matrix onto that of an ideal UCA, where the azimuth estimation is the same as that of ULA, and the 1-D UESPRIT method can be directly conducted. However, it cannot be applied to the method proposed in [23] when the number of antenna elements in the UCA is large enough, as the H matrix used to incorporate all the relevant phase modes into the principal term destroys the Vandermode structure of the beamspace steering vector. Besides, the new method may not be applied as well in the more realistic situation of 2-D DOA estimation with UCA, as the UCA ESPRIT and UCA-ESPRIT-like involved in the algorithms [19,24] are different from the 2-D ESPRIT of URA [17]. Moreover, the elevation-dependence of MCM prevents the application of our approach to the method in [25]. Future work may focus on the real-valued processing of UCA ESPRIT, utilizing the special structure of unitary matrix.