TOA and DOA Estimation for IR-UWB Signals: An Efficient Algorithm by Improved Root-MUSIC

The automatically paired time of arrival (TOA) and direction of arrival (DOA) can be jointly estimated via a high-precision multidimensional spectral peak search(SPS-) based multiple signal classification (MUSIC) algorithm in the impulse radio ultrawideband (IR-UWB) positioning system, while heavy computational burden is required. To tackle this issue, we propose an improved root-MUSIC algorithm for joint TOA and DOA estimation. After modelling the frequency domain form of the received signal, the algorithm first uses the signal subspace to establish the relationship between the two antennas. Then, the MUSIC spatial spectrum function is reconstructed with this relation, which enables it to offer a spectrum function in regard to the one-dimensional (1D) parameter of time delay. For further reducing the complexity, the TOA estimates of one antenna are obtained via 1D polynomial root finding instead of SPS, and the TOA estimates of the other antenna can be calculated by the established relationship. Finally, the DOA estimation can be achieved with the estimated TOAs. Due to the relationship between two antennas with signal subspace, the parameters estimated by the proposed algorithm are autopaired. Numerical simulations substantiate the superiority of the proposed algorithm.


Introduction
Impulse radio ultrawideband (IR-UWB) technique, which yields the advantages of a high data rate, low power consumption, low computational complexity, low power spectrum, and great immunity to multipath fading, has been widely applied in various fields, e.g., radar, location, and wireless communications [1][2][3][4][5]. An IR-UWB signal has great potential in precise positioning applications due to its extremely narrow pulse width [6]. Time of arrival (TOA) estimation-based ranging technology plays an important role in the IR-UWB positioning system, as it can make full use of the extremely high time resolution of the IR-UWB signal. And the acquisition of direction of arrival (DOA) is of benefit to reduce the number of nodes in the positioning system. Therefore, the joint TOA and DOA estimation has attracted extensive interest in recent years.
For TOA estimation, the parameter estimation method based on time domain signal processing mainly includes a coherent detection method using pulse template matched filter [7], incoherent TOA estimation algorithm based on energy detection [8][9][10][11], and two-step TOA estimation algorithm based on energy detection and coherent detection [12]. These algorithms obtain the TOA estimation with the arrival time of the direct path (DP) component. Unfortunately, due to the multipath effect of the UWB signal, especially in nonline-of-sight (NLOS) conditions, the resolution of the aforementioned methods declines because DP is no longer the strongest path (SP) in general. To address this issue, various TOA estimation approaches with superresolution based on frequency domain processing were proposed, such as the minimum norm spectrum estimation algorithm [13], multiple signal classification algorithm [14], and propagation operator method (PM) [15]. These algorithms first model the frequency domain channel impulse response and subsequently employ spectral peak search (SPS) to obtain TOA estimation, which leads to the increase of system complexity and the deterioration of real-time performance.
Numerous DOA estimation approaches have been developed for DOA estimation for UWB signals. A kind of data model similar to the data model of narrowband DOA estimation was presented in [16], which processes the received signal in the frequency domain and then conducts the MUSIC algorithm to achieve DOA estimation. In [17], an iterative quadratic maximum likelihood DOA estimation method is proposed, and a beamforming DOA estimation algorithm with frequency invariant in the frequency domain based on beamspace was studied in [18]. Considering the fact that the UWB signal has high time resolution, extracting DOA information from TOA can lead to high DOA estimation accuracy; many researchers have turned to the study of joint TOA and DOA estimation algorithms [19][20][21][22]. A low-complexity matrix pencil algorithm was proposed in [19,20], whereas the estimation performance is susceptible to noise since only one snapshot data is used, and extra parameter pairing of the estimated TOAs is required. In [21], the coarse TOA estimation is obtained by energy estimation and minimum distance criterion; then, the refined TOAs are estimated by the threshold comparison method based on the signal power delay spectrum, and finally, the DOA estimation is achieved by minimum variance unbiased estimation with the estimated TOAs of each antenna. A three-step joint TOA and DOA estimation algorithm was studied in [22], where the conventional threshold correlation approach is employed to obtain the time delay for the signal to reach each antenna, and then, the least mean square algorithm is employed to estimate the TOA and the time difference of arrival (TDOA) jointly. Finally, more accurate DOAs can be obtained by improving the accuracy of the TDOA. In [23], a successive PM algorithm was proposed to reduce the computational burden of the conventional 2D-PM algorithm, which first obtains coarse initial TOA estimation with low complexity and then performs two one-dimensional (1D) local SPS with the initial TOA estimation. The twodimensional (2D) MUSIC algorithm is able to deal with the TOA pairing problem in [19,20], which can estimate automatically paired TOAs with high resolution, whereas the 2D total SPS leads to a tremendous amount of computation. A root-MUSIC-based algorithm was presented in [24] to tackle this issue, where the closed-form solutions of TOA estimation can be obtained based on the polynomial root principle, and then, the DOAs can be calculated with the estimated TOA. However, since the received signals of the two antennas are processed separately, the cross-correlation information between the two antennas is neglected and extra pairing of TOA estimates is involved.
In this paper, the issue of joint TOA and DOA estimation in the IR-UWB system is investigated and a subspace-based improved MUSIC algorithm is proposed. Specifically, the proposed algorithm first models the received signals in the frequency domain and uses the signal subspace to establish the relationship between the two antennas. Furthermore, the MUSIC spatial spectrum function is reconstructed with this relationship, in which the two-dimensional parameter estimation is transformed into one-dimensional parameter estimation, which can be solved via polynomial root finding with remarkably reduced complexity. The proposed algorithm requires no spectral peak search and can directly give the closed-form solutions of the estimated parameters. Owing to the established relationship with signal subspace, the algorithm can estimate autopaired TOAs. Simulation results demonstrate that the proposed algorithm has significantly better parameter estimation performance than the matrix pencil, PM, and ESPRIT algorithms, and its estimation performance is close to that of the 2D-MUSIC and root-MUSIC algorithms.
The main contributions of our work are summarized as follows: (1) We propose an improved root-MUSIC algorithm for joint TOA and DOA estimation in the IR-UWB positioning system. By constructing and exploiting the relationship between the received signal of the two antennas, the 2D SPS is transformed into 1D SPS and the polynomial root finding is performed to further reduce complexity (2) The proposed algorithm can fully use the signal subspace and noise subspace to obtain autopaired TOA estimates (3) The proposed algorithm possesses close estimation performance to the MUSIC and root-MUSIC algorithms and outperforms the matrix pencil, PM, and ESPRIT algorithms The rest of this paper is organized as follows. The data model is introduced in Section 2, and Section 3 presents the proposed algorithm. In Section 4, we analyze the complexity of the proposed algorithm and list the advantages of the algorithm. Section 5 exhibits the simulation results, and the conclusions are drawn in Section 6.

