Clock Synchronization Over Networks: Identifiability of the Sawtooth Model

In this paper, we analyze the two-node joint clock synchronization and ranging problem. We focus on the case of nodes that employ time-to-digital converters to determine the range between them precisely. This specific design choice leads to a sawtooth model for the captured signal, which has not been studied before from an estimation theoretic standpoint. In the study of this model, we recover the basic conclusion of a well-known article by Freris, Graham, and Kumar in clock synchronization. More importantly, we discover a surprising identifiability result on the sawtooth signal model: noise improves the theoretical condition of the estimation of the phase and offset parameters. To complete our study, we provide performance references for joint clock synchronization and ranging using the sawtooth signal model by presenting an exhaustive simulation study on basic estimation strategies under different realistic conditions. With our contributions in this paper, we enable further research in the estimation of sawtooth signal models and pave the path towards their industrial use for clock synchronization and ranging.


I. INTRODUCTION
C LOCK synchronization across a deployed network is a pervasive and long-standing challenge [1]- [6].Furthermore, new-generation technologies each require more accurate synchronization.To name a few, i) in cellular communications, synchronization between base stations through a backhaul channel is fundamental to maintain frame alignment and permit handover among neighboring cells, and has been identified as a crucial requirement for distributed beamforming, interference alignment, and user positioning [7], [8], ii) in radio-imaging technology [9], accurate clock synchronization between the sparse chips that form an array is critical, and, in active-sensing 3-dimensional cases [10], [11], it results in low-cost wide-aperture ultra-short ultra-wideband (UWB) pulses, increasing both the angular and depth resolutions of the captured images, iii) in wireless sensor networks [12], [13], synchronization is critical to data-fusion, channel-sharing, coordinated scheduling [14], [15], and distributed control [2] and iv) in distributed database solutions that provide external consistency, clock synchronization accuracy regulates latency, throughput, and performance [5].
Consequently with this wide range of application, theoretical insights on the fundamental limitations of clock synchronization over networks are likely to incite radical innovations in a number of fields.In [16], Freris, Graham, and Kumar established the fundamental limitations of the clock synchronization problem in an idealized scenario.Particularly, (Corresponding author: Pol del Aguila Pla) This work was supported by the SRA ICT TNG project Privacy-preserved Internet Traffic Analytics (PITA).
Pol del Aguila Pla is with the Center for Biomedical Imaging, in Switzerland, and with the EPFL's Biomedical Imaging Group in Lausanne, Switzerland (email: pol.delaguilapla@epfl.ch).Pol performed the work while at the KTH Royal Institute of Technology.
Lissy Pellaco, Peter Händel and Joakim Jaldén are with the Division of Information Science and Engineering, School of Electrical Engineering and Computer Science, KTH Royal Institute of Technology, Stockholm, Sweden (e-mail: pellaco@kth.se,ph@kth.se,and jalden@kth.se).
given a network of nodes with noise-less affine clocks and fixed unknown link delays that exchange time-stamped messages, [16] i) showed that clock synchronization was only possible if the link delays were known to be symmetric, and ii) characterized the uncertainty regions of the clock synchronization parameters under different hypotheses.In this paper, we analyze the same problem from a perspective that is closer to real implementation.In short, we analyze the twonode joint clock synchronization and ranging problem [17]- [19] with noisy round-trip time (RTT) measurements without time-stamps [20], for a node design originally proposed in [21] to improve ranging accuracy.The resulting analysis has several advantages.First, because protocols without time stamps require only minimal transmissions of very short pulses carrying no information, the resulting technology minimizes communication overhead, and is beneficial in applications in which the data-rates are critically needed for other uses [6, p. 29].Second, because we consider hardware specifically tailored to ranging accuracy, we reveal how applications that require this accuracy, such as cooperative localization [22], positioning [23], and control [2], can harness the same hardware and protocols for synchronization.Third, the analysis is more realistic, because it takes into account the real-world stochasticity of the measurements.In particular, in our analysis of the problem we i) unveil the need for symmetric delays in RTT-based protocols, in a direct parallel to the discovery in [16], ii) find novel results on the identifiability of sawtooth signal models under diverse conditions, which are of interest by their own right to chaotic system analysis [24], [25] and control, and iii) provide performance references to guide practitioners in their use of this technology.
In summary, in this paper we first derive from basic principles a model for RTT measurements between two nodes equipped with time-to-digital converters (TDC) in a network with fixed, unknown link delays (Theorem 1).Then, we shift our focus towards an encompassing family of signal models, i.e., sawtooth signal models, when one considers different stochastic effects.In this context, we provide results on the identifiability of these models, both negative (Lemma 1) and positive (Theorem 2), under different noise conditions.Here, we obtain the surprising result that the presence of a noise term inside a non-linear model term makes said model identifiable.We then shift the focus again towards clock synchronization and ranging, and we provide a thorough and verifiable empirical evaluation of the basic estimation techniques we propose in [26] to exploit the sawtooth model (Figs.8-11, implementation accessible in [27]).These empirical results, together with the approximated Cramér-Rao lower bounds derived in [26], are clear and simple performance references for clock synchronization and ranging using sawtooth models.Such performance references are of use to both engineers that use this technology and to researchers aiming to develop estimation techniques for sawtooth signal models.

