Joint Ranging and Clock Parameter Estimation by Wireless Round Trip Time Measurements

3.6 MB·English

by  Satyam Dwivedi

Joint Ranging and Clock Parameter Estimation by
Wireless Round Trip Time Measurements
Satyam Dwivedi, Alessio De Angelis, Dave Zachariah, Peter Ha¨ndel
Abstract—In this paper we develop a new technique for estimators to precisely estimate range and clock parameters
estimating fine clock errors and range between two nodes fromthemeasurements.Theestimatedrangecansubsequently
simultaneously by two-way time-of-arrival measurements us-
be used for positioning using standard techniques [12]–[14].
ing impulse-radio ultra-wideband signals. Estimators for clock
Such an approach would simplify physical infrastructure and
parameters and the range are proposed that are robust with
5
respecttooutliers.Theyareanalyzednumericallyandbymeans obviate the need for e.g. laying cables for synchronization in
1
of experimental measurement campaigns. The technique and positioning systems.
0
2 derived estimators achieve accuracies below 1Hz for frequency Similarly, wireless clock synchronization can be useful in
estimation, below 1ns for phase estimation and 20cm for severalotherapplications;whicheitherpresentlyrelyonwired
n range estimation, at 4m distance using 100MHz clocks at both
synchronization or work asynchronously but exhibit perfor-
a nodes. Therefore, we show that the proposed joint approach is
J practical and can simultaneously provide clock synchronization mance loss in the absence of synchronization. For instance,