Data Model
Assume that the UWB positioning system uses the second derivative of the Gaussian pulse as the UWB transmission signal. The transmission signal adopts Direct Sequence Binary Phase Shift Keying (DS-BPSK) modulation and can be expressed by [25,26] where N c represents the pulse repetition number of a single binary data symbol, b i ∈ f−1,+1g is the modulated binary data symbol sequence, and c n ∈ f−1,+1g represents a pseudorandom sequence used to implement multiple access communication. T s is the period of a binary data symbol, and T c is the pulse repetition period. pðtÞ denotes the second derivative of the Gaussian pulse and can be represented by where Γ is the pulse forming factor related to the pulse width.

Wireless Communications and Mobile Computing
UWB signals generally propagate in a multipath environment. According to the well-known Saleh-Valenzuela (S-V) model [27], we assume that a pulse of the transmitted signal generates multiple multipath components after passing through the channel, and these multipath components arrive at the receiver in the form of clusters. Suppose that the signal generates K clusters through the channel, and each cluster has L multipaths; the channel impulse response model of the kth cluster of the UWB channel is where α l ðkÞ stands for the channel attenuation coefficient in the kth cluster of the lth path and obeys Rayleigh distribution; the phase ϕ l ðkÞ is a random variable uniformly distributed on ½0, 2π. δðtÞ denotes the Dirac function. τ l ðkÞ represents the channel delay of the lth path in the kth cluster.
In other words, it can be assumed that the multipath delay of the k-cluster received signal is the same, that is, τ l ðkÞ = τ l . Let α l ðkÞ e jϕ l ðkÞ = β l ðkÞ denote the amplitude of random complex fading; then, (3) can be rewritten as The time domain form in the kth cluster signal received by the UWB system can be also written as where w ðkÞ ðtÞ is the additive white Gaussian noise and * represents convolution.