A. Notation
Discrete random processes will be in uppercase letters and square brackets, such as Y [n], while deterministic sequences, e.g., realizations of said processes, will be lowercase with square brackets, i.e., y[n].For both these sequences, the notation will be simplified by omitting the discrete time index when it can be established by context.Vector random variables will be bold uppercase letters, e.g., Y, while deterministic vectors will be bold lowercase letters, e.g., y ∈ R N .Functions, on the other hand, will be non-italics lower-case letters, e.g., the probability density function (PDF) of a vector random variable Y on a parametric family with vector parameter θ will be f Y (y; θ).Through the paper, we view the modulus equivalence as a function, i.e., we note the mapping x → y ∈ [0, a) | (y = x) mod(a) as mod a : R → [0, a).

II. THE SAWTOOTH MODEL
In applications in which high ranging accuracy with low communication overhead is desired, a low-cost solution in terms of both complexity and power consumption is using node designs that include TDCs to measure RTTs [21].Indeed, such sensors were successfully incorporated in a prototype system to aid firefighters by providing on-site infrastructure-free indoor positioning [23].TDCs, however, induce an asymmetry between the rate at which nodes can measure time and the rate at which they can act upon their environment.This asymmetry generates an unexpected waveform in the sequence of RTT measurements over time.This phenomenon was first reported by [19], where a sawtooth model was proposed and empirically validated, and possible applications to clock synchronization over networks were identified.In this section, we reintroduce the design of [21] and derive the sawtooth model from a few simple assumptions.
Consider now the design of [21], described in Fig. 1.Here, each node or sensor N has a processing unit, a transceiver and a TDC.With this design, a sensor can measure RTTs at the resolution of the TDC, usually in the order of ps, much finer than the period of the processing unit's clock, usually in the order of tens of ns.Besides the clear advantage of this design for ranging through RTT measurements, this creates an interesting asymmetric behavior of the node as an agent and as  Initially proposed in [21], this design uses an independent time-to-digital converter (TDC) to accurately measure round-trip times (RTT).The resulting node N can measure RTTs with a much finer time-resolution than it can react to incoming pulses.Indeed, in order to react to an incoming pulse, N must first access the TDC's memory, process the reading, and decide to send a pulse, all of which require waiting until its next clock cycle.a measuring device.As we will show below (Theorem 1), this asymmetry produces a sawtooth waveform in the measured RTTs that depends on the synchronization parameters, and, under reasonable assumptions, leads to a viable system for clock synchronization over networks.As a consequence, such a design may also be considered for wired networks, where ranging information is usually not relevant.A final by-product of the inclusion of a TDC in the design in Fig. 1 is that we can consider that each node has perfect knowledge of its own clock period, which it measures directly with its TDC.

A. Deterministic model
Consider two nodes, M and S, designed as N in Fig. 1, in a network (wired or wireless).These two nodes execute the RTT measurement scheme illustrated in Fig. 2. In this scheme, M measures the RTT between itself and S by sending pulses (a.k.a.pings) to S and using its TDC to accurately record when a response (a.k.a.pong) is received from S. In particular, M sends a pulse at some of the times at which its clock has upflanks, i.e., at the times t n = KT M n.Here T M > 0 [s] is M's clock period, the sampling factor K ∈ N is designed to determine the sampling period T s = KT M [s], n ∈ N is a discrete-time index, and we assume without loss of generality that M's clock phase offset is zero, i.e., φ M = 0 [rad].If we assume that the delays involved in the pulse traveling from M to S accumulate to a constant value δ → [s], the n-th pulse arrives at S and is recorded in its TDC at time t n + δ → .Nonetheless, S will not be able to access the TDC's memory Fig. 2.
Example of the round-trip time (RTT) measurement scheme in Section II-A.M sends ping pulses to S at its clock upflanks every Ts = 2T M .The ping is recorded in S's time-to-digital converter (TDC).At its next clock upflank, S accesses its TDC and starts a delay of δ 0 = T S before responding with a pong pulse.The pong is recorded in M's TDC as soon as it arrives.before its next clock upflank, and consequently, any action by S will be further delayed until t n + δ → + ∆ n .Here, ∆ n ≥ 0 [s] is the time remaining until S's next clock upflank.If we consider that S has a clock with period T S > 0 [s] and phase φ S [rad], i.e., an offset delay of ϕ S = T S φ S /(2π) [s], then S has its clock upflanks at those times that are at an integer number of periods away from ϕ S , i.e., at the times τ ≥ 0 when mod T S (τ + ϕ S ) = 0. Furthermore, we know that ∆ n ≤ T S , as T S is the time between consecutive upflanks.Consequently, to obtain a closed-form expression for ∆ n , we need to find ) for any a ≥ 0 and b, c ∈ R, the condition for t n + δ → + τ to be the time of one of S's clock upflanks can be rewritten as τ = QT S − mod T S (t n + δ → + ϕ S ) for some Q ∈ Z.Then, because mod T S (t n + δ → + ϕ S ) < T S and τ ∈ (0, T S ), we conclude that We now allow for a known delay δ 0 [s] to be introduced by S, which can account for any processing required to read the TDC's state and prepare the new pulse, and will usually be an integer number of S's clock periods, i.e., δ 0 = K 0 T S .Finally, as we did for the ping pulse, we consider δ ← to express the fixed delay for a pong pulse from S to reach M and be captured by the TDC.In conclusion, if we disregard the effect of the resolution of the TDC, which is usually four orders of magnitude finer than that of the nodes' clocks, the n-th RTT measurement will amount to where δ ↔ = δ → + δ ← .We summarize our result in the following theorem.
Theorem 1 (Deterministic RTT measurement model): Consider two nodes M and S designed as specified in Fig. 1.Then, if M and S follow the RTT measurement protocol specified above, f d = 1/T S − 1/T M and δ ↔ = δ → + δ ← , the n-th RTT measurement y det [n] can be expressed as Proof: From (1) we have that Here, we have used that T s = KT M and ϕ S = T S φ S /(2π) in (5), that mod 1 is periodic with period one in (6), and that f d = 1/T S − 1/T M in (7).Finally, (3) follows from substituting (7) in (2).
In conclusion, by running the RTT protocol specified in Fig. 2, M obtains data intimately related with the parameters it needs to predict S's clock signal, i.e., to synchronize to S. Indeed, having measured its own clock period T M using its TDC, f d reveals T S , which together with φ S characterizes S's clock signal completely.The goal of our study is to establish under which conditions M will be able to simultaneously estimate these parameters.
Other formulations of the model (3) in terms of the usual synchronization parameters for affine clocks, i.e. the clock skew α S = T S /T M and the offset delay ϕ S , including the general expression for when ϕ M = 0, can be found in the supplementary material to this paper.Nonetheless, the expression in (3) remains the most practical, because it expresses the compromise between the sampling period T s and the frequency difference f d of the system, which will prove to be relevant to our analysis.
Incidentally, under the simple assumption that T s ≥ δ 0 +T S , which can be guaranteed under any reasonable f d if K > K 0 + 1, if we assume that S's TDC starts measuring every time S sends a pong and stops measuring when the next ping is received, the n-th measurement x[n] taken by S's TDC can be expressed as As we will see, this will imply that even while M is leading the RTT measurement protocol, S could still perform frequency synchronization.Nonetheless, we will not consider S's TDC measurements for most of the paper, and we will instead focus on determining the conditions under which M can achieve full synchronization and ranging.More details on the derivations of ( 3) and ( 8) can be found in the supplementary material to this paper.Unrealistic parameters were used to obtain a cheap-to-compute, simple, representative figure.For more details on how this simulation was performed see this project's repository at [27].
In this project's repository, accessible at [27], we validate (3) by simulating an ideal physical system as described above and verifying the exact correspondence between the model and the obtained measurements.In Fig. 3, we show the fits of ( 3) and ( 8) on the TDC measurements of M and S throughout a simulated run with noisy clock periods and noisy transmission delays.

