Discrete fourier transform-based TOA estimation in UWB systems

In this paper, we propose two time of arrival estimators for ultra wideband signals based on the phase difference between the discrete Fourier transforms (DFT) of the transmitted and received signals. The first estimator is based on the slope of the unwrapped phase and the second one on the absolute unwrapped phase. We derive the statistics of the unwrapped phase. We show that slope-based estimation almost achieves asymptotically the baseband Cramer-Rao lower bound (CRLB), while the absolute-phase-based estimator achieves asymptotically the passband CRLB. We compare the proposed estimators to the time-domain maximum likelihood estimator (MLE). We show that the MLE achieves the CRLB faster than the DFT-based estimator, while the DFT-based estimator outperforms the MLE for low signal to noise ratios. We describe also how to use the proposed estimators in multipath UWB channels.


I. Introduction
UWB has received increasing attention for many applications like positioning since the FCC (Federal Communications Commission) allowed in 2002 the unlicensed use of the spectrum between 3.1 and 10.6 GHz [1].
Thanks to their ultra wideband (UWB) larger than 500 MHz, UWB signals can be used for highly accurate positioning using the time of arrival (TOA) technique. Many TOA estimators have been proposed in the literature, especially for impulse radio UWB (IR-UWB) signals. Most proposed estimators like the maximum likelihood estimator (MLE), the energy-based estimators, the autocorrelation-based estimators, the threshold-based estimators, and others are based on the time domain [2][3][4][5][6][7][8][9][10]. The drawback of time-domain estimators is that their precision is limited by the sampling frequency being used, and complex interpolation is required in order to improve the performance. Some other estimators for either electromagnetic or acoustic signals are using the discrete Fourier transform (DFT) of the received signal [11][12][13][14][15][16][17].
In this paper, we propose two estimators for the TOA based on the phase of the DFT of the received signal. The first estimator is based on the slope of the phase and the second one, on the absolute phase. For both estimators, we first compute local estimates at the different frequency components, and then we combine them in order to find the global estimates.
The main three contributions of this work are that: • we show that using the DFT, we can achieve asymptotically the CRLBs (Cramer-Rao lower bound) using very simple estimators requiring only few samples and a sampling rate equal to the signal bandwidth. In our approach, the sampling period is much larger than the achieved accuracy, while in time-domain-based estimation, the sampling period must be smaller than the required accuracy. Another advantage of DFT-based estimation is that we do not need to identify the main lobe of the autocorrelation of the used pulses like in time-domain estimation.
• we show that the MLE achieves the CRLB faster than the DFT-based estimator, while the DFT-based estimator outperforms the MLE for low SNRs.
• we compute the statistics of the unwrapped phase of a noisy signal.
The main difference between this work and the previous works using the DFT approach is that in the previous works, the TOA is not estimated based on the phase of the DFT [13,14], or the problem of phase ambiguity is not investigated (by assuming the maximum time delay smaller than the period of the highest frequency component) [15,16], or the problem of phase ambiguity is solved using other approaches (Chinese remainder theorem [11,12] or recursive correction of the TOA estimate [17]). The proposed estimators can be used for IR-UWB signals as well as for multi-carrier UWB (MC-UWB) signals. Note that the main goal of this paper is to give the main ideas about DFT-based TOA estimation. Many improvements can be introduced in order to make the proposed estimators achieve performance closer to the CRLBs.
The paper is organized as follows. In Section II, we describe the system model. In Section III, we consider the MLE of the local phase and compute the statistics of the unwrapped phase. In Section IV, we derive the local slope-based and absolute-phase-based TOA estimators. In Section V, we derive the global slope-based and absolute-phase-based TOA estimators. In Section VI, we show how multipath UWB channel can be handled.

II. System model
We consider a transmitter and a receiver communicating through an additive white Gaussian noise (AWGN) channel.
Denote by s(t), r(t) and n(t) the complex envelopes (baseband) of the transmitted signal, the received signal and the AWGN, filtered around central frequency f c with a bandwidth B ([f c -B/2, f c + B/2]). r(t) can be written as: where a and τ are the gain and the time delay of the channel, and s τ (t) = s(t -τ). After sampling at the rate B, we get: where z[m] denotes the sample of the signal z(t) at t = mT s (T s = 1/B is the sampling period). (n[m]) is a white Gaussian sequence (i.e., the samples n[m] are independent and identically distributed (iid)). The variance of n [m] is given by σ 2 n = 2N 0 B where 2N 0 is the one-sided power spectral density of the AWGN.
Let R[k], (k = -M/2,..., M/2 -1) be the DFT of r[m]: also white Gaussian with a variance equal to σ 2 N = Mσ 2 n = 2MN 0 B [18]. As the Shannon sampling theorem is respected, and by assuming s(t) limited in time and s τ (t) falling in the period of observation, we can write: where S τ (f) and S(f) are the FTs of s τ (t) and s(t), respectively, and S[k] is the DFT of s[m] and: For simplicity reasons, we denote from now S[k], R[k] and N[k] by S k ,R k and N k , respectively. From (1) and (2), we can write R k as: where U k = αe −j2π (f c +f k )τ S k is the DFT of the useful part of the received signal. Denote by r Z ,θ Z , x z and y Z the modulus, phase, real part and imaginary part of any complex number Z. From (4), we can define k as: Given that N k is Gaussian, we can write the probability density function (PDF) of R k as: where σ 2 = σ 2 N /2 = MN 0 B is the variance of x N k and y N k .

