---
abstract: 'The problem of clock offset estimation in a two-way timing exchange regime is considered when the likelihood function of the observation time stamps is exponentially distributed. In order to capture the imperfections in node oscillators, which render a time-varying nature to the clock offset, a novel Bayesian approach to the clock offset estimation is proposed using a factor graph representation of the posterior density. Message passing using the max-product algorithm yields a closed form expression for the Bayesian inference problem.'
address: |
    $^{*}$Department of ECE, Texas A&M University, Texas, TX, 77843 USA.\
    $^{\dagger}$Department of Information Engineering, University of Padova, Padova, Italy.
title: 'Time-varying Clock Offset Estimation in Two-way Timing Message Exchange in Wireless Sensor Networks Using Factor Graphs$^{\ddag}$'
---

Clock offset, factor graphs, message passing, max-product algorithm

Introduction {#sec:intro}
============

Clock synchronization in wireless sensor networks (WSN) is a critical component in data fusion and duty cycling, and has gained widespread interest in recent years [@akyildiz:survey]. Most of the current methods consider sensor networks exchanging time stamps based on the time at their respective clocks [@sadler:survey]. In a *two-way* timing exchange process, adjacent nodes aim to achieve pairwise synchronization by communicating their timing information with each other. After a round of $N$ messages, each node tries to estimate its own clock parameters. A representative protocol of this class is the timing-sync protocol for sensor networks (TPSNs) which uses this strategy in two phases to synchronize clocks in a network [@ganeriwal:tpsn].

The clock synchronization problem in a WSN offers a natural statistical signal processing framework [@serpedin:survey]. Assuming an exponential delay distribution, several estimators were proposed in [@ghaffar:time-out]. It was argued that when the propagation delay $d$ is unknown, the maximum likelihood (ML) estimator for the clock offset $\theta$ is not unique. However, it was shown in [@jeske:ML] that the ML estimator of $\theta$ does exist uniquely for the case of unknown $d$. The performance of these estimators was compared with benchmark estimation bounds in [@eddie:novel]. A common theme in the aforementioned contributions is that the effect of possible time variations in clock offset, arising from imperfect oscillators, is not incorporated and hence, they entail frequent re-synchronization requirements.

In this work, assuming an exponential distribution for the network delays, a closed form solution to clock offset estimation is presented by considering the clock offset as a random Gauss-Markov process. Bayesian inference is performed using factor graphs and the max-product algorithm.

System Model {#sec:sys}
============

By assuming that the respective clocks of a sender node $S$ and a receiver node $R$ are related by $C_R(t)=\theta + C_S(t)$ at real time $t$, the two-way timing message exchange model at the $k$th instant can be represented as [@ghaffar:time-out] [@jeske:ML] $$U_k=d+\theta+X_k, \quad V_k=d-\theta+Y_k\label{U:V}$$ where $d$ represents the propagation delay, assumed symmetric in both directions, and $\theta$ is offset of the clock at node $R$ relative to the clock at node $S$. The network delays, $X_k$ and $Y_k$, are the independent exponential random variables. By further defining $\xi\overset{\Delta}{=}d+\theta$ and $\psi\overset{\Delta}{=}d-\theta$, the model in can be written as $$U_k=\xi+X_k, \quad V_k=\psi+Y_k$$ for $k=1,\ldots,N$. The imperfections introduced by environmental conditions in the quartz oscillator in sensor nodes result in a time-varying clock offset between nodes in a WSN. In order to sufficiently capture these temporal variations, the parameters $\xi$ and $\psi$ are assumed to evolve through a Gauss-Markov process given by $$\xi_k=\xi_{k-1}+w_k, \quad \psi_k=\psi_{k-1}+v_k \quad \text{for}~k=1,\ldots,N$$ where $w_k$ and $v_k$ are $i.i.d$ such that $w_k,v_k\sim\mathcal{N}(0,\sigma^2)$. The goal is to determine precise estimates of $\xi$ and $\psi$ in the Bayesian framework using observations $\{U_k, V_k\}_{k=1}^{N}$. An estimate of $\theta$ can, in turn, be obtained as $$\hat\theta = \frac{\hat\xi-\hat\psi}{2}\label{subs:theta}\;.$$

A Factor Graph Approach {#sec:factor}
=======================

The posterior pdf can be expressed as $$\begin{aligned}
f(\boldsymbol{\xi},\boldsymbol{\psi}|\boldsymbol{U},\boldsymbol{V})&\propto f(\boldsymbol{\xi},\boldsymbol{\psi})f(\boldsymbol{U},\boldsymbol{V}|\boldsymbol{\xi},\boldsymbol{\psi})\nonumber \\
&=f(\xi_0)\prod_{k=1}^{N}f(\xi_k|\xi_{k-1})f(\psi_0)\prod_{k=1}^{N}f(\psi_k|\psi_{k-1})\nonumber \\
&\cdot \prod_{k=1}^{N}f(U_k|\xi_k)f(V_k|\psi_k)\label{pdf:pos}\end{aligned}$$ where uniform priors $f(\xi_0)$ and $f(\psi_0)$ are assumed. Define $\delta_{k-1}^{k}\overset{\Delta}{=}f(\xi_k|\xi_{k-1})\sim\mathcal{N}(\xi_{k-1},\sigma^2)$, $\nu_{k-1}^{k}\overset{\Delta}{=}f(\psi_k|\psi_{k-1})\sim\mathcal{N}(\psi_{k-1},\sigma^2)$, $f_k\overset{\Delta}{=}f(U_k|\xi_k)$, $h_k\overset{\Delta}{=}f(V_k|\psi_k)$, where the likelihood functions are given by $$\begin{aligned}
f(U_k|\xi_k)&=\lambda_{\xi}\exp\left(-\lambda_{\xi}(U_k-\xi_k)\right)\mathbb{I}(U_k-\xi_k)\nonumber\\
f(V_k|\psi_k)&=\lambda_{\psi}\exp\left(-\lambda_{\psi}(V_k-\psi_k)\right)\mathbb{I}(V_k-\psi_k)
\label{pdf:uncons:bayes}
\;.\end{aligned}$$ The factor graph representation of the posterior pdf is shown in Fig. \[fig:fact\]. The factor graph is cycle-free and inference by message passing is indeed optimal. In addition, the two factor graphs shown in Fig. \[fig:fact\] have a similar structure and hence, message computations will only be shown for the estimate $\hat\xi_N$. Clearly, similar expressions will apply to $\hat\psi_N$.

![Factor graph representation of posterior density []{data-label="fig:fact"}](factorGraphFig1.pdf "fig:") ![Factor graph representation of posterior density []{data-label="fig:fact"}](factorGraphFig2.pdf "fig:")

The message updates in factor graph using max-product can be computed as follows $$\begin{aligned}
m_{f_N\rightarrow\xi_N}&=f_N, \quad m_{\xi_N\rightarrow\delta_{N-1}^{N}}=f_N\nonumber \\
%m_{\xi_N\rightarrow\delta_{N-1}^{N}}&=f_N\nonumber& \\
m_{\delta_{N-1}^{N}\rightarrow\xi_{N-1}}&\propto \underset{\xi_N}{\max}~\delta_{N-1}^{N}\cdot m_{\xi_N\rightarrow\delta_{N-1}^{N}}\nonumber\\
&=\underset{\xi_N}{\max}~\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-(\xi_N-\xi_{N-1})^2}{2\sigma^2}\right)
\nonumber \\
&\cdot \exp\left(\lambda_{\xi}\xi_N\right) \mathbb{I}(U_N-\xi_N)\label{m:N}\end{aligned}$$ which can be rearranged as $$\begin{split}
m_{\delta_{N-1}^{N}\rightarrow\xi_{N-1}}\propto \underset{\xi_N \leq U_N}{\max}~\exp &\big( A_{\xi,N}\xi_N^2+B_{\xi, N}\xi_{N-1}^2+ \\
&C_{\xi, N}\xi_N\xi_{N-1} +D_{\xi, N}\xi_N\big)
\end{split}
\label{m:delta:xi:N-1:max}$$ where $$\begin{aligned}
A_{\xi, N}&\overset{\Delta}{=}-\frac{1}{2\sigma^2}, \quad B_{\xi, N}\overset{\Delta}{=}-\frac{1}{2\sigma^2}\nonumber \\
C_{\xi, N}&\overset{\Delta}{=}\frac{1}{\sigma^2}, \quad D_{\xi, N}\overset{\Delta}{=}\lambda_{\xi} \label{const:N}\;.\end{aligned}$$ Let $\bar\xi_N$ be the unconstrained maximizer of the exponent in the objective function above, i.e., $$\begin{split}
\bar\xi_N = \arg \underset{\xi_N}{\max}~\big( & A_{\xi, N}\xi_N^2+B_{\xi, N}\xi_{N-1}^2+C_{\xi, N}\xi_N\xi_{N-1}+  \\
& D_{\xi, N}\xi_N \big) \;.\nonumber \\
\end{split}$$ This implies that $$\bar\xi_N=-\frac{C_{\xi, N}\xi_{N-1}+D_{\xi, N}}{2A_{\xi, N}} \label{bar:N}\;.$$ If $\bar\xi_N>U_N$, then the estimation problem is solved, since $\hat\xi_N = U_N$. However, if $\bar\xi_N \le U_N$, the solution is $\hat\xi_N = \bar\xi_N$. Therefore, in general, we can write $$\hat\xi_N=\min\left(\bar\xi_N,U_N\right)\label{hat:N}\;.$$ Notice that $\bar\xi_N$ depends on $\xi_{N-1}$, which is undetermined at this stage. Hence, we need to further traverse the chain backwards. Assuming that $\bar\xi_{N}\leq U_N$, $\bar\xi_N$ from can be plugged back in which after some simplification yields $$\begin{split}
m_{\delta_{N-1}^{N}\rightarrow\xi_{N-1}}\propto\exp\Bigg\{ &\left(B_{\xi, N}-\frac{C_{\xi, N}^2}{4A_{\xi, N}}\right)\xi^2_{N-1} - \\
& \frac{C_{\xi, N} D_{\xi, N}}{2A_{\xi, N}}\xi_{N-1} \Bigg\} \;.
\end{split}
\label{m:delta:xi:N-1}$$

Similarly the message from the factor $\delta_{N-2}^{N-1}$ to the variable node $\xi_{N-2}$ can be expressed as $$\begin{aligned}
m&_{\delta_{N-2}^{N-1}\rightarrow\xi_{N-2}}\propto  \underset{\xi_{N-1} \leq U_{N-1}}{\max}~\delta_{N-2}^{N-1}\cdot m_{\xi_{N-1}\rightarrow\delta_{N-2}^{N-1}}\nonumber\\
=&\underset{\xi_{N-1}}{\max}~\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(\xi_{N-1}-\xi_{N-2})^2}{2\sigma^2}\right) \nonumber\\
& \cdot\exp\left \{ \left(B_{\xi, N}-\frac{C_{\xi, N}^2}{4A_{\xi, N}}\right)\xi^2_{N-1}-\frac{C_{\xi, N} D_{\xi, N}}{2A_{\xi, N}}\xi_{N-1} \right \}\nonumber\\
&\cdot\exp\left(\lambda_{\xi}\xi_{N-1}\right)\mathbb{I}(U_{N-1}-\xi_{N-1})\nonumber \;.\end{aligned}$$ The message above can be compactly represented as $$\begin{aligned}
m&_{\delta_{N-2}^{N-1} \rightarrow\xi_{N-2}} \propto\nonumber \underset{\xi_{N-1} \leq U_{N-1}} {\max}~\exp(A_{\xi, N-1}\xi_{N-1}^2+ \nonumber \\
&  B_{\xi, N-1}\xi_{N-2}^2+C_{\xi, N-1}\xi_{N-1}\xi_{N-2}+D_{\xi, N-1}\xi_{N-1}) \label{m:delta:xi:N-2:max}\end{aligned}$$ where $$\begin{aligned}
A_{\xi, N-1}&\overset{\Delta}{=}-\frac{1}{2\sigma^2}+B_{\xi, N}-\frac{C_{\xi, N}^2}{4A_{\xi, N}}, \quad \nonumber \\
B_{\xi, N-1}&\overset{\Delta}{=}-\frac{1}{2\sigma^2}, \quad C_{\xi, N-1}\overset{\Delta}{=}\frac{1}{\sigma^2} \\
D_{\xi, N-1}&\overset{\Delta}{=}\lambda_{\xi}-\frac{C_{\xi, N}D_{\xi, N}}{2A_{\xi, N}}\;.\label{const:N-1}\end{aligned}$$ Proceeding as before, the unconstrained maximizer $\bar\xi_{N-1}$ of the objective function above is given by $$\bar\xi_{N-1}=-\frac{C_{\xi, N-1}\xi_{N-2}+D_{\xi, N-1}}{2A_{\xi, N-1}}\label{bar:N-1}$$ and the solution to the maximization problem is expressed as $$\hat\xi_{N-1}=\min\left(\bar\xi_{N-1},U_{N-1}\right)\;.\label{hat:N-1}$$ Again, $\bar\xi_{N-1}$ depends on $\xi_{N-2}$ and therefore, the solution demands another traversal backwards on the factor graph representation in Fig. \[fig:fact\]. By plugging $\bar\xi_{N-1}$ back in , it follows that $$\begin{aligned}
&m_{\delta_{N-2}^{N-1}\rightarrow\xi_{N-2}}\propto\nonumber \\
&\exp\left \{ \left(B_{\xi, N-1}-\frac{C_{\xi, N-1}^2}{4A_{\xi, N-1}}\right)\xi^2_{N-2}-\frac{C_{\xi, N-1} D_{\xi, N-1}}{2A_{\xi, N-1}}\xi_{N-2} \right \}\label{m:delta:xi:N-2}\end{aligned}$$ which has a form similar to . It is clear that one can keep traversing back in the graph yielding messages similar to and . In general, for $i = 1, \ldots, N-1$, we can write $$\begin{split}
A_{\xi, N-i} &\overset{\Delta}{=} -\frac{1}{2 \sigma^2} + B_{\xi, N-i+1} - \frac{C_{\xi, N-i+1}^2}{4 A_{\xi, N-i+1}} \\
B_{\xi, N-i} & \overset{\Delta}{=} -\frac{1}{2 \sigma^2}, \quad C_{\xi, N-i} \overset{\Delta}{=} \frac{1}{\sigma^2} \\
D_{\xi, N - i} & \overset{\Delta}{=} \lambda_{\xi} - \frac{C_{\xi, N - i + 1} D_{\xi, N - i + 1}}{2 A_{\xi, N - i + 1}}
\end{split}
\label{const:N-i}$$ and $$\begin{aligned}
\bar{\xi}_{N-i} &=& -\frac{C_{\xi, N-i} \xi_{N - i - 1} + D_{\xi, N - i}}{2 A_{\xi, N - i}} \label{bar:N-i} \\
\hat{\xi}_{N-i} &=& \min \left( \bar{\xi}_{N-i}, U_{N - i} \right) \label{hat:N-i}\;.\end{aligned}$$ Using and with $i = N-1$, it follows that