The Proposed Joint TOA and DOA Estimation Algorithm
The UWB system model is shown in Figure 1, where only two antennas are required. Assuming that the UWB signals are in the far field of the antenna array, the signals incident to the array can be regarded as a beam of parallel waves. Let Y 1 and Y 2 be the received signals in the frequency domain of the two antennas, respectively. Let τ = ½τ 1 , τ 2 ,⋯, τ L and ς = ½ς 1 , ς 2 ,⋯,ς L denote the vectors of TOAs corresponding to the two antennas, τ l and ς l stands for the TOAs associated with the lth ðl = 1, 2,⋯,LÞ path. Then, according to the signal model in Section 2, Y 1 and Y 2 can be written as where Β = ½βð1Þ, βð2Þ,⋯,βðKÞ ∈ ℂ L×K is a complex fading coefficient of the channel. W 1 = ½w 1 ð1Þ,⋯,w 1 ðKÞ ∈ ℂ N×K and W 2 = ½w 2 ð1Þ,⋯,w 2 ðKÞ∈ℂ N×K is the noise samples received by antenna 1 and antenna 2. E τ and E ς are the delay matrix which can be represented by

Wireless Communications and Mobile Computing
Denote as the difference of the TOAs of the lth path. As can be concluded from Figure 1, Δb τ l can be obtained from where d represents the distance between the two antennas, θ l is the DOA of the lth path, and c denotes the speed of light. From (12), the DOA estimates can be calculated by It can be seen from (13) that TOAs of two antennas should be estimated first to achieve DOA estimation.
3.1. Relationship Conductiond. Before estimating TOA, we first use the signal subspace to establish the relationship between the two antennas. Construct the overall output matrix Z ∈ ℂ 2N×K with the received signal matrix in the frequency domain of two antennas Y 1 ∈ ℂ N×K and Y 2 ∈ ℂ N×K as Let, and Aðτ, ςÞ can be specified by [23] where aðτ l , ς l Þ = ½Sðω 0 Þ, Then, (15) can be represented by The covariance matrix of Z can be approximately calculated byR = ZZ H /K. By performing eigendecomposition,R can be decomposed aŝ whereÛ s ∈ ℂ 2N×L is the signal subspace [28] spanned by the eigenvectors corresponding to the largest L eigenvalues and U N ∈ ℂ 2N×ð2N−LÞ is the noise subspace [29] consisting of the eigenvectors corresponding to the remaining N − L eigenvalues. According to subspace theory, in noise-free condition, the signal subspace U s can be represented by where T ∈ ℂ L×L is a nonsingular matrix. Then, we partition U s into two parts as follows: where U s1 ∈ ℂ N×L , U s2 ∈ ℂ N×L denote the signal subspace of antenna 1 and antenna 2, respectively. As U s1 = A 1 T and U s2 = A 2 T, define a transformation matrix H as where ð⋅Þ + represents the pseudoinverse operation. As T is full rank, then H = A 2 TT −1 A + 1 = A 2 A + 1 and the relationship Wireless Communications and Mobile Computing between A 1 and A 2 can be represented by It can be deduced from (22) and (23) that the relation between A 1 and A 2 can be obtained with the signal subspace U s of the overall output, which can be utilized to simplify the two-dimensional parameter estimation to one-dimensional estimation.