III. Statistics of the unwrapped MLE of the phase
In this section, we consider the MLE of the phase and compute the statistics of its unwrapped version.
The joint log-likelihood function of ρ U k and k can be obtained from (5) and (6): The CRLBs of ρ U k and k are the diagonal elements of the inverse of the Fisher information matrix given by −E{(∂ 2 Λ ρ U k ,ϕ k /∂z i ∂z j )}, z i , z j ∈ {ρ U k , ϕ k }(E{·} denotes the expectation operator). The CRLB of k is given by: where ν k = ρ 2 U k /σ 2 = α 2 ρ 2 S k /σ 2 is the SNR obtained at f k . ν k is called the local (or instantaneous) SNR (corresponding to f k ). The global SNR is defined as: It is obvious that the time delay can be estimated from an estimation of (5) as either: (i) the phase to angular frequency ratio or (ii) the slope of the phase with respect to the angular frequency. For both approaches, the estimated phase must be continuous. With the former it must also be around the true value, while with th e latter a constant offset along the frequency axis is accepted. As in practice the phase is computed modulo 2π (wrapped phase), an unwrapped version of it is needed in order to rebuild the continuous phase.
In practice, the unwrapped phase can be obtained recursively by adding a multiple of 2π to the wrapped phase until the absolute difference between neighboring phases becomes less than or equal to π. Denote byφ k the wrapped MLE of the phase andφ k the unwrapped MLE. We can write the unwrap criterion as: where the non-ambiguity condition (2πΔfτ <π) must be respected. Unwrap procedure described above is well known and "unwrap" MATLAB function can be used to perform unwrapping.
As in practice the true value of the phase is unknown we can start the unwrap procedure from an arbitrary k 0 by takingφ k 0 =φ k 0 , then running the unwrap procedure for k 0 +1,..., M/2 -1 and k 0 -1,..., -M/2. It is obvious that the unwrapped phase may have an offset (almost constant ∀k) with respect to the true phase dependent on the offset at the starting point (2π Let us now consider how to obtain a wrapped estimationφ k of the phase. It can be obtained from (7) using a MLE and taking ∂Λ ρ U k ,ϕ k /∂ϕ k = 0: where {·}* denotes the complex conjugate. The estimatesφ k at different frequencies k are independent because the noise samples N k are independent.
As shown in [19], the PDF ofφ k can be obtained from that of θ R k by integrating (6) with respect to ρ R k : where erfc(z) = (2/ √ π ) +∞ z e −ξ 2 dξ denotes the comple-mentary error function, and the superscript wr the wrapped phase. T wr ϕ k (φ k ) is 2π periodic and can be defined on any interval It is shown in [20] that the distribution of the wrapped phase can be approximated by a normal distribution if the local SNR ν k is sufficiently high, and by a uniform distribution if ν k is very low.
As the unwrap criterion in (10) can be written as As forφ 0 =φ 0 we have D 0 = [-π, π], the domain D k is given by: (11): Note that the domain of Tφ k |φ k−1 (φ k ) depends oñ ϕ k−1 but not its expression. In order to express the marginal PDF ofφ k with respect to that ofφ k−1 , we first compute the joint PDF ofφ k andφ k−1 , and then we integrate with respect toφ k−1 taking into account thatφk ∈ Iφ k−1 (φ k − π ≤φ k−1 ≤φ k + π ): can be computed recursively for k = 1,..., M/2 -1 and k = -1,..., -M/2 using (13). Obviously, the starting point is Tφ 0 (φ 0 ) = T wr ϕ k (φ 0 ). The mean and the variance ofφ k are given by: In Figure 1a, we show the true phase k , a realization of the wrapped MLE of the phaseφ k (φ k ∈ [−π , π ]) , and the corresponding unwrapped MLEφ k versus f k + f c (number of samples M = 16, k = -8, ..., 7). The unwrap procedure is started here from k = -8. The transmitted signal is a cardinal sine (bandwidth B = 2 GHz) modulated by a carrier (f c = 2 GHz). We take τ = 2 ns, and ν = 17 dB (global SNR). We can see thatφ k is almost continuous with a phase offset almost constant with respect to the true phase.
However, some errors multiple of -2π can be introduced during the unwrap procedure as can be seen in Figure 1b, c for two other realizations of the of the wrapped phaseφ k . This happens when the unwrap procedure should add a multiple of 2π to the next phase (for instance at k = -3 in Figure 1b), but does not do it because the absolute difference between the neighboring noisy phases is less than π (|φ −3 −φ −4 | ≤ π ). Every time this phenomenon happens, an additional error of -2π will be introduced.
Note that errors multiple of 2π can also be introduced. This happens when the unwrap procedure should not add a multiple of 2π to the next phase, but does it because the absolute difference between the neighboring noisy phases is greater than π. These errors occur rarely if the slope of the true phase is positive.
Here we have started the unwrap procedure from k = 0. We can see in Figure 2b that for k = 15 (phase corrected at the end of the unwrap procedure), the PDF has three secondary lobes located at -4π, -2π, and 2π from the main lobe. The strongest one is that located at -2π.
As already mentioned, the presence of these secondary lobes is due to errors multiple of ±2π introduced by the unwrap procedure. The main lobe becomes weaker and secondary lobes stronger as the frequency increases which means that we have more chance that such an error occurs. This is due to the fact that the unwrapping is performed recursively for increasing frequencies (see Figure 1a-c), so the ±2π errors accumulate over the course of the procedure. If we increase the number of samples or decrease the global SNR, we will obtain more secondary lobes at · · ·, -4π, -2π, 2π, 4π, · · · from  the main lobe. Errors multiple of -2π (resp. 2π) are more frequent if the slope of the true phase is positive (resp. negative). Obviously, the unwrapped phase is biased, and both the bias and the variance increase with the frequency due to the accumulation of ±2π errors.
We can see in Figure 3a that the simulated variance and MSE closely follow the theoretical ones, which validates our theoretical approach. However, variance and MSE are not following the CRLB, and they increase with the frequency due to the errors multiple of ±2π which are introduced by the unwrap procedure.
In Figure 3b where the local SNRs are sufficiently high (ν k = 28 dB, ∀k), we can see that the derived and simulated variance and MSE are very close to the CRLB. In fact for high SNRs, the wrapped phases are unwrapped correctly because the errors multiple of ±2π become very rare.

IV. Slope-based and absolute-phase-based local TOA estimators
In the last section, we have studied the unwrapped MLẼ ϕ k of the phase k . In this section, we propose two local TOA estimators based onφ k .  In order to overcome the problem of the phase offset mentioned in subsection III, we define the first local estimator of τ based on the slope ofφ k : where the superscript bb denotes that τ is estimated based on the information carried by the baseband frequency components.τ bb k can be named either local slope-based or local baseband (BB) TOA estimator. By assuming (for simplicity reasons) that (φ k ) are independent (not true because of the unwrap procedure), the covariance and variance σ 2 τ bb k ofτ bb k can be written as: In order to benefit from the information on τ carried by the passband frequency components, and as the phase offset betweenφ k and k is multiple of 2π, we can estimate the phase offset by: where "round" denotes the "round to nearest integer" function, andτ bb the global slope-based estimator.τ bb is given in Section V as a linear combination ofτ bb k . As the phase offset is estimated, we can now define the second local TOA estimator from (5) and (18): whereτ pb k is named local absolute-phase-based or local pass-band TOA estimator. By assumingΔϕ equal to the true value (true for high SNRs), the variance of τ pb k can be written as: The local passband CRLB of τ can be obtained from (5) and (8): If we assume in (5) that 2π f c τ is a random phase (if phase uncertainty is introduced during the down-conversion of the signal), the local baseband CRLB can be written as: (22) As for sufficiently high SNRs, the unwrapped phase becomes unbiased and its variance converges to its CRLB (1/ν k ), we can deduce from (20) and (21) (resp. (17) and (22)) that the local passband (resp. baseband) TOA estimator becomes also unbiased and achieves the local passband CRLB (resp. the sum of the local baseband CRLB of f 0 and f k ).
In Figure 4,

V. Slope-based and absolute-phase-based global TOA estimators
In this section, we derive the global TOA estimators based on the local TOA estimators studied in section IV.
The global baseband (resp. passband) TOA estimator τ bb (resp.τ pb ) is defined as the minimum-variance unbiased linear combination of the local estimators τ bb k , k = −M/2, ..., M/2 − 1 (respτ pb k ). Consider M unbiased estimatorsς k of the same parameter ζ. The minimum-variance unbiased linear combination of (ς k ) is given by: where {·} T denotes the transpose operator, zthe vector a the sum of the elements of a, 1 = (1 · · · 1) T , and ς the covariance matrix ofς . The variance ofς is given by: From (17) and (23), we can obtain the global baseband estimator and its variance: Given that the covariance matrix ofτ − pb is diagonal (τ pb k assumed independent), we can write the global passband estimator and its variance as: As the covariances and variances of the local estimatorsτ bb k andτ pb k ) are unknown, we compute the global estimators from (24) and (25) by assuming thatφ k achieves the CRLB c ϕ k , and substituting σ 2 ϕ k by 1/ρ 2 S k (proportional to c ϕ k ). Given that N k in (4) is a white sequence, the global passband and baseband CRLBs of τ can be written as: where ν is the global SNR given in (9) and Letτ ml be the time-domain MLE of τ.τ ml is given by: where s pb (t) and r pb (t) denote the real passband transmitted and received signals and ⊗ the convolution operator.
In Figure 5, take T ml s = 1 ps (125 times smaller than the DFT-based sampling period). The MSEs presented here are obtained by simulation (noise generated 10,000 times).
We can see that the global baseband estimator almost achieves asymptotically the baseband CRLB. We can also see that both the MLE and the global passband estimator achieve asymptotically the passband CRLB. However,τ ml achieves c pb faster thanτ pb . Many improvements can be introduced to our estimators in order to make them achieve the CRLBs faster. Hereafter, we will describe briefly one more baseband estimator and one more passband estimator.
We have already seen that the unwrap procedure introduces sometimes errors multiple of -2π in the unwrapped phase. These errors seriously deteriorate our estimators. In order to overcome this problem, we consider first the following slope-based estimator: where Δf is given in (3) As the unwrapped phase errors described above generate large negative slopes, and as the time delay can be assumed positive by putting the reference pulse at the beginning of the observation period, we can mitigate these errors by keeping only the positive values ofτ sp k . Letτ A new global slope-based estimator can be obtained from (23) and (27): Now, instead of unwrapping the phase recursively, we unwrap eachφ k (wrapped phase) with respect to 2π (f k + f c )τ sp + in order to getφ sp k (new unwrapped phase located around the true phase). The new global baseband (resp. passband) estimatorτ bb + (resp.τ pb + ) is obtained as before from (16) and (24) (resp. (19) and (25)), but after substitutingφ k byφ sp k (resp.φ k +Δϕ byφ sp k ) in (16) (resp. (19)). The MSEs ofτ bb + andτ pb + obtained by simulation are shown in Figure 5. We can see thatτ bb + andτ pb + achieve c bb and c pb faster thanτ bb andτ pb , respectively. Still, the MLE achieves c bb and c pb faster thanτ bb andτ pb . However, for small SNRs (r < 15 dB), the new passband estimator outperforms the MLE.
Fianlly, the main advantage of the MLE is that it achieves the CRLB faster, while the main two advantages of the new estimator are that: i) it requires a sampling rate and a number of samples much smaller than those required by the MLE and that ii) it outperforms the MLE for low SNRs.