$$\bar\xi_{1}=-\frac{C_{\xi, 1}\xi_{0}+D_{\xi, 1}}{2A_{\xi, 1}}, \quad \hat\xi_{1}=\min\left(\bar\xi_{1},U_{1}\right) \label{bar:1}\;.
%\label{hat:1}\;.$$

Similarly, by observing the form of and , it follows that $$m_{\delta_{0}^{1}\rightarrow\xi_{0}}\propto \exp\left \{ \left(B_{\xi, 1}-\frac{C_{\xi,1}^2}{4A_{\xi,1}}\right)\xi^2_{0}-\frac{C_{\xi,1} D_{\xi,1}}{2A_{\xi,1}}\xi_{0} \right \}\;.\label{m:delta:xi:0}$$ The estimate $\hat\xi_0$ can be obtained by maximizing . $$\hat\xi_0=\bar\xi_0=\underset{\xi_{0}}{\max}~m_{\delta_{0}^{1}\rightarrow\xi_{0}}
=\frac{C_{\xi,1}D_{\xi, 1}}{4A_{\xi, 1}B_{\xi, 1}-C_{\xi, 1}^2}\;.\label{hat:0}$$ The estimate in can now be substituted in to yield $\bar\xi_1$, which can then be used to solve for $\hat\xi_1$. Clearly, this chain of calculations can be continued using recursions and .\
Define $$g_{\xi,k}(x)\overset{\Delta}{=}-\frac{C_{\xi,k}x+D_{\xi,k}}{2A_{\xi,k}}\;.\label{g:def}$$

