554
IEEE COMMUNICATIONS LETTERS, VOL. 22, NO. 3, MARCH 2018
Unitary Direction of Arrival Estimation Based on A Second Forward/Backward Averaging Technique Feng-Gang Yan , Member, IEEE, Shuai Liu, Jun Wang , Ming Jin, Member, IEEE, and Yi Shen, Member, IEEE
Abstract— A second forward/backward (SFB) averaging technique is presented to transform the real symmetrical covariance matrix of the popular unitary root-MUSIC (U-root-MUSIC) into a real bisymmetrical one, based on which a SFB-Uroot-MUSIC direction of arrival estimation algorithm with reduced computational complexity is developed. The proposed SFB-U-root-MUSIC algorithm reduces the computational complexity in the eigenvalue decomposition (EVD) stage by a factor about four because it performs real-valued EVDs on two submatrices of about half sizes. The proposed SFB averaging technique is further extended as a generalized dimension reduction method to other unitary DOA estimators for low-complexity EVD computation, and numerical simulations are conducted to demonstrate, such a dimension reduction technique sacrifices statistically nonsignificant root mean square performance that is acceptable. Index Terms— Direction of arrival (DOA) estimation, second forward/backward (SFB) averaging, real-valued computation, unitary root-music (u-root-music).
I. I NTRODUCTION
E
STIMATING the direction of arrivals (DOAs) of multiple narrow-band signals is a fundamental problem that arises in various engineering applications such as radar, sonar, navigation and wireless communication [1]– [5]. This topic has been addressed extensively over the past several decades, and many algorithms have been proposed regularly. As one of the most critical techniques in the field, reducing the computational complexity of various DOA estimators is of great interest in the literature because complexity is generally an important index for many practical systems [6], [7]. Unitary transformation [8] is one of the most representative methods to reduce the complexity, which achieves high computational efficiency via real-valued computations [9]. Compared with conventional complex-valued algorithms such as MUSIC [10], ESPRIT [11], root-MUSIC [12] and their derivations [13]– [15], unitary estimators are able to reduce about 75% computational complexity [16]. Besides, unitary
Manuscript received October 13, 2017; revised November 3, 2017; accepted November 3, 2017. Date of publication November 8, 2017; date of current version March 8, 2018. This work is supported by National Natural Science Foundation of China (61501142), Science and Technology Program of WeiHai, Project Supported by Discipline Construction Guiding Foundation in Harbin Institute of Technology (Weihai) (WH20160107), and the Fundamental Research Funds for the Central Universities (HIT.NSRIF.201725). The associate editor coordinating the review of this letter and approving it for publication was F. Gao. (Corresponding author: Jun Wang.) F.-G. Yan, S. Liu, J. Wang, and M. Jin are with the School of Information and Electrical Engineering, Harbin Institute of Technology at Weihai, Weihai 264209, China (e-mail:
[email protected]). Y. Shen is with the School of Astronautics, Harbin Institute of Technology, Harbin 150001, China. Digital Object Identifier 10.1109/LCOMM.2017.2771489
methods also show improved accuracy and allow reduced complexity with enhanced accuracy [17]. However, it is worth noting that almost all of state-of-the-art unitary methods still involve an eigenvalue decomposition (EVD) step on a covariance matrix of dimensions M × M or on a transformed covariance matrix of the same sizes, where M stands for the number of sensors. Generally, this EVD computation requires about O M 3 flops [18]. When massive arrays are used [19], [20], M can be a very large number and this term of high complexity is unacceptable. Therefore, there is a need for techniques demanding less EVD computations. In this letter, we propose a reduced-complexity formulation of the popular U-root-MUSIC algorithm [21]. We exploit the forward/backward (FB) averaging technique [22], [23] again on the real symmetrical covariance matrix of U-root-MUSIC to obtain a second FB (SFB) averaging covariance matrix, which is found symmetrical about both of its diagonals. Such a bisymmetrical characteristic allows fast EVD on two submatrices of about half sizes, and hence the proposed SFB-Uroot-MUSIC can further reduce the complexity in the EVD stage by a factor about four. We also show that the proposed SFB averaging technique can be regarded as a generalized dimension reduction method for other unitary DOA estimators for the purpose of low-complexity EVD computation. Mathematical Notations: Matrices and vectors are denoted √ by upper- and lower- boldface letters, respectively, j −1, (·)T is transpose, (·) H is Hermitian transpose and E [·] is mathematical expectation. In addition, 0, I and J stand for the zero-, the identity- and the exchange- matrices, respectively. II. S IGNAL M ODEL AND U-ROOT-MUSIC Assume L uncorrelated narrow-band signals sl (t), l ∈ [1, L] impinge from far-field at directions [θ1 , θ2 , . . . , θ L ] on a uniform linear array (ULA) composed of M antenna sensors with inter-sensor spacing d. It is assumed d μ/2 to avoid phase ambiguity, where μ is the wavelength. It is also assumed M > 2L which is reasonable for large ULAs [19], [20]. Let the first sensor be the reference point, the array output at snapshot t, t ∈ [1, T ] can be expressed as [8]–[24] x (t) = As(t) + n (t) , (1) where A = a(θ1 ), a(θ2 ), . . . , a(θ L ) is the M × L steering vector matrix, and (M−1) T (2) a (θl ) p(z l ) = 1, z l , · · · , z l
is the M × 1 steering vector, where z l j (2π/μ)d sin θl ,l ∈ [1, L], s(t) is the M × 1 signal vector, n(t) ∼ CN 0, σn2 I is
1558-2558 © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
YAN et al.: UNITARY DIRECTION OF ARRIVAL ESTIMATION BASED ON A SFB AVERAGING TECHNIQUE
the M ×1 additive noise vector with σn2 denoting the the noise power. The M × M covariance matrix is given by [8]–[24] R = E x (t) x H (t) = ASA H + σn2 I, (3) where S E s (t) s H (t) is the L × L signal covariance matrix. The covariance matrix is estimated using T snapshots of observed data as T 1 x (t) x H (t) . R= T
(4)
t =1
Because R is complex, traditional methods based on the EVD of R require complex-valued computations accordingly. To realize real-valued computations, unitary algorithms [8], [17], [21] exploits the FB averaging covariance matrix [22] 1 R + J R∗ J = AS1 A H + σn2 I, (5) RFB = 2 to obtain a symmetrical real matrix instead of R RUnitary Q H RFB Q = Re Q H RQ , (6)
−(M−1) ,··· , where S1 = 12 S + DS∗ D H , D = diag z 1 −(M−1) [21], and Q is a unitary matrix satisfying Q−1 = zL Q H , which is determined by the parity of M in two cases as
1 I jI Q= √ , M = 2k (7) 2 J − jJ ⎤ ⎡ I √0 jI 1 (8) Q = √ ⎣0 2 0 ⎦ , M = 2k + 1. 2 J 0 − jJ Unitary as Performing real-valued EVD on R S N UTS + UTN , US UN RUnitary =
(9)
where the subscripts S and N stand for the signal- and noisesubspaces respectively, one can exploit the real-valued noise matrix U N for further DOA estimate. For example, we can define the following U-root-MUSIC polynomial UTN p˜ z , (10) f U−root−MUSIC (z) p˜ T z −1 UN where p˜ (z) Q H p (z), and estimate DOAs by finding the L roots of (9) which lie closest to the unit circle as [12], [21] μ (11) z l , l ∈ [1, L]. θl = arcsin 2πd III. T HE P ROPOSED A LGORITHM Let us exploit the FB averaging technique again on RUnitary to obatin the following SFB averaging covariance matrix 1 RUnitary + J RSFB RUnitaryJ . (12) 2 Unitary = A1 S1 AT + σn2 I, where Inserting (5) into (6) gives R 1 H A1 Q A. Therefore, we have 1 RSFB = A1 S1 A1T + JA1 S1 A1T J + σn2 I 2
S1 0 1 A1T = A1 JA1 + σn2 I 0 S1 (JA1 )T 2 = A2 S2 A2T + σn2 I, (13)
555
where S2 diag {S1 , S1 }, A2 [A1 JA1] can be regarded as a virtual steering vector matrix covering A1 and JA1 , and we have span (A2 ) = span (A1 ) ⊕ span (JA1). Performing EVD on RSFB as S N RSFB = VTS + VTN , V VT = VS VN (14) we must obtain span (A2 ) = span V S and span (A2 ) ⊥ span V N . Since span (A1 ) ⊆ span (A2 ), it is clear that span (A1 ) ⊥span (15) VN . Using equation (15), we obtain the following polynomial f SFB−U−root−MUSIC (z) p˜ T z −1 VN VTN p˜ z . (16) By inserting the L roots of (16) that lie closest to the unit circle into (11), DOAs can be estimated immediately. T Notice that RUnitary = RUnitary and J2 = I, we have T RSFB = RSFB , J RSFB J = RSFB .
(17)
Hence, RUnitary is a bisymmetric matrix [18]. Based on such a mathematical fact, we now show that instead of using (14), the noise matrix V N in (16) can be computed by EVDs of two sub-matrices of about half sizes, which is investigated with respect to the parity of M in two cases as follows. 1). M = 2k. According to (17), we can divide RSFB into
R J R21J RSFB = 11 , M = 2k. (18) R21 J R11J Define two k × k sub-matrices R1 R11 − J R21 R2 R11 + J R21.
(19) (20)
Clearly, R1 and R2 two symmetrical matrices, and their EVDs 1 S,1 N,1 R1 = R1T = VTS,1 + VTN,1 (21) V1 V S,1 V N,1 2 = 2 S,2 N,2 V2T = VTS,2 + VTN,2 (22) V2 V S,2 V N,2 R k
require only real-valued computations [18]. Assume λ1,i i=1 k
and λ2,i i=1 are the k eigenvalues of R1 and those of R2 , respectively, we have
λ1,1 , λ1,2 , · · · , λ1,k (23) R1 V1 = diag V1T
T 2 V 2 = diag 2 R (24) λ2,1 , λ2,2 , · · · , λ2,k . V Introduce the following 2k × 2k matrices
V2 V1 Z1 √1 , 2 −J V1 J V2 and another two 2k × 2k matrices
1 I I V1 X1 √ , Y1 −J J 0 2
(25)
0 . V2
Because Z1 = X1 Y1 , Y1−1 = YT and X1−1 = XT , we = Y1−1 X1−1 = Y1T X1T = Z1T . Therefore, Z1 Z−1 1 unitary matrix [18]. Using (19), (20), (23) and (24), straightforward to show that
Z1T RSFB Z1 = diag λ1,k , λ2,1 , · · · , λ2,k . λ1,1 , · · · ,
(26) have is a it is (27)
556
Therefore, the 2k eigenvalues of RSFB are composed of
k k λ1,i i=1 and λ2,i i=1 such that
2 , M = 2k. = diag 1 , (28)
IEEE COMMUNICATIONS LETTERS, VOL. 22, NO. 3, MARCH 2018
TABLE I D ETAILED S TEPS FOR I MPLEMENTING THE P ROPOSED A LGORITHM
In addition, the columns of Z1 are the 2k eigenvectors of RSFB , which are jointed by the k eigenvectors of R1 and those of V N can be given by R2 . Consequently,
N,2 N,1 V V , M = 2k. (29) VN = −J V N,1 J V N,2 2). M = 2k + 1. Based on (17), we can divide RSFB in this case into ⎤ ⎡ R11 d J R21 J (30) RSFB = ⎣ dT a dT J ⎦ , M = 2k + 1. R21 J d J R11 J To simplify the notations, we reuse the definitions for R1 and R2 in (19) and (20). In addition, we further define another (k + 1) × (k + 1) matrix √ T a 2d . (31) R3 √ 2d R2 Since R3 is also a real symmetrical matrix, its EVD 3 = 3 S,3 N,3 R R3T = VTS,3 + VTN,3 V3 V S,3 V N,3 (32)
k+1 requires only real-valued computations. Assume λ3,i i=1 are 3 , we have the k + 1 eigenvalues of R
V3T (33) R3 V3 = diag λ3,1 , λ3,2 , · · · , λ3,k+1 . Introduce the following (2k + 1) × (2k + 1) real matrices ⎡ ⎤ V√ V1 3 (2 : k + 1) 1 ⎣ Z2 √ (34) 0 2· V3 (1) ⎦ , 2 −J V3 (2 : k + 1) V1 J where V3 (1) is an 1 × (k + 1) vector, denoting the first row of V3 , V3 (2 : k) is a k × (k + 1) matrix composed of the last k rows of V3 . Define another two (k + 1) × (k + 1) matrices ⎤ ⎡
I I √0 1 V1 0 ⎦ ⎣ X2 √ . (35) 0 2 0 , Y2 0 V3 2 −J 0 J
V3 (1) By dividing V3 as V3 = , it can be easily V3 (2 : k + 1) −1 −1 −1 verified that Z2 = X2 Y2 and Z2 = Y2 X2 = Y2T X2T = Z2T . Therefore, Z2 is a unitary matrix [18]. Using (23), (30), (33) and (34), we can similarly obtain
Z2T RSFB Z2 = diag λ1,1 , · · · , λ1,k , λ3,1 , · · · , λ3,k+1 . (36) eigenvalues of RSFB are composed of
Hence, k the 2k + 1 k+1 λ1,i i=1 and λ3,i i=1 such that
3 , M = 2k + 1. = diag 1 , (37) In addition, the columns of Z2 are the 2k + 1 eigenvectors of SFB , which can be jointed by the k eigenvectors of R R1 and V N can be given by the k + 1 those of R3 , and ⎡ ⎤ V N,1 V√ N,3 (2 : k + 1) VN = ⎣ 0 2· V N,3 (1) ⎦ , M = 2k + 1. (38) V N,3 (2 : k + 1) −J V N,1 J
Fig. 1. RMSE versus the SNR respect to the source at θ1 = 30◦ , M = 10 sensors, T = 100 snapshots, L = 2 sources at θ1 = 30◦ and θ2 = 40◦ .
In summary, detailed steps for implementing the proposed algorithm fast DOA estimation are given in Table I. Remark 1: Because R1 , R2 and R3 are all of about half sizes as compared to RUnitary , the proposed method provides an obvious complexity reduction in the EVD stage by a factor about four as compared to U-root-MUSIC. This is achieved by the price of enhanced caching requirements. Because the proposed method needs to store the M × M SFB matrix and two of the M/2 × M/2 sub-matrices, about 3/2M 2 additional memories are needed by the proposed method. Remark 2: Although we elaborated the SFB averaging technique with U-root-MUSIC, it is clear that by using the EVDs 2 and R3 instead of that of two sub-matrices among R1 , R of RUnitary, this SFB averaging technique can be similarly extend to other unitary DOA estimators for the purpose of low dimensional EVD computations. IV. S IMULATION R ESULTS Numerical simulations with 500 independent Monte Carlo trials are conducted to assess the performance of the SFB-U-root-MUSIC algorithm and compare it with U-rootMUSIC [21]. Since the proposed SFB averaging technique can be directly extended to reduce the complexity of other unitary DOA estimators including U-MUSIC [9] and U-ESPRTI [17], the modifications based on the proposed SFB averaging technique of those methods, namely, SFB-U-ESPRIT and SFB-U-MUSIC, are also considered. For the RMSE performance comparison, the unconditional Cramér-Rao Lower Bound (CRLB) [24] is applied for a common reference. In the first simulation, we compare the RMSE performances of different algorithms. We plot the RMSEs as functions of the
YAN et al.: UNITARY DIRECTION OF ARRIVAL ESTIMATION BASED ON A SFB AVERAGING TECHNIQUE
557
R EFERENCES
Fig. 2. RMSE versus the number of snapshots respect to the source at θ1 = 30◦ , SNR = 0dB. M = 10 sensors, T = 100 snapshots, L = 2 sources at θ1 = 30◦ and θ2 = 40◦ . TABLE II C OMPARISON OF CPU T IME IN S ECOND
SNR in Fig. 1 and plot those as functions of the number of snapshots in Fig. 2. The simulation parameters are given in the figure captions. It is seen clearly from the two figures that the SFB-based techniques provide similar performances very close to their original versions all along different SNRs and T ’s. Considering that the proposed techniques requires only low-dimension EVD computations with reduced complexity, the new technique makes a sufficiently efficient trade-off between complexity and accuracy. Next, we compare the computational efficiency in terms of CPU times in Table. II. The simulation parameters are set as SNR = 0dB, T = 100, L = 2, θ1 = 30◦ and θ2 = 35◦. The CPU times are given by running the Matlab codes on a PC with the same configurations in the same environment. It is seen that the SFB-based techniques cost much lower CPU time than their original versions all along different numbers of sensors. Therefore, the proposed technique has an obviously enhanced efficiency. V. C ONCLUSIONS We have investigated a SFB-based U-root-MUSIC algorithm for fast DOA estimation based on a novel developed SFB averaging technique. The proposed SFB averaging technique allows efficient EVDs on two sub-matrices in about half sizes, and can be extended to reduce the complexity of othere unitary DOA estimators. Simulations show that SFB-based unitary methods has similar performances to their original versions while the former has an obviously enhanced computational efficiency than the latter.
[1] H. Krim and M. Viberg, “Two decades of array signal processing research: The parametric approach,” IEEE Signal Process. Mag., vol. 13, no. 4, pp. 67–94, Jul. 1996. [2] J. Zhao, F. Gao, W. Jia, S. Zhang, S. Jin, and H. Lin, “Angle domain hybrid precoding and channel tracking for millimeter wave massive MIMO systems,” IEEE Trans. Wireless Commun., vol. 16, no. 10, pp. 6868–6880, Oct. 2017. [3] H. Lin, F. Gao, S. Jin, and G. Y. Li, “A new view of multi-user hybrid massive MIMO: Non-orthogonal angle division multiple access,” IEEE J. Sel. Areas Commun., vol. 35, no. 10, pp. 2268–2280, Oct. 2017. [4] D. Fan, F. Gao, G. Wang, Z. Zhong, and A. Nallanathan, “Angle domain signal processing-aided channel estimation for indoor 60-GHz TDD/FDD massive MIMO systems,” IEEE J. Sel. Areas Commun., vol. 35, no. 9, pp. 1948–1961, Sep. 2017. [5] A. Khabbazibasmenj, A. Hassanien, S. A. Vorobyov, and M. W. Morency, “Efficient transmit beamspace design for searchfree based DOA estimation in MIMO radar,” IEEE Trans. Signal Process. vol. 62, no. 6, pp. 1490–1500, Mar. 2014. [6] K. Yu, R. E. Hudson, Y. D. Zhang, K. Yao, C. Taylor, and Z. Wang, “Low-complexity 2D direction-of-arrival estimation for acoustic sensor arrays,” IEEE Signal Process. Lett., vol. 23, no. 12, pp. 1791–1795, Dec. 2016. [7] V. V. Reddy, M. Mubeen, and B. P. Ng, “Reduced-complexity superresolution DOA estimation with unknown number of sources,” IEEE Signal Process. Lett., vol. 22, no. 6, pp. 772–776, Jun. 2015. [8] K.-C. Huarng and C.-C. Yeh, “A unitary transformation method for angle-of-arrival estimation,” IEEE Trans. Signal Process., vol. 39, no. 4, pp. 975–977, Apr. 1991. [9] J. Selva, “Computation of spectral and root MUSIC through real polynomial rooting,” IEEE Trans. Signal Process., vol. 53, no. 5, pp. 1923–1927, May 2005. [10] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol. AP-34, no. 3, pp. 276–280, Mar. 1986. [11] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” IEEE Trans. Signal Process., vol. 37, no. 7, pp. 984–995, Jul. 1989. [12] B. D. Rao and K. V. S. Hari, “Performance analysis of root-music,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 12, pp. 1939–1949, Dec. 1989. [13] F. Yan, M. Jin, and X. Qiao, “Low-complexity DOA estimation based on compressed MUSIC and its performance analysis,” IEEE Trans. Signal Process., vol. 61, no. 8, pp. 1915–1930, Apr. 2013. [14] F. Yan, M. Jin, and X. Qiao, “Source localization based on symmetrical MUSIC and its statistical performance analysis,” Sci. China Inf. Sci., vol. 56, no. 6, pp. 1–13, Jun. 2013. [15] F.-G. Yan, B. Cao, J. J. Rong, Y. Shen, and M. Jin, “Spatial aliasing for efficient direction-of-arrival estimation based on steering vector reconstruction,” EURASIP J. Adv. Signal Process., vol. 2016, no. 1, p. 121, Dec. 2016. [16] F.-G. Yan, M. Jin, S. Liu, and X.-L. Qiao, “Real-valued MUSIC for efficient direction estimation with arbitrary array geometries,” IEEE Trans. Signal Process., vol. 62, no. 6, pp. 1548–1560, Mar. 2014. [17] M. Haardt and J. A. Nossek, “Unitary ESPRIT: How to obtain increased estimation accuracy with a reduced computational burden,” IEEE Trans. Signal Process., vol. 43, no. 5, pp. 1232–1242, May 1995. [18] G. H. Golub and C. H. Van Loan, Matirx Computations. Baltimore, MD, USA: The Johns Hopkins Univ. Press, 1996. [19] R. Cao, B. Liu, F. Gao, and X. Zhang, “A low-complex one-snapshot DOA estimation algorithm with massive ULA,” IEEE Commun. Lett., vol. 21, no. 5, pp. 1071–1074, May 2017. [20] G.-T. Pham, P. Loubaton, and P. Vallet, “Performance analysis of spatial smoothing schemes in the context of large arrays,” IEEE Trans. Signal Process., vol. 64, no. 1, pp. 160–172, Jan. 2016. [21] M. Pesavento, A. B. Gershman, and M. Haardt, “Unitary root-MUSIC with a real-valued eigendecomposition: A theoretical and experimental performance study,” IEEE Trans. Signal Process., vol. 48, no. 5, pp. 1306–1314, May 2000. [22] D. A. Linebarger, R. D. DeGroat, and E. M. Dowling, “Efficient direction-finding methods employing forward/backward averaging,” IEEE Trans. Signal Process., vol. 42, no. 8, pp. 2136–2145, Aug. 1994. [23] P. Stoica and M. Jansson, “On forward–backward MODE for array signal processing,” Digit. Signal Process., vol. 7, no. 4, pp. 239–252, Oct. 1997. [24] P. Stoica and A. Nehorai, “Performance study of conditional and unconditional direction-of-arrival estimation,” IEEE Trans. Acoust., Speech Signal Process., vol. 38, no. 10, pp. 1783–1795, Oct. 1990.