VI. TOA estimation in multipath channels
Assume now that we have a multipath UWB channel. The baseband channel impulse response can be written as: where a (l) and τ (l) are the gain and the delay of the lth MPC. The baseband signal received through the multipath channel can be written as:

r MP (t) = s(t) ⊗ h(t).
Let Γ r MP, s (t) be the cross-correlation function of the modulus of the baseband transmitted and received signals. Γ r MP, s (t) can be written as: Γ r MP, s (t) = |r MP (t)| ⊗ |s(−t)| where we have considered the modulus in order to get only one peak per MPC (the used baseband pulse must have only one lobe). The coarse estimates of τ (l) can be obtained as locations of the peaks of Γ r MP, s (t) crossing a given threshold. Once the coarse estimates are obtained, we can apply our DFT-based estimators by taking a window around each MPC slightly larger than the pulse width. The final estimates of τ (l) are expected to have the same characteristics shown throughout this paper if the MPCs are not overlapping.

VII. Conclusion
Two TOA estimators are proposed based on the absolute phase and the slope of the unwrapped phase of the DFT of the received signal. The slope-based TOA estimation is used as a coarse estimation in order to rebuild the absolute unwrapped phase and to compute the absolute-phase-based estimator. The statistics of the unwrapped phase are computed. It has been shown that the slope-based estimator almost achieves asymptotically the baseband CRLB, while the absolute-phase-based estimator achieves asymptotically the passband CRLB. The proposed estimators are compared to the time-domain MLE estimator. It has been shown that the MLE achieves the CRLB faster than the DFT-based estimator, while the DFT-based estimator outperforms the MLE for low SNRs. It has also been also described how the proposed estimators can be used in multipath UWB channels. The main theoretical results are validated by simulation.