For real \[lem:1\]numbers $a$ and $b$, the function $g_{\xi,k}(.)$ defined in satisfies $$g_{\xi,k}\left(\min(a, b)\right)=\min\left(g_{\xi,k}(a), g_{\xi,k}(b)\right)\;.$$

Proof: The constants $A_{\xi,k}$, $C_{\xi,k}$ and $D_{\xi,k}$ are defined in and . The proof follows by noting that $\frac{-C_{\xi,k}}{2A_{\xi,k}}>0$ which implies that $g_{\xi,k}(.)$ is a monotonically increasing function.\
Using the notation $g_{\xi,k}(.)$, it follows that $$\begin{aligned}
\bar\xi_1&=g_{\xi,1}\left(\hat\xi_0\right), \quad \hat\xi_1=\min\left(U_1,g_{\xi,1}\left(\hat\xi_0\right)\right)\nonumber \\
%\hat\xi_1&=\min\left(U_1,g_{\xi,1}\left(\hat\xi_0\right)\right)\nonumber \\
\bar\xi_2&=g_{\xi,2}\left(\hat\xi_1\right), \quad \hat\xi_2=\min\left(U_2,g_{\xi,2}\left(\hat\xi_1\right)\right)\nonumber
%\hat\xi_2&=\min\left(U_2,g_{\xi,2}\left(\hat\xi_1\right)\right)\end{aligned}$$ where $$\begin{aligned}
g_{\xi,2}\left(\hat\xi_1\right)&=g_{\xi,2}\left(\min\left(U_1,g_{\xi,1}\left(\hat\xi_0\right)\right)\right)\nonumber \\
&=\min\left(g_{\xi,2}\left(U_1\right),g_{\xi,2}\left(g_{\xi,1}\left(\hat\xi_0\right)\right)\right)\label{g:nest}\end{aligned}$$ where follows from Lemma \[lem:1\]. The estimate $\hat\xi_2$ can be expressed as $$\begin{aligned}
\hat\xi_2&=\min\left(U_2,\min\left(g_{\xi,2}\left(U_1\right),g_{\xi,2}\left(g_{\xi,1}\left(\hat\xi_0\right)\right)
\right)\right)
\nonumber\\
&=\min\left(U_2,g_{\xi,2}\left(U_1\right),g_{\xi,2}\left(g_{\xi,1}\left(\hat\xi_0\right)\right)\right)\;.\end{aligned}$$ Hence, one can keep estimating $\hat\xi_k$ at each stage using this strategy. Note that the estimator only depends on functions of data and can be readily evaluated. For $m \geq k$, define $$G_{\xi, k}^m(.)\overset{\Delta}{=}g_{\xi,m}\left(g_{\xi,m-1}\ldots g_{\xi,k}\left(.\right)\right)\;.\label{G:def}$$ The estimate $\hat\xi_N$ can, therefore, be compactly represented as $$\hat\xi_N=\min\left(U_N,G_{\xi,N}^N\left(U_{N-1}\right),\dots,G_{\xi,2}^N\left(U_1\right),
G_{\xi,1}^N\left(\hat\xi_0\right)\right) \label{xi:1}\;.$$ By a similar reasoning, the estimate $\hat\psi_N$ can be analogously expressed as $$\hat\psi_N=\min\left(V_N,G_{\psi,N}^N\left(V_{N-1}\right),\dots,G_{\xi,2}^N\left(V_1\right),
G_{\xi,1}^N\left(\hat\psi_0\right)\right)$$ and the factor graph based clock offset estimate (FGE) $\hat\theta_N$ is given by $$\hat\theta_N=\frac{\hat\xi_N-\hat\psi_N}{2}\;.\label{theta:fg}$$ It only remains to calculate the functions of data $G(.)$ in the expressions for $\hat\xi_N$ and $\hat\psi_N$ to determine the FGE estimate $\hat\theta_N$. With the constants defined in , it follows that $$G_{\xi,N}^N(U_{N-1})=-\frac{C_{\xi,N}U_{N-1}+D_{\xi,N}}{2A_{\xi,N}}=U_{N-1}+\lambda_{\xi}\sigma^2\;.$$ Similarly it can be shown that $$G_{\xi,N-1}^N(U_{N-2})=U_{N-2}+2\lambda_{\xi}\sigma^2$$ and so on. Using the constants defined in for $i=N-1$, it can be shown that $\hat\xi_0=\frac{C_{\xi,1}D_{\xi,1}}{4A_{\xi,1}B_{\xi,1}-C_{\xi,1}^2}=+\infty$. This implies that $G_{\xi,1}^N(\hat\xi_0)=+\infty$. Plugging this in yields $$\hat\xi_N=\min(U_N, U_{N-1}+\lambda_{\xi}\sigma^2,%U_{N-2}+2\lambda_{\xi}\sigma^2,\nonumber \\
\ldots,U_1+(N-1)\lambda_{\xi}\sigma^2)\;.\label{map:E}$$ Similarly, the estimate $\hat\psi_N$ is given by $$\hat\psi_N=\min(V_N, V_{N-1}+\lambda_{\psi}\sigma^2,%V_{N-2}+2\lambda_{\psi}\sigma^2,\nonumber \\
\ldots,V_1+(N-1)\lambda_{\psi}\sigma^2)\;\label{map:E:psi}$$ and the estimate $\hat\theta_N$ can be obtained using as $$\begin{aligned}
\hat\theta_N=\frac{1}{2} \min(U_N, U_{N-1}&+\lambda_{\xi}\sigma^2,U_{N-2}+2\lambda_{\xi}\sigma^2,\nonumber \\
&\ldots,U_1+(N-1)\lambda_{\xi}\sigma^2)-\nonumber\\\frac{1}{2}\min(V_N, V_{N-1}&+\lambda_{\psi}\sigma^2,V_{N-2}+2\lambda_{\psi}\sigma^2,\nonumber \\
&\ldots,V_1+(N-1)\lambda_{\psi}\sigma^2)\;.\label{example:exp:theta:N}\end{aligned}$$ As the Gauss-Markov system noise $\sigma^2 \rightarrow 0$, yields $$\hat\theta_N\rightarrow \hat\theta_{\textrm{ML}}=\frac{\min\left(U_N, \ldots, U_1\right)-\min\left(V_N, \ldots, V_1\right)}{2}\label{ml}$$ which is the ML estimator proposed in [@jeske:ML].