B. Stochastic model
In Fig. 4, we show real RTT data obtained in [19] from the ultra-wide band testbed of [21] using this RTT scheme, accompanied by an example model fit.Given the observed signal and its expected shape, a simple observation is that S's clock was faster than that of M in the specific experimental set-up, because the ramps in the sawtooth signal have negative slope, which implies that f d > 0. Fig. 4 also exemplifies the two distinct effects that random deviations of the physical parameters can produce on the data.On one hand, large jumps of approximately T S in the measured RTT are observed (effect i]) if a random deviation influences the specific clock period at which S reads the arrival of a ping pulse from its TDC, i.e., it changes which is the first up-flank in S's clock after the ping pulse arrives.On the other hand, if this does not happen, random deviations appear directly in the signal as additive noise (effect ii]).From a modeling perspective, these two effects are not easily represented distinctively.Indeed, variations of the transmission time from M to S, δ → , or jitter in any of the two clock periods, T M or T S , could lead to any of the two described effects, while variations of the transmission time from S to M, δ ← , can only ever lead to effect ii].In this paper, we will consider the effect of random variations on the physical parameters, as well as the quantization by the TDC, in the form of two additive white noise processes W [n] and V [n], inside and outside the nonlinearity, respectively.In short, our stochastic model for the RTT measurements taken by M is For simplicity, we will assume that W [n] and V [n] are zeromean Gaussian processes with respective standard deviations σ w and σ v and we will consider them independent.Analyzing the effect of the existing dependence between them, or evaluating the magnitude or effect of this dependence, is outside of the scope of this paper.Only adding V [n], i.e., setting σ w = 0, could explain the two effects explained above for most sample indices n.Nonetheless, as shown by the indicators of the maximum and minimum of the model fit in Fig. 4, the experimental RTT measurements are not bounded in the range indicating that the noise term outside the non-linearity W [n] is necessary.Furthermore, random variations in δ ← or δ 0 cannot be meaningfully represented by V [n], since variations of these parameters of any magnitude will never affect which upflank of S detects the ping pulse.Fig. 5 exemplifies the effect of each of the noise terms by showing two realizations of our stochastic model (9), one in which σ w = 0 and σ v > 0, and one in which σ w > 0 and σ v = 0.
In order to simplify the notation for the rest of the paper and abstract some of our theoretical results, we will express the stochastic sawtooth model in terms of four generic parameters, an offset α ∈ R, a non-zero amplitude ψ ∈ R\{0} with known sign sign(ψ), a normalized frequency β ∈ [−1/2, 1/2), and a normalized phase offset γ ∈ [0, 1).In other words, we will express the sawtooth signal model as (10) with W [n] and V [n] independent additive white Gaussian noise processes.An empirical analysis verifying this model (10) on real data from the testbed of [21] can be found in [19].
Here, the restriction of the β and γ parameters simply reflects the maximum ranges that we can expect to distinguish, given the periodicity of mod 1 (•) as a function.Indeed, adding any integer factor of n inside the modulus, or any integer by itself, will not change Y [n], and establishes an equivalence of period  9), exemplifying the effects of additive noise inside and outside mod 1 (•).On one hand, σw = 0 and σv > 0 leads to many high jumps around the wrapping points.On the other hand, σw > 0 and σv = 0 leads to a signal that is not bounded by the minimum and maximum values of the deterministic model (3) (shown in dash-dotted horizontal lines).For further examples of the effects of noise in measurements following our stochastic model ( 9), as well as the effects of randomness in the physical quantities described above, see this project's GitHub repository at [27].
Second, if we consider this restriction and examine the relation between the parameters of both models, we observe that and, incorporating that ) and that T M , T s and δ 0 are known, Clearly, then, unless further constraints relating δ → , δ ← and f d are given, it is impossible to recover δ → , δ ← , f d , and φ S from α, ψ, β and γ.In the context of clock synchronization over networks, this is equivalent to the impossibility result of [16], which studied the uncertainty sets where the synchronization parameters are known to lie given time-stamped message exchanges under different conditions.An analysis similar to that in [16] under idealized, noise-free conditions could be reproduced for (3), but is outside of the scope of this paper.
In contrast, we will provide an analysis of identifiability when every physical parameter can be subject to noise.In fact, this analysis will reveal that synchronization with the sawtooth signal model requires a certain level of randomness, i.e., it is impossible without it.Consequently with the discussion above, then, we will assume that δ → is given when one knows δ ↔ and f d , as it happens in a number of applications.For example, in wireless sensor networks, one may generally consider that all nodes are equal and the channels between any two of them are symmetric, and thereby one can assume is the unknown range between M and S in the communication medium and c [m/s] is the speed of light in the medium.Even in this context, true line-of-sight communication is not a requirement, and one only needs to assume that the multipath is not dense and the direct path is not fully blocked, as in such a case the TDC will trigger on the first pulse, corresponding to the shortest path.These assumptions are consistent with ultra-wide band pulses such as the ones used in [19].When convenient in the paper, we will use this assumption combined with δ 1 ≈ 0, and consider the ranging problem of [21], [23] jointly with clock synchronization [19].In the following section, we will characterize the model (10) statistically, providing conditions for its identifiability.Our aims in doing that are 1) to present novel results on the sawtooth signal model, and 2) to provide guarantees for the design of practical synchronization systems using nodes modeled by the design in Fig. 1.

