Energy-Based TOA Estimation

838 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 7, NO. 3, MARCH 2008 Transactions Papers Energy-Based TOA Estimation Antonio A. D’Amico, Umberto Mengali, and Lorenzo Taponecco Abstract— This paper investigates time of arrival estimation via impulse radio ultra-wideband technology. A dense multipath channel is assumed and frequency-selective propagation effects are taken into account that make the individual received pulses different in shape from the transmitted monocycles. In these conditions a conventional correlation-based estimator cannot be used and other schemes must be resorted to. Assuming only an approximate knowledge of the received pulses’ duration, a TOA estimator is derived making use of least mean square techniques. The algorithm operates on energy measurements made on a folded version of the received signal. Its performance is investigated by simulation under implementation conditions of different complexity. In particular, sub-Nyquist sampling rates are considered and analog-to-digital conversion with limited resolution is explored. Comparisons are made with other available schemes. Index Terms— Synchronization, frame frequency estimation, ultra-wideband communications. I. I NTRODUCTION N the last few years ultra-wideband (UWB) communications have been the focus of intense research efforts aimed at the exploitation of the vast unlicensed spectrum made available by the FCC Report and Order dated April 2002. As a result, the conviction has been reached that UWB technology is a viable solution for high-rate short-range wireless communications, as well as for low-rate moderaterange communications with ranging capabilities. In particular, it is believed that the impulse radio version of UWB (IRUWB) can provide centimeter accuracy in ranging and can be exploited in a variety of applications aimed at localizing valuable assets in hospitals, industrial areas, airports and so on [1]. The key toward this goal is in the ultrashort pulses (in the nanosecond scale) employed in IR-UWB radio. On a white additive Gaussian noise (AWGN) channel their time of arrival (TOA) can be measured within the accuracy of the Cramer-Rao bound by correlating the received signal with a template waveform shaped as the transmitted pulse and looking for the time shift of the template that yields the largest I Manuscript received August 4, 2006; revised November 14, 2006; accepted January 19, 2007. The associate editor coordinating the review of this paper and approving it for publication was Y. R. Zheng. The authors are with the Department of Information Engineering, Via Caruso, 56100 Pisa, Italy (e-mail: {antonio.damico, umberto.mengali, lorenzo.taponecco}@iet.unipi.it). Digital Object Identifier 10.1109/TWC.2008.060545. correlation [2]. Unfortunately, an AWGN channel is a poor model for a typical indoor environment where propagation takes place with hundreds of distinct trajectories [3] and the received pulses are distorted because of frequency-selective effects due to scattering [4] and frequency variations of the antenna pattern [5]. Distortions may affect even the directpath pulse (DPP) if building materials are present between transmitter and receiver [6]. In summary, the transmission of a single pulse (monocycle) generates a train of echoes at the receiver, with different shapes and random positions over tens of nanoseconds. In these conditions the TOA estimation problem is much more demanding for the following reasons1 : (i) On a multipath channel the goal is the location of the earliest pulse, not of an isolated pulse. To see the difference, suppose that there are no distortions and the correlator’s template is the transmitted monocycle. Then, the TOA of the DPP is found as the instant where the correlation overcomes a given threshold for the first time. On the other hand, with a single-path channel, the TOA of the (unique) pulse is found as the position of the highest correlation peak, wherever it happens to occur over the observation interval. (ii) As the arriving pulses are distorted and the estimation of their shape is computationally intensive [9], the implementation of a correlation-based TOA estimator is practically impossible. Significant degradations are incurred taking the monocycle as a template [10]. (iii) Pulse overlaps are likely to occur. If the DPP is partially overlapped with the next pulse, a correlation-based method gives wrong results even with perfect knowledge of the pulse shapes. In conclusion, TOA estimation is a challenging task and is currently the focus of intense research. The following is a representative sample of the work published in this area. References [11]-[12] investigate optimal correlation-based TOA estimators. The received pulses are undistorted and their shape is known. It is shown that a particular acquisition strategy, referred to as “bit reversal” search, provides quite shorter acquisition times than a “linear” search where the cells are tested in order. Reference [13] approaches TOA estimation adopting a simplified version of the generalized maximum likelihood 1 In this paper we assume that the straight-line from transmitter to receiver is a viable propagation path. Visual obstructions are allowed provided that the received signal is still detectable. Techniques to identify non-line-of-sight propagation are discussed in [7]-[8] c 2008 IEEE 1536-1276/08$25.00  D’AMICO et al.: ENERGY-BASED TOA ESTIMATION criterion. Again, the received pulses are undistorted and their common shape is used for correlation with the incoming signal. As the direct path need not to be the strongest, the estimation process is split in two parts. The TOA of the strongest pulse is found first, looking for the highest peak in the correlator output. If such a pulse is not the earliest that arrives at the receiver, then TOAs and amplitudes of the prior pulses are computed in succession, proceeding backward in time until the last pulse (the one from the direct path) is reached. A two-part strategy is also adopted in [14]. In the first part, a coarse TOA estimate is derived from signal energy measurements. Then, a hypothesis testing approach is used to refine the estimate. As in [13], pulse distortions are ignored. The approach in [1], [15]-[16] is based on the idea of a “noisy template.” No assumption is made on the pulse shapes and distortions are implicitly taken into account. The template is derived from delay-and-average operations on the received signal r(t). Correlations with r(t) are taken on suitable windows at symbol distance from each other. Next, the correlations are squared and their arithmetic mean is computed. TOA is estimated as the windows’ position that corresponds to the maximum mean. TOA estimation based on signal energy measurements is pursued in [17]. The incoming signal is squared and integrated over intervals comparable with the pulse width. The location of the DPP is computed as the index of the first interval where the energy overcomes a suitable threshold. Assessing pros and cons of the above methods is a complex task that involves performance comparisons [18]-[19] and implementation considerations. In most applications the goal is to achieve ranging accuracies of a meter with low-power and low-complexity algorithms. When dealing with complexity a common warning is that a digital approach involves analogto-digital (ADC) converters operating in the multi-GHz range. Their usage in low-power operations is problematic, especially if high-bit resolutions are required. On the other hand, an analog approach may be as challenging. For example, the algorithm in [1] requires delaying the incoming signal with errors smaller than a nanosecond. This cannot be done in analog form as the accuracy of current chips is not sufficient [20]. In this paper we provide a structured approach to TOA estimation based on least squares (LS) techniques. Similar to many previous works, we adopt a two-step procedure which first leads to a coarse estimate and then to a finer result. To this end a monocycle with pseudo-noise (PN) polarity is transmitted every Tf seconds, so generating a periodic sequence of channel responses (CRs) at the receiver. If each CR is shorter than Tf , the incoming signal is constituted of separate CRs and the aim of the coarse estimator is to compute their positions relative to a local clock running with period Tf . The fine estimator looks for the position of DPP within each CR. We assume that the monocycles propagate through a multipath channel with unknown path gains and delays and with frequency-selective characteristics. Accordingly, the shapes of CR and DPP are not available, except for some knowledge of their duration. In particular, the CR duration is approximated 839 with the channel delay spread while the DPP support is taken equal to that of the transmitted monocycle. The coarse estimator attempts to minimize the Euclidean distance between the received signal and a train of PN modulated CRs of unknown shape and time origin. A similar problem is solved by the fine estimator except that the focus is now on the DPP arrival time. In both cases the solution involves a preprocessing of the received signal followed by suitable energy measurements. Thus, the paper gives mathematical foundation to energy-based methods proposed in literature but, at the same time, it points out the importance of the pre-processing phase to enhance the estimation accuracy. The estimation scheme is digital and operates on signal samples taken at rate 1/Ts , less than or equal to the Nyquist rate. As mentioned earlier, the value of 1/Ts as well as the ADC resolution are critical parameters in view of the power consumption. We show that good performance is achieved with sampling rates of 2 GHz or less and one-bit resolution. Performance comparisons are made with some TOA estimators proposed in literature. The rest of the paper is organized as follows. The signal model is described in the next section while the coarse estimator is illustrated in section III. Section IV deals with the fine estimator and section V with implementation issues and approximations that must be made to limit the complexity of the algorithm. In section VI simulation results are provided to assess the estimator performance and make comparisons with other schemes. Conclusions are drawn in section VII. II. S IGNAL M ODEL The proposed TOA estimator is based on the transmission of a training sequence of the type sT (t) = f −1 ∞ N  ak w(t − iTb − kTf ) (1) i=−∞ k=0 where ak = ±1 is a maximal-length PN sequence [21] of (true) length Nf , w(t) is a pulse of duration Tp , Tf is the frame (true) ) and Tb = Nf Tf is the period (much longer than Tp symbol period. The propagation occurs on an L-path channel whose response to w(t) is modelled as c(t) = L−1  pℓ (t − τℓ′ ) (2) ℓ=0 pℓ (t) being the pulse from the ℓth path and τℓ′ the associated delay. Note that the waveforms pℓ (t) have different shapes in general. With no loss of generality we assume τ0′ ≤ τ1′ ... ≤ ′ so that p0 (t) represents the DPP and τ0′ the TOA we want τL−1 to estimate. For simplicity, the abbreviated notation τ := τ0′ is used henceforth. Shifting c(t) backward by τ seconds gives h(t) := c(t + τ ), which starts in the origin and is referred to as the CR. From (2) we have L−1  pℓ (t − τℓ ) (3) h(t) = ℓ=0 with τℓ := τℓ′ − τ . Next, setting c(t) = h(t − τ ) and replacing 840 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 7, NO. 3, MARCH 2008 w(t) by c(t) in (1), yields the noise-free received signal s(t) = f −1 ∞ N  ak h(t − iTb − kTf − τ ) (4) Minimizing J with respect to h̃(t), m̃ and ε̃, and calling (m̂c , ε̂c ) the pair (m̃, ε̃) corresponding to the minimum gives the coarse TOA estimate i=−∞ k=0 τ̂c = m̂c Tf + ε̂c as a train of amplitude modulated CRs. They are assumed to be separated since in most ranging applications their duration is much shorter than Tf . Their shape is unknown while their duration is approximately equal to the channel delay spread2 . Our goal is the estimation of τ based on the observation of s(t) in the presence of noise. As s(t) is periodic of period Nf Tf , τ can only be measured within multiples of that period. In other words, we must limit τ within the range 0 ≤ τ < Nf Tf . Thus, letting ⌊x⌋ the integer part of x and writing τ as a multiple of Tf plus a fraction thereof, i.e., τ = mTf + ε (5) with m := ⌊τ /Tf ⌋, we have 0 ≤ m ≤ Nf −1 and 0 ≤ ε < Tf . In the next section we propose a TOA estimator based on a presumed duration Th of h(t). In practice Th does not coincide with the true duration and, as we shall see, the discrepancy deteriorates the estimation accuracy. In a typical residential environment, errors on the order of 20-30 ns are incurred. On the other hand, as the symbol period may be as large as 1000 ns or more [22], the algorithm still provides a good reduction of the original TOA uncertainty. (10) To proceed, we substitute (8) in (9). Discarding a factor N Nf and an additive constant independent of h̃(t), m̃, ε̃ and, finally, assuming N >> 1, we get  Th  Th [h̃(t) − rf old (t + ε̃, m̃)]2 dt − rf2 old (t + ε̃, m̃)dt J= 0 0 (11) where 1 rf old (t, m̃) := N Nf N Nf −1  a|i−m̃|Nf r(t + iTf ) (12) i=0 is a folded version of r(t). Its significance is seen letting t vary between 0 and Tf . Correspondingly, r(t + iTf ) describes the ith slot of r(t) (in the interval iTf ≤ t < (i + 1)Tf ) and rf old (t, m̃) results in the average of all the slots, each weighted with an element of the sequence a|i−m̃|Nf . As h̃(t) is varied, from (11) we see that J achieves a minimum (equal to the last term) for h̃(t) = rf old (t + ε̃, m̃) 0 ≤ t ≤ Th (13) Thus, the pair (m̂c , ε̂c ) we look for is computed as III. C OARSE E STIMATOR m̂c , ε̂c = arg max To begin, we rewrite s(t) in a more convenient form. Setting Tb = Nf Tf in (4) and combining the result with (5) yields s(t) = f −1 ∞ N  (14) ε̃ (6) Next, letting j := iNf + k + m and denoting |x|Nf := x mod Nf , one has k = |j − m|Nf so that (6) may be rearranged as ∞  a|j−m|Nf h(t − jTf − ε) (7) j=−∞ The estimation of m and ε is carried out with LS techniques. To this end we denote by h̃(t) a hypothetical realization of the CR, by Th its duration, and by m̃ and ε̃ trial values of m and ε, respectively. Next, we form the hypothetical signal s̃(t) = rf2 old (t, m̃)dt The following comments are in order: ak h[t − (iNf + k + m)Tf − ε] i=−∞ k=0 s(t) = m̃,ε̃ ε̃+T  h ∞  a|j−m̃|Nf h̃(t − jTf − ε̃) (8) j=−∞ and consider its Euclidean distance (over N symbols) from the received waveform r(t) = s(t) + n(t), a mixture of signal plus noise. We have  N Tb [r(t) − s̃(t)]2 dt (9) J= 0 2 The duration of h(t) depends on the positions of the transmit and receive antennas and is difficult to measure since the operation requires estimating pℓ (t) and τℓ′ in (2). (i) The estimator (14) exploits energy measurements on rf old (t, m̃). This is in contrast with the scheme in [17] where the measurements are directly made on r(t). The advantage of operating on rf old (t, m̃) rather than r(t) is that the noise component is reduced while the signal component keeps the same timing information. To see why let us replace r(t) by s(t) + n(t) in (12). Bearing in mind (7) we get rf old (t, m̃) = sf old (t, m̃) + nf old (t, m̃) (15) with sf old (t, m̃) = ⎤ ⎡ N Nf −1 ∞   1 ⎣ a|i−m̃|N a|i−j−m|N ⎦ h(t + jTf − ε) f f N N f i=0 j=−∞ (16) nf old (t, m̃) = 1 N Nf N Nf −1  i=0 a|i−m̃|N n(t + iTf ) (17) f For example, if n(t) is white Gaussian noise with spectral density N0 /2, then nf old (t, m̃) is white as well but its power spectral density is reduced to N0 /(2N Nf ). As for sf old (t, m̃), let us recall that ak is a maximal-length PN sequence, which implies that the following relation holds D’AMICO et al.: ENERGY-BASED TOA ESTIMATION 841 ~ 2 sfold (t,m) ∼ ε ε Th Th(true) t 3Tf 2Tf Tf 0 Fig. 1. Waveform s2f old (t, m̃) for |m − m̃|Nf = 2. A channel response with nine multipath components is assumed. ε^ c Th Th(true) ε (a) t 1 in the interval 2Tf  t  3Tf ). The only chance for ε̂c to equal ε is that the strongest multipath components are packed in the first Th seconds of the CR, which is not necessarily (true) true. Finally, for Th > Th it is seen from Fig. 2b that (true) the integral in (14) exhibits a plateau of length Th − Th as ε̃ is varied. In the presence of noise any point within the plateau is likely to correspond to a maximum of the integral. In conclusion, the accuracy of the estimation algorithm depends on our guess about the CR duration. As we shall see in section VI, satisfactory results are obtained taking Th equal to the channel delay spread. (iii) The maximization in (14) must be carried out numerically, letting m̃ vary in the interval [0, Nf − 1] and changing ε̃ by multiples of some step-size △ε . As the estimation accuracy is necessarily rough, it would be useless to take △ε too small as it would only increase the computational load. In the simulations described later we set △ε = Th . 3Tf 2Tf IV. F INE E STIMATOR ^ε c Th ε Th(true) (b) t 3Tf 2Tf (true) Fig. 2. (a) Th is smaller than Th and the integral in (14) achieves (true) a maximum for ε̂c > ε in general. (b) Th is larger than Th and the integral exhibits a plateau. In this section we discuss a TOA estimator with improved performance relative to the previous section. It performs a search over a tentative delay (and other parameters) letting it vary in a small interval around τ̂c = m̂c Tf + ε̂c . In this way, the information from the coarse estimator is exploited to restrict the range of the fine search. As before, we adopt an LS approach in which the Euclidean distance of the received signal r(t) to a hypothetical realization of s(t) is minimized. To begin we substitute (2) into (4) to obtain s(t) = L−1  f −1 ∞ N  ak pℓ (t − iTb − kTf − τℓ′ ) (20) ℓ=0 i=−∞ k=0 Next we write τℓ′ as a multiple of Tf plus a fractional part εℓ [21] 1 Nf Nf −1  a|i−p|N a|i−q|N = i=0 f f  τℓ′ = mℓ Tf + εℓ 1 − N1f p = q mod Nf otherwise (18) Then, for Nf >> 1 the quantity in square brackets in (16) is either unity (for j = m̃ − m mod Nf ) or nearly zero (otherwise) and sf old (t, m̃) becomes sf old (t, m̃) ≃ ∞  (21) with mℓ := ⌊τℓ′ /Tf ⌋. Finally, plugging (21) into (20) and rearranging as in (7) we get s(t) = L−1  ∞  a|j−mℓ |N pℓ (t − jTf − εℓ ) (22) f ℓ=0 j=−∞ In view of (22) we model s̃(t) as h(t + jTb + |m̃ − m|Nf Tf − ε) (19) j=−∞ i.e., a periodic train of h(t)-pulses of period Tb . In particular, the pulse in the interval 0 ≤ t ≤ Tb starts at t = |m − m̃|Nf Tf + ε. Figure 1 illustrates s2f old (t, m̃) for |m − m̃|Nf = 2. The individual components of the CR (see (3)) are shaped as the second derivative of a Gaussian function. For simplicity, only nine multipath components are shown. (ii) A physical interpretation of (14) is useful to explain where the limitations of the TOA estimator are. Assume for (true) the true simplicity that the noise is negligible and call Th duration of h(t), as distinguished from the assumed duration (true) Th . For Th = Th it is clear from Fig. 1 that the integral in (14) achieves a maximum for m̂c = m and ε̂c = ε, which is (true) we have ε̂c > what we expected. Vice versa, for Th < Th ε in general, as illustrated in Fig. 2a (an enlargement of Fig. s̃(t) = L−1  ∞  a|j−m̃ℓ |N p̃ℓ (t − jTf − ε̃ℓ ) ℓ=0 j=−∞ (23) f where, as usual, x̃ denotes a tentative value of x. For simplicity, we take the waveforms p̃ℓ (t) with a common duration Tp . In the absence of channel distortions we would set Tp equal (true) , the monocycle duration. Next we compute the to Tp Euclidean distance of s̃(t) from r(t) as given in (9). In doing so we make the simplifying assumption that the multipath components in (3) are separated, which allows us to write p̃ℓ1 (t − τ̃ℓ′1 )p̃ℓ2 (t − τ̃ℓ′2 ) = 0 for ℓ1 = ℓ2 (24) Condition (24) may not be valid in general and serves only to facilitate the derivation. The performance of the resulting estimator in realistic circumstances is assessed later by simulation. 842 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 7, NO. 3, MARCH 2008 Under the above hypotheses, we get after some manipulations (in which a factor N Nf and an irrelevant additive constant are dropped) J= − Tp L−1  ℓ=0 0 Tp L−1  a⏐i–m⏐ ~ Fig. 3. rf2 old (t + ε̃ℓ , m̃ℓ )dt p̃ℓ (t) = rf old (t + ε̃ℓ , m̃ℓ ) 0ℓL−1 (26) and its further minimum as a function of ε̃ℓ and m̃ℓ occurs for rf2 old (t, m̃ℓ )dt 0ℓL−1 ε̃ℓ (27) The following remarks are useful: (i) Equation (27) indicates that L separate searches are needed to compute the pairs (m̂ℓ , ε̂ℓ ) and, ultimately, the TOAs τ̂ℓ′ := m̂ℓ Tf + ε̂ℓ of all the arriving pulses. The operations can be carried out sequentially by repeating L times the basic algorithm m̂, ε̂ = arg max m̃,ε̃ ε̃+T  p ∧ 2 τc I&D FINE SEARCH Nf ∧ τ Conceptual block diagram of the energy-based estimator. (25) 0 ≤ t ≤ Tp ε̃ℓ+Tp COARSE SEARCH (•) FOLDING 2 Clearly, the minimum of J as a function of p̃ℓ (t) is obtained for m̃ℓ ,ε̃ℓ r(t) [rf old (t + ε̃ℓ , m̃ℓ ) − p̃ℓ (t)] dt ℓ=0 0 m̂ℓ , ε̂ℓ = arg max ~ rfold (t,m) rf2 old (t, m̃)dt be reached by the DPP peak (missed detection). The problem is further exacerbated by the strength of the DPP, which varies in a random way with the channel realizations. This makes it difficult to design λ so as to neglect missed detections. The choice of λ is discussed in the next section. (iv) Another important issue is the search region. Its size should be as small as possible to limit the computational complexity. Although the coarse search indicates τ̂c as the arrival time of the strongest multipath components, the DPP might be in advance of them. In particular, with the channel models CM1 (residential LOS) described in [23], the DPP may arrive up to 40 ns earlier [17]. Thus, bearing in mind that τ̂c may be in error up to 30 ns or so, in the following the search region is set as τ̂c − 70 ns  τ̃ ′  τ̂c + 30 ns. (v) Figure 3 gives a conceptual scheme of the algorithm which, henceforth, is referred to as energy-based estimator (EBE). The integrals in (14) and (29) are computed in the integrate-and-dump (I&D) circuit. Some implementation issues are discussed in the next section. (28) ε̃ ′ In other words, letting τ̃ =: m̃Tf + ε̃, the L highest peaks of the integral ε̃+T  p ′ rf2 old (t, m̃)dt (29) I(τ̃ ) := ε̃ provide the TOAs of all the multipath components while the earliest peak gives the TOA of the DPP. It is worth noting that (28) coincides with (14) except for the replacement of Th with Tp . (ii) Unfortunately, as L is unknown, the above procedure cannot be applied. Indeed, if we assume L̂ paths instead of L, two adverse situations may occur. For L̂ < L, the set of the highest peaks may not include the DPP peak, meaning that a later pulse is mistaken for the DPP. Vice versa, for L̂ > L, the set of the highest peaks include noise artifacts and, again, an error occurs if an artifact precedes the DPP peak. Note that the highest peak need not be due to the DPP as the latter may be attenuated by obstructions. (iii) A way to overcome this problem is to adopt a thresholdbased approach in which the DPP TOA is computed as the first time I(τ̃ ′ ) overcomes a suitable threshold λ. Clearly, the threshold value is an important design parameter. If it is too low, i.e., comparable with the peaks I(τ̃ ′ ) exhibits in the noiseonly region (NOR) (i.e., where no signal pulses are present, as shown in Fig. 1), the threshold crossing is likely to be a noise artifact (false alarm). Vice versa, if it is too high, it may not V. I MPLEMENTATION I SSUES We start with the folding operation in (12) which involves delayed versions of r(t). It is easily recognized that the implementation of the delays must be very accurate, say within a fraction of the monocycle width, otherwise pulses from adjacent frames would not add coherently. As this accuracy cannot be achieved with analog delay lines, we must proceed digitally. Thus, we compute the samples of rf old (t, m̃) from those of r(t). Calling 1/Ts the sampling rate, assuming Tf /Ts an integer Mf , and letting r[k] := r(kTs ) and rf old [k, m̃] := rf old (kTs , m̃), from (12) we have rf old [k, m̃] = 1 N Nf N Nf −1  a|i−m̃|N r [k + iMf ] (30) f i=0 Next, we concentrate on the coarse search (14). As pointed out in section III, the variable ε̃ is varied by multiples of Th . Thus, assuming Th /Ts an integer Mh , letting ε̃ = ñMh Ts (ñ = 0, 1, 2, ..), and approximating the integral by a sum yields (ñ+1)Mh −1 m̂c , n̂c = arg max m̃,ñ  rf2 old [k, m̃] (31) k=ñMh Note that ε̂c in (14) is computed as ε̂c = n̂c Mh Ts . A similar approach is adopted with the integral in (29). Here we let ε̃ vary by multiples of Ts , i.e., we set ε̃ = ñTs (ñ = 0, 1, 2, ..). Also, we assume Tp /Ts an integer Mp . Finally, D’AMICO et al.: ENERGY-BASED TOA ESTIMATION 843 ~ I(τ´) (M p−1)Ts threshold λ ~ τ´ ~ τ´ first Fig. 4. Typical shape of peak location estimate PF A = exp − I(τ̃ ′ ). we approximate the integral by a sum. Dropping an irrelevant factor Ts we get ñ+Mp −1 I(τ̃ ′ ) =  rf2 old [k, m̃] (32) k=ñ with τ̃ ′ = (ñ + m̃Mf )Ts . Another issue is the location of the first peak of I(τ̃ ′ ) above a given threshold. To conceive a suitable procedure it is convenient to temporarily neglect the noise and think of rf2 old [k, m̃] as a sequence of zeros (before the DPP) followed by a group of Mp positive values (during the DPP). In these conditions I(τ̃ ′ ) is zero before the DPP arrival. As the DPP comes in, the first nonzero sample rf2 old [k, m̃] enters the sum in (32) and I(τ̃ ′ ) takes a positive value, as indicated in Fig. 4. In the next Mp − 1 steps I(τ̃ ′ ) increases, then it decreases toward zero. The above description suggests the following peak location strategy. Let τ̃f′ irst be the first time I(τ̃ ′ ) overcomes a given threshold λ (see Fig. 4). If I(τ̃ ′ ) is still above λ in the next Mp − 1 steps, i.e., I(τ̃f′ irst + iTs ) > λ i = 1, 2, ..., Mp − 1 (33) then the maximum of I(τ̃ ′ ) is sought î = arg nor ε are exactly known and we must content ourselves with an approximate delimitation of NOR obtained from (35) by (true) replacing Th and ε with Th and ε̂c , respectively. It is worth noting that in NOR the samples r[k] are independent zero mean Gaussian random variables (RVs) whose variance σ 2 can be easily measured. Thus, PF A can be expressed as a function of σ 2 as follows. From (30) and (32) it can be shown that I(τ̃ ′ ) is a chi-square RV with Mp degrees of freedom, mean value Mp σ 2 /(N Nf ) and variance 2Mp (σ 2 /(N Nf ))2 . Thus, for Mp even we have [24] max i=0,1,..,Mp −1 I(τ̃f′ irst + iTs ) (34) and the peak location is estimated as τ̂ ′ = τ̃f′ irst + îTs − Ts /2. Otherwise, if (33) does not hold, the original overcoming of λ at τ̃f′ irst is viewed as a noise artifact and is ignored. The threshold value has an impact on the estimation performance. We have found by simulation that satisfactory results are obtained choosing λ such that the false alarm probability PF A is about 10−8 . To see how this condition is achieved let us consider the samples r[k] in the NOR. As is now argued, such a region is the set of intervals (true) 0  t < ε and iTf + ε + Th  t  (i + 1)Tf + ε i = 0, 1, 2, ... (35) This is recognized bearing in mind that the channel re(true) sponses have a duration Th and, as indicated in (7), they (true) begin at iTf + ε (i = 0, 1, 2, ..). In practice neither Th λN Nf 2σ 2 Mp /2−1  m=0 1 m! λN Nf 2σ 2 m (36) from which λ is computed setting the right hand side equal to 10−8 . In practice the map σ 2 → λ can be stored in a readonly memory and the threshold value can be picked up as a function of the noise level. In the above discussion an infinite ADC resolution has been implicitly considered. To assess the impact of this assumption, in the next section we also report on simulations with onebit quantization. In this case the samples r[k] take values ±1 with the same probability and are independent. Then, assuming N Nf >> 1 and invoking the central limit theorem in (30), rf old [k, m̃] is approximated as a zero-mean Gaussian RV with variance 1/N Nf . Thus, I(τ̃ ′ ) is still chi-square with Mp degrees of freedom, but its mean value is Mp /(N Nf ) and its variance 2Mp /(N Nf )2 . In conclusion PF A is still expressed as in (36) except that σ 2 = 1. Sampling rates in the range 0.5 to 8 GHz, as considered in the next section, are so high that one wonders whether they can be achieved with current technology. Flash converters are the architecture of choice in these cases and low-power consumption is a constraint. The state-of-the-art for ADCs is surveyed in [25] which is dated 1999. More recent results are found in [26], [27] and references therein. In particular, [26] reports on a 6 bit converter with 1.25 GHz sampling rate and 2.5 mW power consumption, while in [27] the sampling rate is 4 GHz, the resolution is 4 bit and the consumption 89 mW. Bearing in mind that higher sampling rates can be traded with lower resolutions [25], we expect that the whole sampling range of interest can be covered at least with 1 bit ADCs. VI. S IMULATION R ESULTS Computer simulations have been run to assess the EBE performance and make comparisons with other algorithms available in literature. The following setting has been chosen. The channel response is generated according to the CM1 model in [23]. Here, the rays arrive in clusters with a Poisson distribution of rate Λ = 0.047 ns−1 . The arrival times of the returns within each cluster are a blend of two Poisson processes with rates λ1 = 1.54 ns−1 and λ2 = 0.15 ns−1 and mixture probability β = 0.095. The amplitude αk,l of the kth ray in the lth cluster follows a Nakagami-m distribution pαk,l (α) = 2 Γ(mk,l ) mk,l Ωk,l mk,l α2mk,l −1 exp − mk,l 2 α Ωk,l (37) 844 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 7, NO. 3, MARCH 2008 where Γ(mk,l ) is the gamma function and mk,l has a lognormal distribution with mean m0 = 0.67 and standard deviation m0 = 0.28. Parameter Ωk,l reflects the power decay over the clusters and within a given cluster and is given by Ωk,l = Ωl exp(−τk,l /γ) γ[βλ1 + (1 − β)λ2 + 1] (38) where γ = 12.53 ns is the intra-cluster decay constant and Ωl is the energy of the lth cluster. The latter has an exponential decay as a function of the arrival time Tl of the cluster 10 log(Ωl ) = 10 log(exp(−Tl /Γ)) + Mc (39) where Γ = 22.61 ns and Mc is normally distributed with standard deviation σc = 2.75. The assumed duration Th of the channel response equals 20 ns, about the channel delay spread. There are Nf = 15 frames in a period, each of duration Tf = 200 ns, and the symbol period is Tb = Nf Tf = 3000 ns. The TOA estimate is derived from the observation of N = 100 symbols. Three different monocycles are considered, corresponding to alternative signal spectra and receive filter characteristics. The type-1 monocycle is shaped as the second derivative of a Gaussian function and has duration 4 ns. The receive filter is low-pass and has a bandwidth of 960 MHz. Nyquist sampling in this range requires approximately 2 GHz ADC clocking. As we shall see, however, a ranging accuracy of a meter is achieved even with lower sampling rates. The type-2 monocycle is shaped as before but its duration is reduced to 1 ns. The receive filter is low-pass, with a bandwidth of 4 GHz. The pulse shortening improves the estimation accuracy but it also increases the sampling rate by a factor 4. The type-3 monocycle is a modulated raised-cosine pulse compliant with IEEE 802.15.4a standard [22]. Its shape is w(t) = p(t) cos(2πf0 t) where f0 = 3.75 GHz and p(t) = sin(πBt) cos(παBt) πBt 1 − (2αBt)2 |t|  2 ns (40) with α = 0.6 and B = 1 GHz. So, w(t) has still a duration of 4 ns but the receive filter is band-pass, with center frequency f0 and bandwidth (1 + α)B. The sampling rate varies from 1 to 3 GHz. The EBE performance with monocycles of type 1 and 2 is compared with that of the schemes discussed in [13], [1] and [17]. Some assumptions made with these schemes are now summarized. As mentioned in the introduction, the algorithm in [13] correlates the received signal with the transmitted monocycle and compares the result with a threshold. In our simulations the received pulses are taken undistorted and the threshold value is chosen (by trial and error) so as to get the best estimation accuracy at any given SNR. The sampling rate is Nyquist, i.e., 2 GHz with type-1 and 8 GHz with type-2 monocycles. The “noisy template” in [1] is computed as wN T (t) = N1 −1 1  r(t + iTb ) N1 i=0 0  t  Tb (41) where r(t) is the received signal when a sequence of all ones is transmitted. Next, wN T (t) is periodically extended with period Tb ∞  wN T (t − jTb ). (42) wN T (t) := j=0 Then, the sequence pattern is changed to symbols of alternating sign and the statistics yk (τ̃ ) := τ̃ +kT  b +TR r(t)wN T (t)dt k = 0, 1, . . . , N2 − 1 τ̃ +kTb (43) are computed, with TR a design parameter (less than Tb ) that is derived via channel sounding. Finally, the arithmetic average of yk2 (τ̃ ) is calculated over N2 symbols, say z(τ̃ ) = N2 −1 1  yk2 (τ̃ ) N2 (44) k=0 and the τ̃ corresponding to its maximum is chosen as the TOA estimate. In our simulations we have set N1 = N2 = 50, which makes the overall observation interval N1 + N2 equal to 100 symbols. The estimation accuracy varies with TR . Experimentally we have found that the best results are obtained with TR = Tb − 170 ns. The noisy template and the statistics in (43) are computed from samples of r(t) taken at Nyquist rate. The TOA estimator in [17] assumes that symbol timing is already available and confines the delay τ (see (4)) within 0  τ < Tf . This interval is divided into segments of size ∆, say (n − 1)∆  τ < n∆ (n = 1, 2, ...), and the energy of r(t) on each segment is measured. The energy from the nth interval of the kth frame of the ith symbol is given by iTb +kT  f +n∆ z[n, k, i] = r2 (t)dt (45) iTb +kTf +(n−1)∆ Next, energies with the same index n are taken from N symbols and are put together to produce Z[n] = f −1 N −1 N  i=0 z[n, k, i] (46) k=0 Finally, Z[n] is compared with a threshold λ and the index n̂ of the first crossing is computed. The TOA is estimated as τ̂ = (n̂−0.5)∆. The threshold λ is a function of the minimum and maximum values of Z[n]: λ = min {Z[n]} + K[max {Z[n]} − min {Z[n]}] (47) and the parameter K is experimentally chosen so as to optimize the estimator performance. In our simulations we have set N =100, ∆ = Tp and K equal to 0.1 or 0.2. Figures 5-6 illustrate the EBE performance and make comparisons with [13], [1] and [17]. The ordinates represent the TOA root-mean-square error (RMSE) in meters (i.e., the RMSE in seconds multiplied by the speed of light) while the abscissas give the ratio of the received signal energy per bit, Eb , to the noise spectral density N0 . A type-1 monocycle is adopted and the assumed DPP duration is Tp = 4 ns. The D’AMICO et al.: ENERGY-BASED TOA ESTIMATION 845 Eb/N0, dB Eb/N0, dB Fig. 5. Comparisons between EBE and the algorithms in [13] and [1] with type-1 monocycle and infinite ADC resolution. Fig. 7. Comparisons between EBE and the algorithms in [13] and [1] with type-2 monocycle and infinite ADC resolution. Eb/N0, dB Eb/N0, dB Fig. 6. Comparisons between EBE and the algorithm in [17] with type-1 monocycle and infinite ADC resolution. Fig. 8. Comparisons between EBE and the algorithm in [17] with type-2 monocycle and infinite ADC resolution. ADC has infinite resolution and the sampling rate fs varies from 0.5 GHz to 2 GHz. As mentioned before, the sampling rate with [13] and [1] is 2 GHz while the I&D circuit that computes the integrals (45) is clocked at 1/∆ = 250 MHz. It is seen that EBE outperforms the algorithms in [1] and [17] at any Eb /N0 . For fs = 2 GHz the method in [13] is a little better at low SNRs but its advantage is questionable in view of the optimistic assumptions made on the template and threshold setting. Figures 7-8 show analogous results with a type-2 monocycle. Everything is as before except that the pulse duration is reduced to 1 ns. As expected, all the algorithms have improved performance but the ranking remains the same. Degradations due to mismatch between true and assumed pulse duration are discussed in Fig. 9. It is seen that, while taking Tp < Tp has only marginal consequences, the op(true) posite is true with Tp > Tp , especially at high SNR. This suggests a strict attitude when guessing the pulse duration. Figure 10 illustrates the effects of a finite ADC resolution. The monocycle is of type 1 and the parameters are as in Fig. 5. Dashed lines indicate the EBE performance with a 1-bit ADC resolution. As we see, the losses due to quantization are not severe. Further simulations (not shown for reasons of space) show that they reduce to a fraction of dB with a 2-bit resolution. Finally, Fig. 11 shows the EBE performance with type-3 monocycle. The sampling rate varies from 1 to 3 GHz. Again, the quantization to 1-bit has a limited impact on the estimation accuracy. Interesting enough, with a sampling rate of 2 GHz the performance at very low and very high SNR is about as (true) 846 IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS, VOL. 7, NO. 3, MARCH 2008 Eb/N0, dB Fig. 9. Effects of mismatch between true and assumed pulse duration. (true) Sampling rate fs = 2GHz and pulse duration Tp = 4 ns. Eb/N0, dB Fig. 11. EBE performance with type-3 monocycle. Several approximations to the original scheme have been made to facilitate implementation. In particular, the effects of sub-Nyquist sampling have been explored along with the impact of a reduced ADC resolution. It has been shown that estimation accuracies of a meter are achievable even with sampling rates of 500 MHz and 1-bit resolution. Comparisons with other estimation methods available in literature show that, with comparable implementation complexity, our scheme exhibits much better performance. ACKNOWLEDGMENTS The authors would like to express their appreciation to the anonymous reviewers for their comments and suggestions that have improved the quality of the paper. Eb/N0, dB Fig. 10. Effects of finite ADC resolution. in Fig. 10, although monocycles of the same duration but different form are used. This suggests a weak dependence of the estimation performance on the monocycle shape. VII. C ONCLUSIONS A new TOA estimation algorithm has been proposed that operates on energy measurements made on a folded version of the received signal. Its motivation is based on the idea that, because of frequency-selective propagation effects, the shape of the received pulses is not known at the receiver and, in consequence, a correlation-based estimator is not feasible. We have revisited the TOA estimation problem assuming only an approximate knowledge of the pulse duration. Making use of LS techniques we have come up with an algorithm that exploits the pulse duration information as a window size in energy measurements. R EFERENCES [1] S. Gezici, Z. Tian, G. B. Giannakis, H. Kobayashi, A. F. Molish, H. V. Poor, and Z. Sahinoglu, “Localization via ultra-wideband radios,” IEEE Signal Processing Mag., vol. 22, no. 4, pp. 70–84, July 2005. [2] S. M. Kay, Fundamentals of Statistical Signal Processing: Detection Theory. Upper Saddle River, NJ: Prentice-Hall, Inc., 1998. [3] M. Z. Win and R. A. Scholtz, “Characterization of ultra-wide bandwidth wireless communications channels: A communication theoretic view,” IEEE J. Select. Areas Commun., vol. 20, no. 9, pp. 1613–1627, Dec. 2002. [4] R. C. Qiu, C. Zhou, and Q. Liu, “Physics-based pulse distortion for ultra-wideband signals,” IEEE Trans. Veh. Technol., vol. 54, no. 5, pp. 1546–1555, Sept. 2005. [5] R. A. Scholtz, D. M. Pozar, and W. Namgoong, “Ultra-wideband radio,” EURASIP J. Applied Signal Processing, vol. 2005, no. 3, pp. 252–272, March 2005. [6] A. F. Molisch, “Ultrawideband propagation channels-theory, measurements, and modeling,” IEEE Trans. Veh. Technol., vol. 54, no. 5, pp. 1528–1545, Sept. 2005. [7] S. Gezici, H. Kobayashi, and H. V. Poor, “Non-parametric non-lineof-sight identification,” in Proc. 58th Vehicular Technology Conference (VTC 2003 Fall), vol. 4, Oct. 2003, pp. 2544–2548. [8] S. Venkatraman and J. Caffery, “A statistical approach to non-lineof sight BS identification,” in Proc. 5th International Symposium on Wireless Personal Multimedia Communications (WPMC 2002), vol. 1, Oct. 2002, pp. 296–300. D’AMICO et al.: ENERGY-BASED TOA ESTIMATION [9] A. Taha and K. M. Chugg, “On designing the optimal template waveform for UWB impulse radio in the presence of multipath,” in Proc. IEEE Conf. on Ultra Wideband Systems and Technologies, May 2002, pp. 41–45. [10] R. D. Wilson and R. A. Scholtz, “Template estimation in ultra-wideband radio,” in Proc. of the 37th Asilomar Conf. on Signals, Systems and Computers, vol. 2, Nov. 2003, pp. 1244–1248. [11] E. A. Homier and R. A. Scholtz, “Rapid acquisition of ultra-wideband signals in the dense multipath channel,” in Proc. IEEE Conf. on UltraWideband Systems and Technologies, May 2002, pp. 105–109. [12] ——, “A generalized signal flow graph approach for hybrid acquisition of ultra-wideband signals,” International J. Wireless Inf. Networks, vol. 10, no. 4, pp. 179–191, Oct. 2003. [13] J.-Y. Lee and R. A. Scholtz, “Ranging in a dense multipath environment using an uwb radio link,” IEEE Trans. Select. Areas Commun., vol. 20, no. 9, pp. 1677–1683, Dec. 2002. [14] S. Gezici, Z. Sahinoglu, A. F. Molish, H. Kobayashi, and H. V. Poor, “A two-step time of arrival estimation algorithm for impulse radio ultra wideband systems,” in Proc. 13th European Signal Processing Conference (EUSIPCO), Sept. 2005. [15] L. Yang and G. B. Giannakis, “Ultra-wideband communications: An idea whose time has come,” IEEE Signal Processing Mag., vol. 21, no. 6, pp. 26–54, Nov. 2004. [16] ——, “Timing ultra wideband signals with dirty templates,” IEEE Trans. Commun., vol. 53, no. 11, pp. 1952–1963, Nov. 2005. [17] I. Guvenc and Z. Sahinoglu, “Threshold-based TOA estimation for impulse radio UWB systems,” in IEEE Int. Conf. on UWB, Sept. 2005, pp. 420–425. [18] D. Dardari and M. Z. Win, “Threshold-based time-of-arrival estimators in UWB dense multipath channels,” in Proc. IEEE Int. Conf. on Communications, June 2006. [19] D. Dardari, C.-C. Chong, and M. Z. Win, “Analysis of threshold-based TOA estimators in UWB channels,” in Proc. 14th European Signal Processing Conference (EUSIPCO), Sept. 2006. [20] “Active delay lines,” RCD Components Inc., Tech. Rep., [online report], URL: http://www.rcd-comp.com/rcd/rcdpdf/A08-A14-SA08SMA14.pdf [cited 6 July 2006]. [21] J. J. Stiffler, Theory of Synchronous Communications. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1971. [22] G. M. Maggio, “TG4a UWB-PHY layer overview,” IEEE P802.15 Working Group for Wireless Personal Area Network, IEEE P802.1505/707r1-TG4a, Tech. Rep., November 2005. [Online]. Available: http://www.802wirelessworld.com [23] A. F. Molisch, K. Balakrishnan, C.-C. Chong, S. Emami, A. Fort, J. Karedal, H. Schantz, U. Schuster, and K. Siwiak, “IEEE 802.15.4a channel model-final report,” Tech. Rep., Sept. 2004, [online report], URL: http://www.ieee.org/15/pub/TG4a.html, document 04/535r0. [24] M. K. Simon, Probability Distributions Involving Gaussian Random Variables. Boston: Kluwer Academic Publishers, 2002. [25] R. H. Walden, “Analog-to-digital converter survey and analysis,” IEEE J. Select. Areas Commun., vol. 17, no. 4, pp. 539–550, April 1999. 847 [26] G. Van der Plas, S. Decoutere, and S. Donnay, “A 0.16pJ/conversionstep 2.5 mW 1.25 GS/s 4 bit ADC in a 90 nm digital CMOS process,” in Proc. 2006 IEEE International Solid-State Circuits Conference, Feb. 2006. [27] S. Park, Y. Palaskas, and M. P. Flynn, “A 4 GS/s flash ADC in 0.18 µm CMOS,” in Proc. 2006 IEEE International Solid-State Circuits Conference, Feb. 2006. Antonio A. D’Amico received the Dr. Ing. Degree in Electronic Engineering in 1992 and the Ph. D. degree in 1997, both from the University of Pisa, Italy. He is currently a Research Fellow at the Department of Information Engineering of the University of Pisa. His research interests are in digital communication theory, with emphasis on synchronization algorithms, channel estimation and detection techniques. Umberto Mengali (M69-SM85-F90) received his education in Electrical Engineering from the University of Pisa. In 1971 he obtained the Libera Docenza in Telecommunications from the Italian Education Ministry. Since 1963 he has been with the Department of Information Engineering of the University of Pisa where he is a Professor of Telecommunications. In 1994 he was a Visiting Professor at the University of Canterbury, New Zealand as an Erskine Fellow. He has served in the technical program committee of several international conferences and has been technical Co-Chair of the 2004 International Symposium on Information Theory and Applications (ISITA). His research interests are in digital communications and communication theory, with emphasis on synchronization methods and modulation techniques. He has co-authored the book Synchronization Techniques for Digital Receivers (Plenum Press, 1997). Professor Mengali is a member of the Communication Theory Committee and has been an editor of the IEEE Transactions on Communications from 1985 to 1991 and of the European Transactions on Telecommunications from 1997 to 2000. Lorenzo Taponecco was born in Sarzana, Italy, in 1977. He received the Dr. Ing. degree (cum laude) in telecommunications engineering from the University of Pisa, Pisa, Italy, in 2002, where he is currently working towards the Ph.D. degree. Since 2003, he has been with the Department of Information Engineering at the University of Pisa. His research interests include the area of digital communication theory, with special emphasis on ultrawideband (UWB) communications, parameter estimation and synchronization techniques.