![Comparison of MSE of $\hat\theta_N$ and $\hat\theta_{\textrm{ML}}$.[]{data-label="Bayes_theta_eps"}](Bayes_theta_vs_k_exp_distr.pdf)

Simulation Results
==================

With $\lambda_{\xi}=\lambda_{\psi}=10$ and $\sigma=10^{-2}$, Fig. \[Bayes\_theta\_eps\] shows the MSE performance of $\hat\theta_N$ and $\hat\theta_{\textrm{ML}}$, compared with the Bayesian Chapman-Robbins bound (BCHRB). It is clear that $\hat\theta_N$ exhibits a better performance than $\hat\theta_{\textrm{ML}}$ by incorporating the effects of time variations in clock offset. As the variance of the Gauss-Markov model accumulates with the addition of more samples, the MSE of $\hat\theta_{\textrm{ML}}$ gets worse. Fig. \[Bayes\_theta\_vs\_sigma\_GM\] depicts the MSE of $\hat\theta_N$ in with $N = 25$. The horizontal line represents the MSE of the ML estimator . It can be observed that the MSE obtained by using the FGE for estimating $\theta$ approaches the MSE of the ML as $\sigma < 10^{-3}$.

![MSE in estimation of $\theta_N$ vs $\sigma$.[]{data-label="Bayes_theta_vs_sigma_GM"}](Bayes_theta_vs_sigma_GM.pdf)