III. IDENTIFIABILITY OF THE SAWTOOTH MODEL
Identifiability is a basic requirement on any statistical model that relates to the minimal conditions that make parameter estimation a reasonable goal [28].
Def.  (2) , ∀y ⇔ θ (1) = θ (2) .( 14) That is, the mapping between the parameter θ and the distribution specified by Y θ is one-to-one.If (14) is not met, the data observed when the parameter value is θ (1) and the data observed when the parameter value is θ (2)  have the same distribution, and therefore, distinction between these two parameters from observed data is impossible.Unintuitively, even if (14) is not given, one could possibly design good estimators for θ.Specifically, as long as the selected metric in the space of parameters Ω does not assign much importance to the difference between the pairs θ (1) and θ (2)  that do not fulfill (14), estimation could remain a sensible objective.In this paper, the data model Y θ is defined by (10), and the considered parameters are θ = [α, ψ, β, γ] T .Hence, this section will be dedicated to establishing under which conditions, in terms of the values of σ w and σ v in (10), a one-to-one relation between θ and the distribution of the data T can be ensured.

A. Unidentifiability without inner noise
In order to analyze the relation between θ and f Y (y; θ), we will first consider the simplifying assumption σ v = 0.This is an unrealistic assumption under most applicable uses of the sawtooth model (10), including that of clock synchronization, but it will be useful for our analysis.We will show that under this assumption, (10) yields an unidentifiable model Y θ in which the effect of α and γ cannot be fully distinguished in the observed data Y ∼ Y θ .
Consider first (10) with σ v = 0 and observe that then, Y θ = N µ θ , σ 2 w I N , where I N is the N × N identity matrix and with the modulus operation mod 1 (•) applied component-wise, T ∈ R N , and n = [0, 1, . . ., N − 1] T .Because the normal distribution is fully characterized by its location and scale parameters, we know that changes in θ will only affect the distribution in terms of its location, controlled by its mean µ θ .Consequently, the condition for identifiability in Def. 1, i.e., (14), can be restated as µ θ (1) = µ θ (2) ⇔ θ (1) = θ (2) .Lem. 1 establishes that, when σ v = 0, there are changes in α and γ that violate this condition.
Lem. 1 (Unidentifiability of (10) when given by (10) 1), when the inner noise is disregarded, i.e., σ v = 0.Then, Y θ is unidentifiable.In particular, there are different combinations of α and γ that yield the same distribution of Y under Y θ .
Fig. 6 illustrates the idea of our proof of Lem. 1 when ψ > 0, and shows, in Fig. 6(b), changes in α and γ that cannot be distinguished in the mean of Y for a simple example with N = 3.Consequently, Fig. 6(b) serves as a straightforward counterexample to the identifiability of Y θ when σ v = 0. Note that both the result in Lem 1 and the counter-example of Fig. 6(b) are valid when σ w = 0, i.e., when the model is deterministic as in (3).In contrast with our result in Theorem 2, this implies arbitrarily accurate synchronization could be impossible in an idealized scenario without noise.
In our proof of Lem. 1, we exploit the formal definition of mod 1 (•) to claim that its value will always be strictly less than one, and therefore, we obtain the margin + under which changes in α of the same sign as ψ and positive changes in γ are not distinguishable.However, the real limitation on identifiability is given by the points closest to the discontinuity from both sides, and, in most cases (i.e., γ = 0 and for most βs), a similar margin − can be obtained under which changes  I).To access the fully reproducible code to generate this figure, see this project's repository [27].For comparison, the figure shows the Cramér-Rao lower bounds (CRLB) for estimating an offset in white noise, which here serves as a reference for the estimation of α and γ (see the supplementary material of this paper for details).Here, we observe that while the MSE in estimating initially decays as predicted by the CRLB, it stabilizes around the variance of a uniform noise of width * + .Furthermore, the MSE in estimating γ becomes worse and stabilizes around the same value when the estimation of α reaches that level. in α of the sign opposite to ψ and negative changes in γ are not distinguishable.
While our analysis is concerned with a fixed value of , the lack of identifiability stated in Lem. 1 may be less problematic in an asymptotic regime.In particular, if increasing the sample size tends to reduce the segment at the left of the non-linearity without any sample, i.e., + → 0 when N → +∞, the size of the changes in α and γ that cannot be distinguished in the data would decrease with N , making the model identifiable in an asymptotic regime, or at least taking away importance from our proof of non-identifiability for large N s.In particular, if we consider the sequence of elements of the vector (15), µ θ n , we obtain what is known in dynamical systems as the orbit of a rotation of the circle.Then, if β ∈ R \ Q we have that the orbit is minimal [30, ch.1.3., proposition 1.3.3.],i.e., that the set {µ θ n } N 1 when N → +∞ is dense in [0, 1), and thus, + → 0. This contradicts the intuitive notion of finite-sample identifiability as a necessary condition for consistent estimators to exist, seen, for example, in [31, p. 62].In contrast, when β ∈ Q, (15) is periodic and hence In particular, if β = ±M/Q with M and Q two co-prime naturals, then (15) is periodic with minimal period Q, and increasing the sample size beyond Q will not result in any further reduction of + , i.e., any further improvement from an identifiability perspective.In the case of clock synchronization, this specific case corresponds to coherent sampling, in which T s f d = ±M/Q.The effect of this specific case in the estimation error of a global grid search (GGS) strategy on the prediction mean square error of the model (3) (see [26]) when ψ and β are known is illustrated in Fig. 7.
Our analysis has assumed that α was part of the parameter vector θ, and that one wanted to recover it.Although this can be the prominent case in many applications of the sawtooth model, e.g., synchronization in wireless sensor networks or networks of autonomous vehicles, other applications may consider α to be known.Within synchronization, this would be the case of base-station synchronization in cellular networks, in which the backhaul links will most likely have a known and stable transmission delay.This would invalidate the identifiability analysis in Lem. 1, and under some additional conditions, Y θ could be shown to be identifiable.Regardless, in the next section we analyze the full model in the presence of noise inside the mod 1 (•) non-linearity, i.e., with θ = [α, ψ, β, γ] T when σ v > 0, and show its identifiability.
Because the proof of Theorem 2 is rather technical, we place it in the Appendix A. However, it is worthwhile to note here that it is not limited to the case in which W [n] and V [n] are Gaussian processes.Indeed, the statements in there apply mutatis mutandis under a wide variety of distributions for W [n] and V [n].In particular, any W [n] consisting of independent and identically distributed (IID) samples from any location-scale family with some reference PDF ϕ(w) with unbounded support will allow for the conclusion in (19).Furthermore, such a W [n] together with any V [n] consisting of IID samples from a location family that accepts a PDF and leads to a monomodal distribution after wrapping with mode equal to the location parameter, e.g., IID Cauchy distributed samples [32, p. 51], will also preserve all the statements therein.Nonetheless, to our knowledge, the literature mostly considers timing errors to be Gaussian (see, among others, [5], [7], [8], [17], [18], [33], [34]), with little empirical incentive to consider other models.
The contrast between Lem. 1 and Theorem 2 is initially non-intuitive.Indeed, it implies that the presence of noise inside the non-linearity improves the theoretical condition of the estimation problem.This result recalls the popular theories of stochastic resonance for testing and estimation [35]- [39] and of dithering for improving the signal-independence of quantization noise [40], but is, in fact, less expected.In summary, both these theories delve into using noise to improve the performance of knowingly suboptimal strategies.In contrast, our identifiability result reveals how the inclusion of noise makes the data more informative with respect to the underlying parameters.To understand this result, one must first consider that it is very specific to sawtooth models, as it relies on that for any x ∈ R \ Z, there is an > 0 such that mod 1 (x + ) = + mod 1 (x).In other words, small changes in phase and offset are equivalent around almost every x.In order to have any hope to distinguish them, one needs to guarantee that the points in which their effect differs play a role in shaping the distribution of the data.A phase noise V [n] with long enough tails provides this guarantee, because changing the phase will alter the wrapping of the tails of the distribution through these non-linear points, while changing the offset will not.