TOA Estimation Based on Root-MUSIC Algorithm.
According to the orthogonal relationship between each column vector of A and noise subspace U N , the 2D-MUSIC spatial spectrum function can be constructed as where aðτ, ςÞ, a 1 ðτÞ, and a 2 ðςÞ are the column vectors of A ðτ, ςÞ, A 1 , and A 2 , respectively. By locating the largest L peaks of P 2D-MUSIC ðτ, ςÞ, τ and ς that correspond to the peaks are exactly to be estimated, whereas expensive computational cost is required due to the two-dimensional global spectral peak search. In addition, the polynomial rootfinding technique is difficult to be implemented directly since a two-dimensional parameter to be estimated is involved. To tackle this issue, we use the relationship established between A 1 and A 2 to transform the two-dimensional SPS into one-dimensional SPS, and then, the polynomial root-finding technique is exploited to replace SPS, which can reduce the computational complexity significantly. According to (23), we can derive the relation a 2 ðςÞ = Ha 1 ðτÞ; then, (24) can be reconstructed as where Q = ½I N , H H Û N , a 1 ðτÞ = SbðτÞ, and bðτÞ = ½1, e −jΔωτ ,⋯,e −jðN−1ÞΔωτ T represent the column vector of delay matrix E τ . As can be seen from (25), only 1D SPS for τ is involved to obtain the TOA of the first antenna, where the complexity can be reduced substantially compared with 2D SPS. Subsequently, the polynomial root-finding technique is employed to further reduce complexity. Denote z = e −jΔωτ , then we have bðτÞ = bðzÞ = ½1, z,⋯,z N−1 T as well as a 1 ðτÞ = a 1 ðzÞ = SbðzÞ, the delay estimation of (25) can be thereby transformed into the following polynomial rooting problem: The L rootsẑ 1 ,ẑ 2 , ⋯,ẑ L of (26) nearest to the unit circle are exactly required, and z l = e −jΔωτ l . Accordingly, the TOA estimates of antenna 1 can be calculated by Once the rootsẑ l = e −jΔωτ∧ l are obtained,â 1 ðτ l Þ = Sbðz l Þ = S½1, z∧ l ,⋯,z∧ l N−1 T can be calculated and subsequentlŷ a 2 ðς l Þ can be estimated bŷ By usingÂ 2 ðςÞ = ½â 2 ðς 1 Þ,â 2 ðς 2 Þ,⋯,â 2 ðς L Þ,Ê ς can be estimated viaÊ According to (12), define E a and E b as matrices composed of the first and the last N − 1 rows of E ς , respectively, and we can get ⋯,e −jΔως L Þ. Correspondingly, defineÊ a andÊ b as the matrices composed of the first and the last N − 1 rows ofÊ ς in (29), the TOAs of antenna 2 can be obtained by where ð⋅Þ + represents pseudoinverse operation. It is worth noting that the TOA estimates b τ l and b ς l are automatically paired, as the relationship between A 1 and A 2 is exploited in (28) with H.
3.3. The Procedure of the Proposed Algorithm. The major steps of the proposed improved root-MUSIC algorithm in the IR-UWB system are as follows: (1) CalculateR of the total received frequency domain signal Z in (15), then obtainÛ s andÛ N with eigenvalue decomposition ofR (2) Compute H and reconstruct the MUSIC spatial spectral function according to (22)- (25) (3) Construct the polynomial in (26) and conduct polynomial rooting to find the L largest rootsẑ 1 ,ẑ 2 , ⋯, z L to estimate TOA estimates corresponding to antenna 1 b τ l ðl = 1, 2,⋯,LÞ via (27) (4) Calculateâ 2 ðς l Þ andÊ ς via (28) and (29) (5) Obtain the TOA estimates corresponding to antenna 2 b τ l through (30) 5 Wireless Communications and Mobile Computing (6) Calculate Δb τ l = b ς l − b τ l and estimate the DOAs b θ l according to (14) Remark 1. It is assumed that the multipath rays L are known in this paper. If L is unknown, it can be estimated through the method in [30].

Remark 2.
We assume that the number of rays in each cluster is the same, since it is generally believed that the channel remains stable and the rate of change of the channel is very slow during the observation time.

Performance Analysis
4.1. Complexity Analysis. The complexity of the proposed algorithm is mainly caused by the transformation of the signals into the frequency domain, calculation of covariance, eigenvalue decomposition of covariance, and polynomial rooting. Specifically, the complexity of transforming the transmitted and received signals into the frequency domain model is OfN 2 + 2KN 2 g. Calculating covarianceR needs O f4KN 2 g, and the complexity of eigenvalue decomposition is Of8N 3 g. Calculating the transformation matrix H = U s2 U + s1 requires OfN 3 + 2LN 2 + L 2 Ng, and construct Q needs Of2N 2 ð2N − LÞg. The procedure of polynomial rooting costs OfðN − 1Þ 3 g. Consequently, the total complexity of the proposed algorithm is OfðN − 1Þ 3 + 13N 3 + ð6K + 1ÞN 2 + L 2 Ng, while the 2D-MUSIC algorithm costs Of8N 3 + ð6K + 1ÞN 2 + m 2 ð2N + 1Þð2N − LÞg, where m is the peak search times. And the root-MUSIC algorithm [24] requires Of10 Figure 2 shows the complexity of the above methods, where m = 2000, K = 100, and L = 3. It can be concluded from Figure 2 that the proposed algorithm requires the lowest computational complexity. And the complexity of the polynomial rooting-based algorithm is much lower than that of the SPS-based algorithm.  (1) The proposed algorithm can achieve automatically paired two-dimensional TOA estimation (2) The proposed algorithm has an estimation performance close to that of the 2D-MUSIC algorithm, but remarkably lower computational burden is required, since the 2D SPS is transformed into 1D polynomial rooting (3) The proposed algorithm possesses approximate estimation performance compared with the root-MUSIC algorithm and outperforms the matrix pencil [19,20], ESPRIT, and PM [15] algorithms