Conclusion
==========

The estimation of a possibly time-varying clock offset is studied using factor graphs. A closed form solution to the clock offset estimation problem is presented using a novel message passing strategy based on the max-product algorithm. This estimator shows a performance superior to the ML estimator proposed in [@jeske:ML] by capturing the effects of time variations in the clock offset efficiently.

[1]{} I. F. Akyildiz, W. Su, Y. Sankarasubramanium, and E. Cayirci, “A survey on sensor networks,” *IEEE Commun. Mag.*, vol. 40, no. 8, pp. 102-114, Aug. 2002. B. Sadler and A. Swami, “Synchronization in sensor networks: An overview,” in *Proc. IEEE Military Commun. Conf. (MILCOM 2006)*, pp. 1-6, Oct. 2006. S. Ganeriwal, R. Kumar, and M.B. Srivastava, “Timing-sync protocol for sensor networks,” in *Proc. SenSys*, Los Angeles, CA, pp. 138-149, Nov. 2003. Y.-C. Wu, Q.M. Chaudhari, and E. Serpedin, “Clock synchronization of wireless sensor networks,” *IEEE Signal Process. Mag.*, vol. 28, no. 1, pp. 124-138, Jan. 2011. H. S. Abdel-Ghaffar, “Analysis of synchronization algorithm with time-out control over networks with exponentially symmetric delays,” *IEEE Trans. Commun.*, vol. 50, no. 10, pp. 1652-1661, Oct. 2002. D. R. Jeske, “On the maximum likelihood estimation of clock offset,” *IEEE Trans. Commun.*, vol. 53, no. 1, pp. 53-54, Jan. 2005. K.-L. Noh, Q. M. Chaudhari, E. Serpedin, and B. Suter, “Novel clock phase offset and skew estimation using two-way timing message exchanges for wireless sensor networks,” *IEEE Trans. Commun.*, vol. 55, no. 4, pp. 766-777, Apr. 2007.