IV. NUMERICAL RESULTS
In [26], we propose two estimators for the parameters of the sawtooth model (10).In particular, we first introduce an heuristic technique that is computationally light, intuitive, and surprisingly robust, which we call periodogram and correlation peaks (PCP).In this approach, one uses peaks in the power of the discrete Fourier transform to estimate the normalized frequency, and peaks in a circular correlation of the first estimated period to estimate the phase parameter.Using PCP as a starting point, we also introduce the local grid search (LGS), a computationally heavy method that improves the final performance by exploring a grid of normalized frequencies and phase parameters around the PCP estimate, and picks the combination that minimizes the model's prediction mean square error, i.e., how well the model fits the observed signal.In both techniques, the offset parameter is chosen as the least squares estimate for each pair of frequency and phase estimates.For more detail on these techniques, and for the values of their parameters in the experimental results presented here, see [26], [27].
Here, we evaluate these estimators in a series of realistic simulation studies of clock synchronization and ranging.In the exposition of these results we intend to aid 1) engineers that aim to apply this technology, by providing reasonable expectations on its current possibilities, and 2) theoreticians that aim to develop estimators for the parameters of the sawtooth signal model (10), by revealing the strengths and pitfalls of the techniques that are currently available.With respect to 1), we include in all our results references that make it easier to identify different standard performance measures, regardless of the scale and aspect of each figure.In particular, a) in figures reporting ranging performance, i.e., MSE(ρ), horizontal lines corresponding to standard errors of 1 cm or 0.1 cm are shown, b) in figures reporting frequency-difference estimation performance, i.e., MSE( fd ), horizontal lines corresponding to standard errors of 10 ppb (1 Hz) and 1 ppb (0.1 Hz) of M's frequency, 1/T M = 100 MHz (see Table I), are shown, and c) in figures reporting S's phase estimation performance, i.e., MSE( φS ), horizontal lines corresponding to standard errors in ϕ S that are one or two orders of magnitude below T S are shown.In reference to c), we intentionally report performance on the estimation of the phase parameter φ S instead of the absolute time delay ϕ S = T S φ S /(2π).In our opinion, this is a better and fairer measure of how useful a specific clock synchronization technique is, because the scale of the errors in ϕ S will always be dominated by the order of magnitude of T S .
In other words, if T S ≈ 10 ns, even guessing φ S at random between 0 and 2π achieves errors in ϕ S that are on the order of ns.In order to streamline the exposition of this section and avoid unnecessary repetitions of the experimental conditions, we detail the default values for the physical and simulation parameters in Table I.Note here that we chose to randomize the most critical parameters, so that the dependences shown in Figs.8-11 can be considered the average behavior within the selected ranges.