2 and positioning in an experimental system. many packet-based radio systems operate asynchronously and
2 Keywords: clock synchronization, range estimation, robust esti- suffer packet collisions. Joint positioning and time synchro-
mators, ultra-wideband nization can help reduce such losses in a network. Tradi-
]
I tionally, time synchronization is a necessity for a variety
N I. INTRODUCTION
of networks. Network time protocol (NTP) is a well-known
. A common time reference among nodes along with node mechanism to synchronize clocks over the internet [15].
s
c positioning enables coordination in time and space domains. Similarly, time synchronization protocol for sensor networks
[ Bysuchcoordinationthenodescanexecutetasksefficientlyin (TPSN) is in use [16]. In the reference broadcast synchro-

1 internet-of-things(IOT),machine-to-machine(M2M)commu- nization (RBS) protocol, time synchronization is achieved by
v nication,orsimilarscenarios[1]–[3].Incertainscenariostime broadcast beacon and receiver timestamp exchanges [17]. In
0 synchronization among nodes is necessary for positioning. all these network synchronization mechanisms, the clocks are
5
Specifically, it is needed in positioning techniques where time eventually synchronized between two nodes in a network. In
4
ofarrival(TOA)measurementsarerequiredandmeasurements this paper, by contrast, we avoid timestamp exchanges, and
5
0 are performed using a clock. As an example, in the commer- hence communication overhead, using round-trip time (RTT)
. cially available Ubisense system, sensor nodes are synchro- measurements. This enables simultaneous estimation of the
1
nizedusingcableinordertoperformtimedifferenceofarrival clock parameters and range between two nodes.
0
5 (TDOA)measurementsofultra-wideband(UWB)pulsesfrom To achieve high accuracy in wireless clock synchronization
1 tag nodes [4]. In these types of systems, replacing the cabling it is necessary to estimate the delay arising from the signal
v: with wireless synchronization will reduce the installation cost propagation time between nodes. Here we consider using
i instationaryenvironmentsandincreaseinflexibilitywillopen UWBasthewirelesstechnologyanddemonstrateresultsusing
X
up new mobile applications. Positioning and effects of clock a UWB-based radio with a very accurate time measuring
r parametermismatchhavebeenaddressedjointlybeforein[5]– device or chronometer. Our previous work has demonstrated
a
[7]andthereferencestherein.In[5],optimalalgorithmswere time measurement precision up to 50ps using a time-to-
derived for joint positioning and synchronization with anchor digital converter (TDC) [18]. Furthermore, usage of UWB in
uncertainty in position and time. In [6], the authors proposed very accurate range and position estimation is already well
distributed algorithms for positioning and synchronization in known [6], [12], [13], [19]. Under line of sight conditions
UWB ad-hoc networks. In [7], the effect of clock skew was the signal propagation delay is directly proportional to the
discussed and estimators for TDOA-based positioning were distance between the nodes. Estimating this delay therefore
formulated. A few recent works on joint ranging and clock enables ranging as a by product of clock synchronization.
synchronizationarementionedin[8]–[11].In[9],authorshave Very few clock parameter estimation results have been
proposed a network wide clock synchronization and ranging verifiedexperimentallyintheliterature.In[17],measurements
mechanism by pairwise timestamp measurements aided with were captured using Berkeley motes and logic analyzers and
passive listening. In [10], authors have proposed clock syn- verified offline by synchronization algorithms. The timescale
chronization and localization along with a low cost testbed and accuracy in these experiments are in microseconds.
for its demonstration. Whereas,inourproposedsystemthetimescaleandaccuracyis
Ourmaincontributioninthispaperliesinsuggestingamea- intheorderofanano-second,asshowninsubsequentsections.
surement mechanism which provides information about clock In [20], a timestamping system for temporal information
parametersandrangebetweennodes.Further,weproposejoint disseminationisproposedachievingsub-nanosecondaccuracy
with UWB pulses. The timestamping system is integrated in
a ZigBee platform for wireless sensor networks. However, in Ping Slave
Master
[20], one-way time-of-arrival (TOA) is considered, and range Respond
estimation is not supported, differently from our proposed
technique. Transmitter Receiver
In the literature various theoretical models and methods Transmitter Receiver
δ
have been discussed for clock synchronization over wire-
slave
less networks, cf. [21]–[23]. In [22], authors have suggested
t
message-passing methods for network clock estimation. The
Start Stop
relative clock model between two clocks adopted in next
section is discussed in [23] in greater detail. In [24], the TimetoDigitalConverter master
(TDC) RTT
Crame´r-Rao bounds were derived on the parameters for this
clockmodel.Insubsequentsectionswewillproposetheusage Fig.1. Principleofoperationoftheconsideredsystem:themasternodesends
of a clocked delay in the slave node. In our previous work, aPINGsignalandtheslavenoderesponds.Themasternodeisequippedwith
atime-to-digitalconvertercapableofmeasuringtheround-trip-time.
the generated delay is assumed to be fixed and analog [14],
[25], [26]. Here, by contrast, we turn the clock into a delay-
generating device producing informative measurements.
PING signal and the reception of a RESPOND signal, sent by
Since the proposed technique estimates the range between
the slave, after a predetermined delay. A conceptual overview
a master and a slave node, it can be used as a fundamental
is given in Fig. 1. This time interval is denoted as the round-
building block of scalable and cooperative distance-based
trip-time (RTT) which we express as a function y˜ of the
techniques for network positioning. Specifically, an example
unknown range between the master and slave ρ and the delay
ofawidely-usedtechniqueexploitingnode-to-noderangingis
δ
basedonmulti-dimensionalscaling(MDS)[27].Furthermore, 2
y˜=δ+ ρ, (1)
cooperative message-passing methods based on Bayesian net- c
works have been developed in [12], and experimentally eval-
wherecisthespeedoflightinmeterspersecond.Convention-
uated on a UWB ranging platform in [28]. An analysis and
ally, the delay is an ideal analog implementation and known
implementationofsuchmethodsisbeyondofthescopeofthe
at the master, cf. [14], [18]. Here, however, we assume that
present paper.
all timed events are based on local clocks at each node. As
Theremainderofthepaperisorganizedasfollows:Section
we will show below, this enables the joint estimation of both
II describes the considered system model for the paper. In
range and parameters for clock synchronization.
section II-B, the clock parameter measurement mechanism
Time is measured using a clock by counting a number of
is introduced. This section also introduces the measurement
periods of the clock signal
model. In section III, estimators are proposed for estimating
n
clock parameters and range from the measurements. Section C = +ϕ
0
f
IV is composed of a numerical evaluation of the model 0
proposed in section II. Section V provides results from exper- where n is the number of clock cycles, f is the nominal
0
imental setup and performance of estimators on the acquired clock frequency and ϕ [0,1/f ) is the initial phase of the
0
∈
data. In section VI we conclude the paper. clock when the measurement begins. Any fractional deviation
Notation: 1 denotes a column vector of 1s. a b is from the nominal clock frequency is termed as skew of the
(cid:12)
the element-wise, or Hadamard, product between vectors a clock, represented usually as α > 0. Time is measured using
and b. x1/2 is the element-wise square-root of vector x. a skewed clock, expressed in the nominal frequency f as
0
x W =√x(cid:62)WxwhereWisapositivesemidefinitematrix. n
(cid:107)Th(cid:107)e modulus function is written compactly as mody(x) (cid:44) C = αf0 +ϕ,
modulo(x,y).
and the actual frequency f deviates from nominal clock
frequency given by α = f /f. Equivalently the above model
0
II. SYSTEMMODELANDPROBLEMFORMULATION
is widely used as [7], [29], [30]:
In this section, we first describe the principle of operation
C = αt+ϕ,
oftheconsideredsystem,thenweillustratetheproposedclock
measurement mechanism. Subsequently, we define the related
For simplicity and to establish the notation, let us first
signal model and present a basic estimator of the parameters
consider a standard linear clock model of the observed times
of interest, based on existing literature.
at master and slave nodes,
C = α t+ϕ
A. Principle of Operation: Round-Trip Time m m m
Weconsidertwonodesequippedwithawirelesstransmitter and
and a receiver. One being called the master and the other
C = α t+ϕ ,
called slave [18]. We assume that the master node is capable s s s
of measuring the time interval between the transmission of a where α and α denote their clock skew with respect to
m s
δ
C 0
s
αs = 1
αm ̸ 1
ρ
c
h(t)
φ
2π(fm fd)
−
C
m y˜
Fig. 2. Relative clock error models, (Cm,Cs) in (2) αs/αm denotes the
clockskewandφinitialphaseinradians.
Fig.3. Aspace-timediagramdepictingtheclockedges(verticalarrows)at
theslave(toprow)andmaster(bottomrow).Thediagramrelatesthequantities
in(1)and(3).Specifically,theremaindertermh(t)in(3)ishighlighted.
nominal frequency f and ϕ is the initial phase offset.
0
The clock parameters are assumed time-invariant [24]. If
fm and fs denote the clock frequencies at the master and nism, which we will describe in more detail below, the time
slave, then we define the frequency difference fd (cid:44)fm fs. measuredatthemasterisgivenbyanintegernumberofslave
−
The relative clock skews will then be parameterized by fd as clockcycles(i.e.δ )plusaremainderterm,whichisnaturally
0
αs/αm =fm/fs =fm/(fm fd). The time observed at the modeled by a periodic modulus operation. In the following
−
slave Cs can then be expressed relative to a reference, i.e. the subsectionwepresentamodelofthisterm,denotedh(t;f ,φ)
d
master node clock Cm, which is dependent on the unknown clock parameters f and
d
(cid:18)C ϕ (cid:19) φ. Therefore, we can write the delay as
m m
C =α − +ϕ
s s s
α
m δ =δ (f )+h(t;f ,φ) (3)
(cid:18) (cid:19) 0 d d
α α
= s Cm+ ϕs s ϕm (2) andincombinationwith(1)theround-trip-timemeasurements
α − α
m m
(cid:18) (cid:19) will be informative of the range ρ to the slave but also of the
f φ
= m C + , relative clock parameters of the slave, f and φ, via δ.
m d
f f 2π(f f )
m d m d
− −
where the second term is the relative clock offset in seconds
B. Clock Measurement Mechanism
and we correspondingly define the relative clock offset in
radians φ as Assuming the frequency difference is small, i.e. fd fm,
(cid:28)
φ(cid:44)2π(f f )(cid:18)ϕ αs ϕ (cid:19) [0,2π). wfe cafn neagnldecthtutshewdeecpaenndapenpcroexoimf aδt0eoδnafsdasiknncoewfnsc=onfsmtan−t,
m d s m d m 0
− − α ∈ (cid:39)
m i.e.
(C(cid:48) C )f
This model is illustrated in Fig. 2. Given that Cm is the δ (cid:98) s− s m(cid:99).
0
reference clock, fm is known and to perform synchronization (cid:39) fm
of master and slave clocks, only estimation of fd and φ is Similarly, (3) can be approximated at a specific time instant t
necessary [31]. as
If we assumethat the slave generates thedelay by counting (C(cid:48) C )(f f ) 1
a fixed and known number of clock cycles, equivalent to δ = (cid:98) s− s m− d (cid:99) + mod2π(fdt+φ)
f f 2π(f f )
m d m d
quantizationinthetimedomain,thenominaldelayδ becomes − −
0 1
dependent on the frequency difference and can be written as δ0+ mod2π(fdt+φ)
(cid:39) 2πf
m
(C(cid:48) C )(f f ) (4)
δ (f )= (cid:98) s− s m− d (cid:99),
0 d
fm fd where periodic modulus operation represents the remainder
−
where Cs and Cs(cid:48) denote the time of reception of the PING andisdependentontheclockoffsetandfrequencydifference.
and transmission of the RESPOND signal at the slave. We now illustrate how the remainder term appears at the
On the other hand, the master is equipped with a TDC, master node. Consider the following example based on the
and is therefore able to measure time intervals with a much events described in Fig. 1:
higher resolution than its clock period. Based on this mecha- 1) Master sends a PING pulse on a positive edge of clock,
TDC receives the START signal instantly and starts 4984
measuring time.
2) On receiving the pulse, slave starts the delay generation 4982
from subsequent positive edge of its own clock.
3) SlavesendsRESPONSEafterthedelayδ0(say,twoclock ns]4980
cycles in this example) [
)
)
4) Master receives the RESPONSE, TDC gets the STOP t(i4978
signal and measures the RTT. y(
T
This sequence of events is also illustrated in the space-time RT4976
diagram of Fig. 3 where the remainder term h(t) is high-
lighted. The master repeats the above transmission procedure
4974
withafixedperiodicity.TheresultisillustratedinFig.4which
shows the role of the clock in repeated pings. At instance
‘5’ PING from master is received by slave at a clock edge 49720 100 200 300 400 500
and hence results in smallest RTT. Whereas, at instance ‘6’ sampleindex [i]
slave receives PING narrowly after the clock edge and hence
Fig. 5. Experimental data record showing the “sawtooth” behavior of the
waits nearly a full clock duration to start its delay count of a
RTTmeasurements.Thedatawasobtainedwiththemasterandslaveantennas
clockcycle.ThisresultsinlargestRTTmeasurementatmaster. placed0.5mapart,δ0=4900ns.
As can be seen, the clock offset and frequency results in a
perodictime-varyingremaindertermh(t)withasawtooth-like
waveform corresponding to the modulus function that follows
up to 50ps [32]. As shown in Fig. 1, TDC measures the time
from rounding.
difference between two events and thus measures RTT. The
In addition to the effect of the rounding operation, the
TDC has its own clock for coarse counting [32].
remainder term in (4) is also subject to zero-mean random
Here we wish to highlight a few important features:
clock jitter, which we denote v. At measurement time instant
t we therefore model the remainder as • Measurement model (7) arises from sampling the slave
i
clock over a wireless medium and measuring the relative
T
h(t ;f ,φ)= mmod (f t+φ+v(t )), (5) clock drift using a TDC at master node.
i d 2π d i
2π • No timestamping and hence no data transfer is required
where T =1/f . Finally, inserting (5) into (3) under f between master and slave devices.
m m d
(cid:28)
fm, (1) results in the following the RTT measurement model • Traditionally, clock parameters are estimated as shown
at time t : in Fig. 2, where measured time monotonically increases.
i
T 2 Long time measurements using a clock require high
y(ti)= mmod2π(fdti+φ+v(ti))+δ0+ ρ+n(ti), (6) accuracy. Whereas, RTT measurements as seen in Fig. 5
2π c
require short time measurements, hence clock accuracy
where n(t ) denotes zero-mean measurement noise. Figure 5
i andstabilityisalesserconcern.Whenmeasuringsmaller
shows a measurement capture of above phenomena on our
time duration the TDC clock skew is assumed to be
flexible UWB radio test-bed [18]. In Section V-E we validate
negligible.
the model (6) by means of residual analysis, justifying the
• Clock parameter information is available in the time
approximations used above for the considered experimental
domain as well as the frequency domain. The latter fact
setup.
facilitates extracting information using frequency-based
methods.
C. Signal Model and Measurements
Let y = [y(t ) y(t )](cid:62) RN and t = [t t ](cid:62) Next, we introduce a conventional estimator for estimating
1 N 1 N
RN be the vector··o·f N meas∈urements and time··i·nstance∈s, fd, φ and ρ by unwrapping the saw-tooth measurement (7);
an early work is [29], cf. discussion in [24], [33].
respectively. Then model (6) can be written compactly in
vector form:
2ρ D. Unwrapped Least Squares (ULS)
y=h(f ,φ,v)+δ 1+ 1+n, (7)
d 0
c
We consider the model in (6) and formulate a three-step
where unwrapped least-squares method (ULS). First, by approximat-
T ing the arithmetic average of the modulus term in (6) as the
h(f ,φ,v)(cid:44) mmod (f t+φ1+v).
d 2π 2π d amplitude of the sawtooth-like waveform, Tm/2, we obtain
The vectors v and n contain the random noise terms with N
1 (cid:88)
unknown variances σ2 and σ2, respectively. The goal is to y¯= y(t )
v n N i
estimate ρ, fd and φ from y. i=1
CentraltoRTTmeasurementproceduresisthedevicecalled N
Tm 2ρ 1 (cid:88)
time-to-digital converter (TDC) as is shown in Fig. 1. Com- +δ + + n(t ).
0 i
(cid:39) 2 c N
mercially available TDCs can measure time with resolution i=1
RTT
Clock period at slave
1 2 3 4 5 6 7 8 9 10 11
Clock period at master
Fig. 4. A diagram depicting the clock mechanism. Note that the RTT measurements, plotted below, vary depending on the shift of the slave clock edge
relativetothemasterclock.
From this we form an approximate least-square estimate of periodogram [36] we can compute an estimate up to an
the range unknown sign,
(cid:18) (cid:19)2 (cid:12) (cid:12)2
ρˆ=argρmin y¯− T2m −δ0− 2cρ fˇd =argmax (cid:12)(cid:12)(cid:12)(cid:88)N y(ti)e−j2πfdti(cid:12)(cid:12)(cid:12) ,
(cid:18) (cid:19) fd∈F (cid:12)i=1 (cid:12)
c T
m
= 2 y¯− 2 −δ0 . where F [0,fmax] is a uniform grid of frequencies. Under
⊂
uniform sampling, the estimate can be computed efficiently
Then, we normalize the measurements y(ti) to the [ π,π)- using the fast Fourier transform. Note that f determines the
− d
interval as follows:
slope of the periodic sawtooth waveform in (6), cf. Fig. 5.
z(t )(cid:44) 2π (y(t ) y¯) . The sign of fd can be resolved jointly with phase estimation
i i
T − by correlating the observed signal with sawtooth waveforms
m
Next, we apply a ‘phase-unwrapping’ algorithm, denoted by based on fˇd and fˇd, corresponding to positive and negative
−
(z) where z = [z(t ) z(t )](cid:62) [24], [29], [34]. This slopes.
1 N
Ualgorithmaddsmultiplesof···2πradtothephasewhentheab- Then using the estimate fˇd a waveform p(φ,s) =
solutedifferencebetweenco±nsecutiveelementsisgreaterthan mod2π(sfˇdt + φ1) is generated, for a fixed frequency sign
π rad. Finally, we obtain φˆ and fˆ by solving the following s = 1,1 and nominal phase φ. The phase estimate, and
d {− }
linear least-squares problem on the unwrapped samples: the resolved sign of the frequency, is then obtained by the
correlation peak,
[fˆd,φˆ]=arfgdm,φin (cid:107)U(z)−fdt−φ1(cid:107)22. (8) [φˆ,sˆ]= argmax y(cid:62)p(φ,s),
| |
Whilst the ULS estimator is computationally efficient, erro- φ∈Φ,s∈{−1,1}
neously unwrapped samples by () which occur in poor where Φ [0,2π) is a uniform grid of phases. The frequency
signal conditions lead to performaUnce· degradation. estimate i⊂s fˆ = sˆfˇ. Finally, the range estimate is obtained
d d
Furthermore, as formulated it is not robust with respect to by a direct least-squares fit using (7) which results in a
noise outliers and its estimates may deteriorate significantly. particularly simple form. Using vector notation, it can be
The presence of outliers is a common issue in practical written as
wireless measurement systems. As an evaluation example, the (cid:13) (cid:13)2
wireless localization scenario is considered in [35]. In such ρˆ=argmin(cid:13)(cid:13)y Tmp(φˆ,sˆ) δ01 2ρ1(cid:13)(cid:13)
(cid:13) − 2π − − c (cid:13)
a scenario, outliers in the time-of-arrival measurements are ρ 2
(cid:18) (cid:19)
present due to propagation channel conditions, and robust re- = c 1(cid:62) y Tmp(φˆ,sˆ) δ 1 .
0
gression techniques are employed for the purpose of reducing 2N − 2π −
the adverse impact of outliers on localization performance.
We summarize this low-complexity three-step approach as
computingtheperiodogramandcorrelationpeaks(PCP).This
III. PROPOSEDESTIMATORS
approach circumvents the need for phase unwrapping the
In view of the limitations of the unwrapped least squares dataandutilizesfrequencydomaincharacteristicsofmeasure-
approachdescribedabove,wedeveloptwoalternativemethods ments.
to estimate the offset frequency f , offset phase φ and range
d
ρ in (7).
B. Weighted Least Squares (WLS)
A. Periodogram and Correlation Peaks (PCP)
In the first approach we exploit the periodicity of the In the second approach we address the issue of robust
signal model (6) to obtain a frequency estimate. Using the estimation in the presence of noise outliers, using a weighted
squared-error criterion A. Simulation Setup
J(fd,φ,ρ)(cid:44)(cid:13)(cid:13)y−h(fd,φ,0)−δ01−2c−1ρ1(cid:13)(cid:13)2W, (9) othTehrweisneuminedriiccaatleds,imuusilnatgiotnhse hfoavlleowbienegnfipxeerdfovrmaluedes, oufnltehses
and find the minimizing arguments. Here W = diag(w) parameters of interest: f = 32Hz, ρ = 2m. For each
0 RN×N is a diagonal weight matrix. The weights w ca(cid:23)n Monte Carlo iteration, thedinitia−l phase φ is set to a random
∈ +
be chosen to mitigate the effect of outliers by downweighting valueaccordingtoauniformdistributioninthe[0,2π)interval.
the samples likely to be outliers, as discussed in the next Unless otherwise noted, a record comprised of N = 100
subsection. This makes WLS robust with respect to outliers. samplesisgeneratedforeachMonteCarloiterationaccording
Setting uniform weights w 1 reduces (9) to a standard to the model in (7). The true value of the master clock
∝
least-squares criterion. frequency is T = 10−8s, the value of δ is 5 10−6s and
m 0
×
We begin by minimizing (9) with respect to ρ. This gives a measurement update period of T = 10−3s has been used.
s
Taking into account the model in (7), we define
cw(cid:62)(y h(f ,φ,0) δ 1)
ρˆ= − d − 0 . (10)
21(cid:62)w T2 (2π)2
SNR (cid:44)10log m , SNR (cid:44)10log .
Inserting the estimate (10) back into (9) and defining c 10 σ2 j 10 σ2
n v
r(f ,φ) (cid:44) y h(f ,φ,0) δ 1, we can write the mini-
d d 0 SNR is signal-to-noise ratio over the wireless channel and
− − c
mization with respect to f and φ as a two-dimensional grid
d SNR is signal to noise ratio of the clock jitter noise [39].
j
search
Thechannelnoiseforsimulationsisgeneratedinsecondsand
w(cid:62)r(f ,φ)2 thejitternoiseisthedimensionlessphasenoiseascanbeseen
[fˆ,φˆ]= argmin w1/2 r(f ,φ) 2 | d | ,
d (cid:107) (cid:12) d (cid:107)2− 1(cid:62)w from(7).Theperformancesareevaluatedusingtherootmean-
fd∈F,φ∈Φ
(11) square error (RMSE) and computed by averaging over 1000
Monte Carlo iterations.
where F [ f ,f ] and Φ [0,2π) denote the grids.
max max
⊂ − ⊂
See the appendix for the derivation.
B. Simulation Results
A complete characterization of the Monte Carlo simulation
C. Choosing the Weights for Outlier-prone Data results is provided in Fig. 6. Due to the fixed resolution of
the search grid, WLS and PCP exhibit quantization effects,
For applying the WLS method derived in Section III-B to
as shown e.g. in the asymptotic constant behavior at high
outlier prone data, we propose to choose the i-th element w
i SNR in Fig. 6a. Potentially, these effects can be overcome by
of the weight vector w as follows:
refining estimates using standard interpolation, zero padding,
(cid:26) (cid:12) (cid:12)
wi = 10 iofth(cid:12)eyr(wtiis)e−, µˆ1/2(y)(cid:12)≤3σˆmad(y) (12) oisrsatdilalpstaivtiesfgarcitdorsye.aIrnchfamcte,thaolldcs.onHsoidweerveder,esthtiemiratpoerrsfoarrmeaanbclee
to obtain accuracy of the order of 1Hz or better, when SNR
where µˆ () denotes the sample median operator, and σˆ c
1/2 mad
· and SNR are greater than 20dB.
denotes the normalized Median Absolute Deviation (nMAD), j
Furthermore,westudiedthebehavioroftheproposedmeth-
which is defined, cf. [35], as
odsasthenumberofsamplesN varies.Thesimulationresults
(cid:0) (cid:1)
σˆmad(y)=1.483 µˆ1/2 y 1µˆ1/2(y) . areshowninFig.6g-6i.TheproposedWLSestimatorshows
· −
improved accuracy in small-sample conditions with respect
As the signal model (6) has no trend and a constant offset,
to ULS and PCP, which is advantageous in scenarios with
this choice of weights (12) takes large deviations from the
slowly time-varying parameters. The performance of ULS at
median to be an indicator of a likely outlier and downweights
high SNR and for large N is generally better than PCP and
the corresponding samples.
WLS,asshowninFig.6a,d,andg.Thisisconsistentwiththe
We note that other outlier-rejection strategies are possible
expected behavior. As noise vanishes, in fact, the unwrapping
[37], but are not presented here; see [38] for an extensive
procedure results in a linear trend, which can be estimated
survey. An advantage of the proposed strategy is that it is
statistically efficiently by the linear least-squares method.
automatic, since it does not require any input or parameter
Finally, in order to analyze the effect of outliers, we
tuning by the user. Moreover, in the considered context,
performed a set of Monte Carlo simulations substituting a
another desirable feature of the strategy is its tight integration
fraction of the data with outliers, uniformly distributed in the
in the WLS framework. Finally, the proposed strategy gives (cid:2)3.5 10−6,4.9 10−6(cid:3)s interval. The outliers are inserted at
good results in a practical outlier-prone scenario, as we show · ·
randomly-selected locations, yielding the results in Fig. 6j -
in the experimental evaluation of Section V.
6l. The proposed WLS estimator has considerably improved
robustnesswithrespecttotheotherconsideredmethods.More-
over, it exhibits a threshold behavior, reaching a breakdown
IV. NUMERICALRESULTS
point when the fraction of outliers is approximately equal to
In this section, we explore the properties of the model in 40%. On the other hand, ULS shows a degraded performance
(7) and evaluate the performance of the proposed estimators evenatalowfractionofoutliers,inFig.6gand6h,therefore
in realistic scenarios. validating the proposed WLS approach. PCP operates first in
ULS PCP WLS
5
104 10 102 102
z]
ˆMSEf[Hz]d 110002 ˆfRMSE[Hd 100ˆMSEφ[ns] 110001 MSEˆρ[cm] 101
R −R5 R
10
0 5 10 15 20 25 30
10−2 10−1 SNRc [dB] 100
0 10 20 30 0 10 20 30 0 10 20 30
SNRc [dB] SNRc [dB] SNRc [dB]
(a) (b) (c)
104 102 102
ˆf[Hz]d 102 ˆφ[ns] 101 ˆρ[cm] 101
MSE 100 MSE 100 MSE
R R R
10−2 10−1 100
10 20 30 40 10 20 30 40 10 20 30 40
SNRj [dB] SNRj [dB] SNRj [dB]
(d) (e) (f)
102 101 102
ˆf[Hz]d 100 ˆφ[ns] ˆρ[cm] 101
MSE 10−2 MSE 100 MSE 100
R R R
10−4 10−1
101 102 103 101 102 103 101 102 103
numberofsamplesN numberofsamplesN numberofsamplesN
(g) (h) (i)
105
ˆf[Hz]d 102 ˆφ[ns] 102 ˆρ[cm] 110034
MSE MSE MSE 102
R 100 R 100 R 101
100
10−2 10−1 100 10−2 10−1 100 10−2 10−1 100
fractionofoutliers fractionofoutliers fractionofoutliers
(j) (k) (l)
Fig.6. Simulationresults,averagedover1000MonteCarloiterations.(a)-(c)varyingSNRc,withSNRj =40dB.(d)-(f)varyingSNRj,withSNRc=30dB.
(g)-(i)varyingthenumberofsamplesN,withSNRj =SNRc=40dB.(j)-(l)inthepresenceofoutliers,varyingthefractionofrandomly-insertedoutliers,
withSNRj =SNRc=40dB.
the frequency domain which gives it a degree of insensitivity it was stored and used for offline processing. The apparatus
to outliers. and components used in the experiments are relatively low
cost and easily available. This brings the proposed idea closer
to practice.
C. Figure of Merit
ThenumericalresultsinFig.6showthat,usingtheconsid- As can be seen from Fig. 8, while capturing measurements,
ered RTT clock measurement mechanism and the proposed wecaptureground-truthvaluesoftheparameters.Theground-
WLS estimator, 100 samples are sufficient to achieve an truth values are used for analyzing RMSE of the estimates
estimation accuracy of Hertz-order for clock frequency dif- over measured data. The ground truth of range (ρ) is mea-
ference, nano-second order for initial phase difference and sured using a handheld laser distance meter [41]. Ground-
decimeter-order for ranging in practical noisy and outlier- truth values of master and slave clock frequencies (f ,f )
m s
prone conditions. This can be verified from the curves in are captured using Agilent 53230A frequency counter [42].
Fig. 6(g)-(i), where an SNR =SNR =40dB is considered. Ground-truth of phase (φ) is obtained by capturing the clock
c j
Then, from Fig. 6(a)-(f) it is possible to observe that such stateusingprimitivedelaysinFPGA.Theclockofslavenode
accuracy level is still achieved for SNR as low as 10dB and ispropagatedthroughthesedelaysandthereceivedfirst PING
c
SNR as low as 20dB. More importantly, this accuracy is inasequence PINGsfrommastertriggerscaptureoftheclock
j
achieved by the proposed WLS estimator when up to 30% phase [43].
of the measured data are outliers, as shown in Fig. 6(j)-(l).
We placed the master and slave nodes on two carts in
Notice that N = 100 samples correspond to an observation
line-of-sight conditions in the office-like indoor environment
time of 0.1s for the considered update period of T =10−3s,
s depicted in Fig. 9. For each considered distance a record of
which we have chosen to be of the same order of magnitude
65356 RTT measurement values was acquired. In the experi-
as our practical experimental setup. Such an observation time
mentalconditions,weobservedanSNR varyingfromapprox-
is feasible for most practical applications of joint wireless c
imately 14dB at a distance of 0.5m to approximately 0dB at
synchronization and ranging.
the maximum operating distance of about 4m. Furthtermore,
the SNR is approximately 40dB, and the fraction of outliers
j
D. Sensitivity to Parameter Variations varies from about 5% up to approximately 20%.
To evaluate the performance of the considered estimators
The amplitude of the observed sawtooth waveform is ap-
when the parameter f varies, we provide the simulation
d proximately10nsascanbeseeninFig.5,whichisconsistent
results in Fig. 7, showing that WLS is insensitive to the
with the nominal clock frequency of 100MHz of the FPGA
value of f . On the other hand, ULS and PCP show con-
d boards used. Furthermore, the period is approximately 30
siderable degradation of the frequency difference estimation
samples which, considering the sampling frequency of 5kHz,
performanceasf increases,seeFig.7(a).Asfarastherange
d corresponds to a frequency difference in the order of 200Hz.
is concerned, we note that the model in (7) is linear in the
Suchadifferenceisconsistentwiththetolerancespecification
range ρ, thus it is expected from the theory that the behavior
of the local oscillators used in the experiments.
does not change significantly for different values of ρ. This
is confirmed numerically in Fig. 7(c), where the observed The raw RTT data from the experimental measurement
performance is approximately constant. platformwasdirectlyusedastheinputfortheWLSestimation
algorithm. Whereas, the raw data was pre-processed by a
V. EXPERIMENTALRESULTS coarse outlier-rejection method before being input into the
ULS and PCP algorithms. In particular, the pre-processing
We applied the proposed estimators to experimental data
outlier rejection method consists of the identification of the
obtained from an in-house test-bed which is based on the IR-
outliers by means of the same criterion in (12), and in the
UWBtechnology[18].Inthefollowing,weprovideadescrip-
followingsubstitutionrule:iftheoutlierisisolated,namelythe
tion of the experimental setup, we evaluate the performance
precedingandfollowingsamplesarenotoutliers,thenalinear
of the considered estimators (ULS, PCP, WLS) on measured
interpolation is performed. Specifically, the outlier sample is
data, and eventually we discuss the validation of the proposed
replacedbytheaverageofthefollowingandprecedingsample.
model in (7).
Conversely, if the outlier is not isolated, i.e. if the preceding
or following samples are also outliers, the outlier sample is
A. Experimental Setup
replaced by the median of the entire record.
Experimental setup is shown in Fig. 8. Each of the two
For validation purposes, we also performed a characteriza-
nodes was running on its own clock which was provided by a
tion of ULS and PCP without using the pre-processing step.
local oscillator in the ML505 evaluation board using a Xilinx
Resultsofsuchcharacterization,notpresentedhereforbrevity,
Virtex-5 FPGA [40]. Since the two local oscillators were not
show a considerable performance degradation of one order of
synchronized, a non-zero frequency difference between the
magnitude or more with respect to WLS.
two oscillators was present. In the experimental setup, the
masternodewasprogrammedtorepeattheRTTmeasurement Therefore, in the following subsections, we provide results
procedure with a fixed sampling rate of 5kHz, and then from experimental tests where the pre-processing method is
transfereachresultfromtheTDCtothehostcomputer,where applied to ULS and PCP.
ULS PCP WLS
105 8 4
ˆMSEf[Hz]d 101 ˆfRMSE[Hz]d 10ˆ0MSEφ[ns] 46 MSEˆρ[cm] 3.35
R R R 2.5
100 10−5
0 5 10 15 20 25 30
2
SNR [dB] 2
20 40 60 80 100 20 40 c60 80 100 20 40 60 80 100
fd [Hz] fd [Hz] fd [Hz]
(a) (b) (c)
Fig.7. Numericalsimulations:RMSEperformanceoftheproposedestimatorsfordifferentvaluesoftheparameterfd.Here,SNRc=20dB,SNRj =20dB,
N =200,andnooutliersarepresent.Resultsareobtainedbyaveragingover1000MonteCarloiterations.
B. Frequency Difference Experimental Results PC SlaveClockPhaseCapture(φ)
We acquired RTT data from the system at four master-
slave distances from 1m to 4m at 1m steps. Simultaneously, Data(RTT,y(ti)) CMleoacskurFermeqeunetnscy FrequencyCounter
Capture (Agilent53230A)
the related ground-truth data, i.e. the true master-slave clock
MasterClockFrequency SlaveClockFrequency
frequency difference, was acquired using the setup above Measure(fm) Measure(fs)
and found to be approximately 30Hz with small variations
among datasets at different dista−nces. RoundTripTime(RTT)measurements,y(ti)
Subsequently,theacquireddatawassegmentedintorecords, ρ
eachconsistingof1000consecutivesamples.Eachrecordwas
Master Slave
processed using the three considered methods. The resulting
estimates were compared with the ground-truth, to obtain the Fig.8. Diagramoftheexperimentalsetup.
results shown in Fig. 11. The proposed WLS method outper-
formsboththeotherconsideredmethods,anditsperformance
is relatively insensitive to variations in range. Moreover, the
global RMSE, computed taking into account data acquired
at all distances, is equal to 0.96Hz for the proposed WLS
estimator. Thus, the capability of sub-Hertz relative frequency
measurement in experimental conditions is demonstrated.
Furthermore, a data record illustrating the presence of
outliers is shown in Fig. 10. The presence of outliers can be
explainedbyanalyzingthebehaviorofthereceiveremployed.
Inourexperimentalsetup,infact,theenergy-detectionreceiver
cannot distinguish between a proper response signal and an
external RF interference signal. Therefore, when external in-
terference bursts cross the pre-defined threshold in the energy
detector, they cause spurious detections. Furthermore, the
wideband noise which is inherently present in UWB systems
can also cause spurious detections. Such spurious detections
cause outliers during the RTT measurement procedure. Fig.9. Pictureoftheexperimentalsetup.
C. Phase Experimental Results D. Range Experimental Results
Figure 12 shows error between measured phase and esti- In order to relate the measured RTT to the master-slave
matedphaseforvariousestimators.Sincethephaseparameter range, we performed a calibration procedure based on that
is same for one set of measurements, computing RMSE of described in detail in [18]. The calibration procedure removes
phaseestimatewouldrequiremanysetofmeasurementswhich various systematic delays, particularly delay in RF front end
is practically difficult and time consuming. Hence, instead of of the transceivers. In particular, 1000 RTT samples were
RMSE plots, a few phase error estimates are shown in the acquired at a set of known distances from 0.5m to 4.5m
figure. As can be seen, the phase estimation error is of the at 0.5m steps. Then, the calibration curve was calculated by
order of 1ns for WLS and PCP estimators. As expected from applying a fifth-order polynomial fitting to the experimental
simulationresults,phaseestimationerrorofULSisverylarge. data. The range estimate characterization was then performed
102
ULS
s] PCP
n WLS
[
r101
5000 o
r
r
E
se100
a
4500 h
s] P
n e
T [ ut10−1
RT 4000 sol
b
A
10−2
0 1 2 3 4
3500 Range [m]
Fig.12. Phaseerrormeasurements.
0 50 100 150 200 250
sample index

Fig. 10. Experimental data record showing outliers denoted with magenta
circular markers. The data was obtained withthe master and slave antennas 0.3 ULS
placed5mapart.Thedashedgreenlinedenotesthesamplemean,whichis PCP
clearlyinfluencedbytheoutliers,andthedottedredlinedenotesthesample m]
medianwhichisrobustagainstoutliers. E [ 0.2 WLS
S
M
R 0.1
0
1 1.5 2 2.5 3 3.5 4
distance [m]

(a)
z] 10 ULS 0.6
H
E [ PCP 0.4
MS 5 WLS m]
R or [ 0.2
err 0
0
1 1.5 2 2.5 3 3.5 4 −0.2
distance [m]
(a) 1 2 3 4
distance [m]
(b)
0
z] Fig.13. Experimentalresultsfortheestimationofrange.(a)RMSEofthe
H −1 rangeestimationerrorvsactualrange,forthethreeconsideredmethods.(b)
or [ boxplotoftherangeestimationerrorvsactualrange,fortheproposedWLS
err −2 ethsetimceantotrra.lTrhedebmluarekeddegneosteosftehaechmebdoixand,eannodteththeeto2p5tahndanbdo7tt5otmhpenerdcsenotfiltehse,
blackdashedlinesdenotethemaximumandminimumvalues,respectively.
−3
1 2 3 4
distance [m]
on an independent dataset, acquired at four distances from
(b) 1m to 4m at 1m steps. We denote the latter dataset as the
validation dataset. For each distance in the validation dataset,
Fig.11. Experimentalresultsfortheestimationoffrequencydifference.(a)
RMSEofthefrequencydifferenceestimatesvsrange,forthethreeconsidered weacquiredmultiple1000-samplerecordsandprocessedeach
methods. (b) box plot of the frequency deviation estimation error vs range, recordusingalltheconsideredmethods.Therangeestimation
fortheproposedWLSestimator.Theblueedgesofeachboxdenotethe25th
results are summarized in Fig. 13. As clearly noticeable
and75thpercentiles,thecentralredmarkdenotesthemedian,andthetopand
bottom ends of the black dashed lines denote the maximum and minimum from Fig. 13a, the difference between the three considered
values,respectively. estimationmethodsislessthan1cmand,therefore,negligible.
TheresultingglobalRMSE,whichiscalculatedovertheentire
validation dataset, is 17 cm for the proposed WLS method.