Simulation Results
In this section, simulations are provided to evaluate the estimation performance of the proposed algorithm. The signal to noise ratio (SNR) is defined by SNR = 10 lg ðkyðtÞk 2 F / kwðtÞk 2 F Þ, where yðtÞ represents the received time domain signal, wðtÞ is the received noise, and k⋅k F stands for Frobenius norm. The root mean squared error (RMSE) is presented for evaluation which is defined as where χ l and b χ l,m denote the theoretical parameter, i.e., TOA or DOA, and the estimate of the lth path in the mth trial, respectively. M denotes the number of Monte Carlo trials.
In the simulations, we assume that the UWB pulse wave function is pðtÞ = exp ð−2πt 2 /Γ 2 Þð1 − 4πt 2 /Γ 2 Þ, where Γ = 0:25 ns represents the shaping factor for the pulse. The repetition of every symbol is N c = 5, the chip duration is T c = 2 ns, and the symbol duration is T s = N c T c = 10 ns. We plot the UWB pulse wave function pðtÞ and the transmitted BPSK-UWB signal sðtÞ in Figures 3 and 4, respectively. We suppose that the UWB channel parameters β ðkÞ l are known, which are random complex fading amplitudes. There are N = 64 frequency samples, K = 100 observed clusters, and L = 3 rays of BPSK-UWB arriving signals whose TOAs and DOAs are 0.2 ns, 0.3 ns, and 0.4 ns and 10°, 30°, and 45°, respectively. Remark 3. Besides the Gaussian pulse, some other types of pulse wave functions for use in IR-UWB communications have been proposed such as Laplacian monopulse waveshape [31] and Hermite pulse shapes [32]. In practice, however, the second derivative of a Gaussian pulse is widely used [33]. The pulse waveform can be changed by regulating the shaping factor Γ.     are random time-varying functions. However, the rate of the variations is very slow compared to any useful signaling rate that is likely considered, e.g., higher than tens of kbit/s. Thus, these parameters can be treated as virtually time-invariant random variables.
Remark 5. We suppose that each cluster has the same number of rays. Normally, we consider that the variation rate of the channel is slow and the channel remains stable in the observing time, so we can assume that the corresponding delays in each cluster are the same. Remark 6. We set the number of frequency samples as an integer power of 2. As a result, fast Fourier transform (FFT) can be used, which is more efficient than DFT.
We provide the parameter estimation results of the proposed algorithm in Figures 5 and 6, where SNR = 20 dB and M = 50. As illustrated in Figures 5 and 6, the proposed algorithm can estimate TOAs and DOAs accurately. Figures 7 and 8 compare the parameter estimation performance with different methods, including the proposed algorithm, 2D-MUSIC, root-MUSIC [24], matrix pencil [19,20], PM [15], and ESPRIT algorithms, where N = 64, L = 3, K = 100, and M = 500. It is observed that the proposed algorithm, root-MUSIC algorithm, and 2D-MUSIC algorithm outperform the other algorithms in parameter estimation. The matrix pencil algorithm performs worst, as this algorithm is proposed for single snapshot data, and its estimation performance is very sensitive to noise. In the case of multiple snapshots, we average the estimated values under each snapshot, while its performance is still poor. The PM algorithm has poor performance at low SNR; this is because the frequency domain division of the received signal and the transmitted signal is used to obtain the frequency domain channel impulse response in [15], which enhances the noise at low SNR. Furthermore, the proposed algorithm possesses the same parameter estimation as the root-MUSIC algorithm which is close to that of the 2D-MUSIC algorithm, while the proposed algorithm has lower computational complexity.
The RMSE results with different L of the proposed algorithm are provided in Figures 9 and 10, where N = 64, K = 100, and M = 500. As Figures 9 and 10 demonstrate, the parameter estimation performance of the proposed algorithm improves with the decrease of L.

Conclusions
In this paper, we have presented a subspace-based improved root-MUSIC algorithm with an IR-UWB system for joint TOA and DOA estimation. The relation between the two received antennas is first established by partitioning the signal subspace, with which the 2D spatial spectral function can be converted into a 1D one. Then, the polynomial root-finding technique is exploited to obtain the TOA of one antenna instead of SPS for further reducing the computational burden, and the TOA of the other antenna can be estimated by utilizing the relation between the two received antennas. Finally, the DOAs can be calculated with the difference of the TOAs. Closed-form solutions to the parameter estimation can be attained, and no extra parameter pairing is required. Simulations demonstrate the validity of the proposed algorithm for the IR-UWB system with regard to computational complexity and parameter estimation performance.

Data Availability
For reasonable requests, data supporting the results of this study can be obtained from the corresponding author.

Conflicts of Interest
The authors declare no conflict of interest.  Wireless Communications and Mobile Computing