A. Sensitivity to the parameter values
One of the weaknesses of the estimation approaches we present in [26] is low performance when f d is small.In particular, because the PCP includes a mean-removal step before the DFT, the low frequencies are supressed.If f d is small, then, the peak we aim to detect in the DFT is most likely dampened and we detect noise instead.This is visualized in Fig. 8, which shows the average performance in the estimation of f d by both PCP and LGS as a function of f d , when all other parameters are randomized according to Table I.Here, we see that the estimators fail when f d ∈ [−10, 10).Furthermore, we also observe the effect of the grid underlying the DFT on the PCP frequency estimate.In particular, for a given N , the PCP will only propose as estimates frequencies that are at one of the points in the DFT (e.g., −200 Hz and 200 Hz in Fig. 8), and so frequencies that are close to those will be better estimated than those that are far.This same phenomenon is observed in the experimental results we present in [26] for varying N and a fixed combination of physical parameters.Remarkably, the performance obtained with LGS remains below 1 ppb of 1/T M for most frequencies f d outside the [−10, 10) area.
A weakness of the measurement protocol described in Section II-A is that, due to the lack of time-stamps in the exchanged packets, the time origin shift between M and S's clock appears only in terms of a phase term inside the mod 1 (•) function, e.g., φ S if we assume that ϕ M = 0.If one desires absolute time synchronization, this can have dire consequences.While φ S = 2π(1 − ξ/2) and φ S = πξ for some small ξ > 0 are only 2πξ away under the non-linear wrapping, their difference implies errors in ϕ S of T S (1−ξ).In order to take this effect into account, we use the conventional Euclidean distance to quantify the error for the phase term φ S , without taking into account the wrapping effect.The ensuing increase of error around the extremes can be seen plainly in Fig. 9.This weakness is characteristic of protocols for clock synchronization without time-stamping.In turn, however, these protocols minimize the communication overhead and are more robust to malicious nodes [6, p. 29].

B. Average consistency
An important property of estimators is consistency, i.e., the convergence of the MSE towards zero as the sample size increases.In order to characterize consistency beyond the example given in [26] with fixed physical parameters, in Fig. 10 we report the average performance for randomized parameter values for PCP and LGS in the estimation of the clock synchronization and ranging parameters as a function of the sample size.However, we sample frequency values in a reduced range, i.e., f d ∈ [−200, −10) ∪ [10,200), in order to avoid the instability of our estimators when f d ∈ [−10, 10).This provides an impression of the estimators' performance within their range of probable use, and significantly reduces the amount of Monte Carlo repetitions necessary to obtain intelligible results.Fig. 10 strengthens the conclusions from the analysis of the empirical results in [26].Indeed, 1) the results suggest that both PCP and LGS are consistent for the estimation of ρ, f d and φ S , 2) LGS seems to have an asymptotically efficient rate of convergence with N for the estimation of both ρ and f d , while the rate of convergence stagnates for φ S .Here, we use as a reference the assymptotic rates of convergence for a linearized version of ( 9), which we derive from the Fisher information matrix in [26].
Finally, from a practical point of view, average estimation errors under 0.1 cm in the range parameter can be expected For reference and comparison, the convergence rates given by the inverse Fisher information matrix we present in [26] are portrayed by lines with the corresponding slope adjusted to fit the data.
for N ≥ 1000, and average estimation errors under 1 ppb in the frequency parameter can be expected for N ≥ 1500.For the phase parameter φ S , estimation errors below 2π/10 can only be expected with LGS and for N ≥ 1000.

C. Sensitivity to the inner and outer noises
All the results in 8 and 9 (as well as those in [26]) were obtained under the optimistic noise conditions SNR in = 40 dB and SNR out = 20 dB.In Fig. 11   Fig. 11 reveals that range estimation performance decays progressively with the decrease of SNR out throughout the investigated range for both PCP and LGS.In contrast, the decrease of SNR in creates a progressive decay of performance only up to a breaking point around SNR in = 10 dB.In this regime (very low SNR in ), one could consider modeling the non-linear term in (10) as a uniform noise term and employ techniques tailored to the estimation of offsets in Gaussian plus uniform noise [41].Practically, both PCP and LGS achieve estimation accuracies below 0.1 cm only when SNR out ≥ 10 dB for SNR in = 40 dB, and only when SNR in ≥ 15 dB for SNR out = 30 dB.
For frequency estimation, we see that both PCP and LGS are very robust to a decrease of SNR out up to a breaking point around SNR out = 10 dB, and LGS maintains an improvement of 15 dB over PCP for any SNR out above the breaking point.
In contrast, when SNR in degrades, the breaking point for both PCP and LGS is SNR in = 20 dB, and the improvement of LGS over PCP increases gradually as SNR in increases, only reaching 15 dB when SNR in = 40 dB.Practically, on one hand, PCP achieves estimation accuracies below 10 ppb only when SNR out ≥ 10 dB for SNR in = 40 dB, and only when SNR in ≥ 20 dB for SNR out = 30 dB.On the other hand, LGS consistently improves on it, reaching estimation accuracies below 1 ppb only when SNR out ≥ 10 dB for SNR in = 40 dB, and only when SNR in ≥ 30 dB for SNR out = 30 dB.Promising directions for improving frequency estimation could come from two different fronts.First, one could generalize the framework in [41] to admit frequency estimation, aiming to protect the resulting method for noises inside the nonlinearity beyond the breaking points in SNR in .Second, one could use the structure of the sawtooth signal in techniques similar to [42] to extract more information from the spectrum of the data.
For phase estimation, our results are not conclusive, because the randomization of both the frequency and phase parameters result in wide variability that would require further Monte Carlo averaging.Furthermore, while the unfavorable region of frequency parameters f d ∈ [−10, 10) has been avoided, the phase parameters are still sampled from the whole range φ S ∈ [0, 2π), which includes the very challenging regions around the wrapping point (see Fig. 9).However, the results seem to suggest a progressive decay of performance for both PCP and LGS when either SNR out or SNR in degrade.Furthermore, we see that phase estimation is much more sensitive to a degradation of SNR in than that of SNR out .Possible improvements of phase estimation could be expected using LGStype techniques complimented with better frequency parameter estimates.

V. CONCLUSIONS
In this paper, we provide practical and theoretical insights on the fundamental limitations of clock synchronization over networks in applications that require high ranging accuracy and low communication overhead.From a practical standpoint, we show from first principles that using TDCs for measuring RTTs enables the use of a mathematical model that leads to very accurate ranging (e.g., accuracies beyond 0.1 cm for N ≥ 1000 samples in Fig. 10) and remarkable frequencysynchronization performance (e.g., accuracies of 1 ppb for N ≥ 1500 samples in the same figure), all with very simple estimation techniques.Furthermore, we point at promising directions of research that could hold the key to the improvement of these performance values, such as the extension of [41] to improve range estimation when SNR in is very low (e.g., under 10 dB in Fig. 11), or the use of techniques similar to [42] to improve frequency synchronization.In fact, in [26] we provide a reference on the potential for improvement in the form of approximated Cramér-Rao lower bounds based on a linearization of the sawtooth model and standard results for Gaussian models.Additionally, we pinpoint the weaknesses of both the model and our estimation strategies.First, we acknowledge that more research is needed to consistently obtain accurate phase synchronization (phase estimation accuracies beyond 2π/100).Second, we identify that schemes that rely on this mathematical model are not best suited for absolute time synchronization because the offset delay only appears as part of a phase term, which leads to wrapping errors.
From a theoretical standpoint, we establish that RTT-based schemes for clock synchronization are characterized by similar fundamental limitations as those relying on time-stamped message exchanges, previously discovered by [2].Namely, synchronization is impossible with unknown path delays δ → and δ ← if one cannot assume some relation between δ → , δ ← and the frequency difference f d , e.g., that they are symmetric, i.e., δ → = δ ← .Furthermore, we have discovered a surprising property of sawtooth signal models, i.e., that adding noise inside the non-linearity allows for the joint identifiability of its offset and phase terms (cf.Lemma 1 and Theorem 2).This result challenges the convention of random variations as a negative component of a model, and is of use beyond our application domain.

APPENDIX A PROOF OF IDENTIFIABILITY
Proof of Theorem 2, Identifiability of (10) when σ v > 0: The backward implication of identifiability, i.e., that the same parameters lead to the same data distribution (see Def. 1), is immediately clear from (10).

Fig. 1 .
Fig. 1.Internal design of any node N considered throughout this paper.Initially proposed in[21], this design uses an independent time-to-digital converter (TDC) to accurately measure round-trip times (RTT).The resulting node N can measure RTTs with a much finer time-resolution than it can react to incoming pulses.Indeed, in order to react to an incoming pulse, N must first access the TDC's memory, process the reading, and decide to send a pulse, all of which require waiting until its next clock cycle.

Fig. 3 .
Fig.3.TDC measurement buffer, in units of the TDC's clock, taken by M and S in a simulation of the protocol described by Fig.2and Section II-A.Unrealistic parameters were used to obtain a cheap-to-compute, simple, representative figure.For more details on how this simulation was performed see this project's repository at[27].

Fig. 4 .Fig. 5 .
Fig.4.RTT measurements from the ultra-wide band testbed from[21], compared to a fit of the deterministic model(3).In dash-dotted horizontal lines, the maximum and minimum values of the model.

12 Fig. 7 .
Fig.7.MSE( αGGS ) and MSE(γ GGS ) against the sample size N when σv = 0, ψ = 1 and β = M/Q with M = 1 and Q = 10.αGGS and γGGS are the result of a global grid search (GGS) on the prediction mean squared error with model (3) (see[26]) with 1000 discretization points for γ ∈ [0, 1) when β and ψ are known.Results obtained from 2000 Monte Carlo repetitions and SNRout = 5 dB (see TableI).To access the fully reproducible code to generate this figure, see this project's repository[27].For comparison, the figure shows the Cramér-Rao lower bounds (CRLB) for estimating an offset in white noise, which here serves as a reference for the estimation of α and γ (see the supplementary material of this paper for details).Here, we observe that while the MSE in estimating initially decays as predicted by the CRLB, it stabilizes around the variance of a uniform noise of width * + .Furthermore, the MSE in estimating γ becomes worse and stabilizes around the same value when the estimation of α reaches that level.

Fig. 8 .
Fig. 8. Result of 300 Monte Carlo repetitions in the conditions specified in Table I evaluating the MSE in the estimation of f d by both PCP and LGS with respect to its actual value.

2 Fig. 9 .
Fig. 9. Result of 300 Monte Carlo repetitions in the conditions specified in TableIevaluating the MSE in the estimation of φ S by both PCP and LGS with respect to its actual value.

2 Fig. 10 .
Fig. 10.Result of 2000 Monte Carlo repetitions evaluating the MSE in the estimation of ρ, f d and φ S by both PCP and LGS with respect to the sample size.The range and phase terms ρ and φ S were randomized as specified in Table I, while f d ∼ U ([−200, −10)∪[10, 200)), in order to avoid the wildly irregular behavior of the estimators when f d ∈ [−10, 10), shown in Fig. 8.For reference and comparison, the convergence rates given by the inverse Fisher information matrix we present in[26] are portrayed by lines with the corresponding slope adjusted to fit the data.
we report the average performance for randomized parameter values for PCP and LGS in the estimation of the clock synchronization and ranging parameters as a function of either SNR in and SNR out , when SNR out = 30 dB and SNR in = 40 dB, respectively.Similarly as in Section IV-B, we sample f d in the range f d ∈ [−200, −10) ∪ [10, 200), in order to avoid the instability of our estimators when f d ∈ [−10, 10).

2 Fig. 11 .
Fig. 11.Result of 1000 Monte Carlo repetitions evaluating the MSE in the estimation of the clock synchronization and ranging parameters ρ, f d and φ S , by both PCP and LGS with respect to the value of the SNRs.Each SNR was fixed to its maximum when the other was varied.The range and phase terms ρ and φ S were randomized as specified in Table I, while f d ∼ U ([−200, −10) ∪ [10, 200)), in order to avoid the wildly irregular behavior of the estimators when f d ∈ [−10, 10), shown in Fig. 8.

Fig. 12 .
Fig. 12. Example on how the same changes of parameters ∆α ≥ 0 and ∆γ ≥ 0 that exemplified in Fig. 6(b) that the model Y θ with σv = 0 was not identifiable yield different supports of P [n], i.e., supp ∆α and supp ∆γ , when σv > 0. Accordingly, changes in α and γ, however small, will lead to different distributions of Y [n].
Assume that if Y ∼ Y θ for some given θ, Y has PDF f Y (y; θ).Then, Y θ is an identifiable model, and θ is an identifiable parameter, if and only if 1 (Identifiability): (From [29, Definition 11.2.2, p. 523]) Let Y θ be a statistical model with parameter θ ∈ Ω.

TABLE I VALUES
FOR THE PHYSICAL AND SIMULATION PARAMETERS IN OUR NUMERICAL RESULTS, UNLESS OTHERWISE STATED.