Computer Science
ComputerScience,1970-2008-Computer Science-2007.jsonl.gz.gz
Computer Science
fdaa156216484ea6062bbd95a904282b8becc47d
doi:10.1155/2007/57054 Research Article Distributed Time Synchronization in Wireless Sensor Networks with Coupled Discrete-Time Oscillators

In wireless sensor networks, distributed timing synchronization based on pulse-coupled oscillators at the physical layer is currently being investigated as an interesting alternative to packet synchronization. In this paper, the convergence properties of such a system are studied through algebraic graph theory, by modeling the nodes as discrete-time clocks. A general scenario where clocks may have different free-oscillation frequencies is considered, and both time-invariant and time-variant network topologies (or fading channels) are discussed. Furthermore, it is shown that the system of oscillators can be studied as a set of coupled discrete-time PLLs. Based on this observation, a generalized system design is discussed, and it is proved that known results in the context of conventional PLLs for carrier acquisition have a counterpart in distributed systems. Finally, practical details of the implementation of the distributed synchronization algorithm over a bandlimited noisy channel are covered.


INTRODUCTION
Distributed timing synchronization refers to a decentralized procedure that ensures the achievement and maintenance of a common time-scale (frequency and phase) for all the nodes of the network [1]. This condition enables a wide range of applications and functionalities of a sensor networks, including complex sensing tasks (distributed detection/estimation, data fusion), power saving (all nodes sleep and wake-up at coordinate times), and medium access control for communication (e.g., time division multiple access and cooperative communications).
Conventional design of distributed algorithms for timing synchronization prescribes the exchange of local time information through packets carrying a time-stamp to be appropriately elaborated by the transmitting and receiving nodes [1]. Packet-based synchronization has been widely studied, especially in the context of wireline networks. However, the specific features and requirements of wireless sensor networks call for alternative methods that improve both the computational complexity (and therefore energy efficiency) and scalability. Toward this goal, physical layer-based synchronization protocols are currently being investigated that exploit the broadcast nature of radio propagation. The idea is to build distributed algorithms based on the exchange of pulses at the physical layer, thus avoiding the need to perform complex processing at the packet level.
Physical layer-based synchronization was studied in [2] using a mathematical framework developed in [3] in order to model the spontaneous establishment of synchronous periodic activities in biological systems, such as the flashing of fireflies. In [2,3], nodes are modeled as integrate-and-fire oscillators coupled through the transmission of pulses. Convergence is proved under the assumption of an all-to-all interconnection among the nodes. The model was later extended in [4], by explicitly including constraints on the transmission range of each node. In particular, the authors derived a bound on the velocity of convergence by using algebraic graph theory [5]. An implementation of distributed synchronization on a real sensor network testbed was reported in [6]. A related work is [7], where a generalization of the model in [3] is proposed and the regime of an asymptotically dense network is investigated. As a final remark, it should be noted that the framework of physical layer-based timing synchronization has been recently interpreted as a means to achieve distributed estimation/detection [8,9] or data fusion [10].

EURASIP Journal on Wireless Communications and Networking
In this paper, we reconsider physical layer-based synchronization by modeling the sensors as coupled discrete-time oscillators. Basically, each node modifies its current clock based on a weighted average of the residual differences of timing phases as measured with respect to other nodes. The synchronization algorithms proposed in [11] in the context of interbase station communication and [12] for intervehicle transmission can be seen as instances of this general model. The analytical framework is at the same time a generalization and an application of the literature on discrete-time consensus problems for networks of agents (see, e.g., [13]). In particular, differently from [13], here we address the case of clocks with generally different free-oscillation frequencies, and account for the specific features of a wireless network, namely channel reciprocity and randomness (fading). Analysis of convergence of the synchronization process is carried out by algebraic graph theory as in [4], allowing to relate global convergence properties to the local connectivity of the network. The results are first derived for a time-invariant scenario, and then extended to the case where the network topology (or fading) varies with time, building on the results presented in [14].
A central contribution of this paper is the observation that the distributed synchronization system at hand can be modeled as a set of coupled discrete-time phase locked loops (PLLs). The system can thus be seen as a discrete-time version of the network synchronization scheme of [15], that is based on continuously-coupled analog PLLs. This fact allows us to generalize the system design by introducing the concept of loop order. Moreover, we prove that known results about the convergence of conventional PLLs for carrier acquisition have a counterpart in distributed systems. In particular, it is shown that, under appropriate conditions on the interconnections between sensors, (i) a system of first-order distributed PLLs is able to recover perfectly a phase mismatch among the clocks; (ii) in case of a frequency error, first-order loops are able to recover the frequency gap, but at the expense of an asymptotic phase mismatch; (iii) this asymptotic phase mismatch can be reduced by considering second-order loops.
Finally, the analysis is complemented by addressing the issue of a practical implementation of the distributed synchronization algorithm over a bandlimited Gaussian channels.

SYSTEM MODEL AND MAIN ASSUMPTIONS
Let the wireless network be composed of K sensors, where each node, say the kth, has a discrete-time clock with period T k . If the nodes are left isolated, the timing clock of the kth sensor evolves as t k (n) = nT k +τ k (0), where 0 ≤ τ k (0) < T k is an initial arbitrary phase and n = 1, 2, . . . runs over the periods of the timing signal. Two synchronization conditions are of interest. We say the K clocks are frequency synchronized if for each k and for sufficiently large n, where 1/T is the common frequency. A more strict condition requires full frequency and phase synchronization 1 : We remark that the network is said to fractionate into, say, two clusters of synchronization if there exist a permutation function on the nodes' labels, π(i) : [1, . . . , n] → [1, . . . , n] such that for n large enough where the number of nodes in the two clusters is r and K − r, respectively. The definition above generalizes naturally to more than two clusters. Towards the goal of achieving synchronization, the clocks of different sensors can be coupled by letting any node radiate a timing signal as the one sketched in Figure 1. A pulse 2 is transmitted at times t k (n) by the kth node and received through independent flat fading channels by the other sensors. It is assumed that all the nodes transmit with the same power, and that the power P ki received on the wireless link between the ith and the kth user reads where C is an appropriate constant that depends on the transmitted power (assumed here to be the same for all nodes), d ki (n) = d ik (n) is the distance between node i and node k at the nth period, G ki (n) is a random variable accounting for the fading process, and γ is the path loss exponent (γ = 2 ÷ 4). Notice that the fading channel is reciprocal (all transmissions use the same carrier frequency), which implies that G ik (n) = G ki (n) and P ik (n) = P ki (n) for i / = k [16]. As detailed in the following, each node (at any period n) processes the received signal in order to estimate the time difference between its own clock t k (n) and the corresponding "firing" instant of other nodes, that is, t i (n)−t k (n), i / = k, and, based on this measure, it updates its own clock.

The synchronization algorithm
In this section, we consider the synchronization procedure under the ideal assumptions that any node, say the kth, is able to measure exactly the time differences t i (n) − t k (n) and the powers P ki (n) of other nodes (i / = k) based on the received signal. This model is elaborated upon in the first part of the Figure 1: Clock t k (n) of the kth node. τ k (n) is the timing phase in the nth period of the clock.
paper. A practical implementation of the system that alleviates the said assumptions (and in particular, does not require estimation of time of arrivals) is then discussed in Section 7. At the nth period, the kth node updates its clock t k (n) according to a weighted sum of timing differences Δt k (n + 1) 3 : where ε is the step-size (0 < ε < 1) and the coefficients α ki (n) are selected so that α ki (n) ≥ 0 and K i=1, i / =k α ki (n) = 1. The updating rule (5) generalizes the algorithms of [11,12] (and the consensus algorithms, see, e.g., [13]) to a frequencyasynchronous scenario. In this paper, we focus on the following choice for the coefficients α ki (n): The selection of the weighting coefficients (6) is inspired by the algorithms proposed in [11,12]. The rationale of this design is that time differences measured over more unreliable (i.e., low-power) channels should be weighted less when updating the clock, thus rendering the algorithm robust against measurement errors over the fading channels (see also Section 7). Notice that by using (5b) we are implicitly neglecting the propagation delays among nodes, that are assumed to be smaller than the timing resolution. A method to handle propagation delays is described in [11]. As a final remark, we notice that the dynamic system (5) updates the clock t k (n+1) as a convex combination of the times {t i (n)} K i=1 [14].
By defining the vector containing the clocks of all nodes as t(n) = [t 1 (n) · · · t K (n)] T and the vector of clock periods T = [T 1 · · · T K ] T , we can express (5) as the difference vector equation where A(n) is a K × K matrix such that we have [A(n)] ii = Notice that even though we assume channel reciprocity, matrix A(n) is not symmetric. Moreover, by construction, matrix A(n) is nonnegative and stochastic since the sum of the elements on each row sums to one, or equivalently

TIME-INVARIANT FREQUENCY-SYNCHRONOUS NETWORK
In this section, we study the convergence properties of the distributed synchronization algorithm (5) under the following assumptions: (i) frequency-synchronous network, that is, all the clocks share the same period T = T 1 = · · · = T K ; (ii) the network is time-invariant, that is, P ki (n) = P ki for any n and k / = i. From assumption (i), the clock of the kth node can be expressed as where τ k (n) is the timing phase 0 ≤ τ k (n) < T of the kth node in the nth period (see Figure 1). Moreover, by substituting (9) into (5a) and using assumption (ii), it easily follows that the synchronization algorithm (5) can be written in terms of the phases τ k (n) as τ k (n + 1) = τ k (n) + ε · Δτ k (n + 1), (10a) with coefficients α ki : Finally, by defining the vector containing the timings of all nodes as τ(n) = [τ 1 (n) · · · τ K (n)] T , the vector model (7) becomes where A is a K × K matrix such that we have [A] ii = 1 − ε on the main diagonal and [A] i j = ε · α i j for i / = j. Model (12) resembles the one considered in the literature on multiagent coordination (see, e.g., [13]). The goal of this section is to determine the conditions under which the system (12) converges to a unique cluster or to multiple clusters of synchronization for a fixed realization of the fading variables G ki in (4), that is, matrix A is assumed to be deterministic. We will define the conditions of convergence in terms of the properties of the graph associated to the wireless network under study, or equivalently in terms of the system matrix A.

The associated graph and useful definitions
The synchronization algorithm defines a weighted directed graph G = (V, E , A) of order K on the sensor network, where V ={1, . . . , K} is the set of nodes and E ⊆ V × V is the set of edges weighted by the off-diagonal elements of the K × K adjacency matrix [A] i j = α i j . The edge connecting the ith and the jth nodes, i / = j, belongs to E if and only if α i j > 0. Notice that the graph is directed (α i j / = α ji for i / = j), even though fading links are reciprocal (P i j = P ji for i / = j). Moreover, notice that the system matrix reads where L is the graph Laplacian of the network that is defined as [13]: The main result of this section (Theorem 1) relates the convergence properties of the distributed synchronization procedure in (10) with the connectivity of the graph G associated to the sensor network. We need the following definitions.

Definition 1.
A graph G is said to be strongly connected if there exists a path (i.e., a collection of edges in E ) that links every pair of nodes.
It can be proved that strong connectivity of graph G is equivalent to the irreducibility of matrix A [18].

Definition 2.
A K × K matrix A is said to be reducible if there exists a K × K permutation matrix P and an integer r > 0 such that where The degree of irreducibility of a matrix A, or equivalently of strong connectivity of the associated graph G, can be measured by the following quantity where the minimum is taken over all nonempty proper sub- It can be shown that σ = 0 if and only if the matrix A is reducible, or the associated graph G is not strongly connected [19]. 4

Convergence properties
The main result of this section can be now stated as follows.

Theorem 1. (i)
The distributed synchronization (10) converges to a unique cluster of synchronized nodes, τ 1 (n) = · · · = τ K (n) = τ * for n → ∞, if and only if the associated weighted directed graph G is strongly connected, or equivalently if system matrix A is irreducible. (ii) In this case, the system (12) converges to (for n → ∞) An immediate consequence of Theorem 1 is that the timing vectors converge to the average of their initial values τ(0) if and only if the system matrix A is doubly stochastic (i.e., if A T is stochastic as well). In fact, in this case A T 1 = 1 and vector v in (16) reads v = 1/K·1. This condition occurs in balanced networks [13], where i / = j α i j = 1 = i / = j α ji . In sensor networks, this result is of interest in applications where the steady state value of synchronization is used in order to infer the status of the process monitored by the sensor [8,9,20].
Proof. The proof of part (i) of Theorem 1 is available in the literature for applications where the graph G associated to the dynamic system (12) is undirected [5]. In the case of a directed graph, strong connectivity can generally be proved to be only a sufficient condition for synchronization. However, in a wireless fading case with reciprocal channels, the result can be proved as shown in the following. The second part (ii) of Theorem 1 follows from a result derived, among the others, in [13].
As explained above, in order to prove Theorem 1, we only need to show that strong connectivity is also a necessary condition for synchronization. As a by-product, the proposed proof brings insight into the formation of multiple clusters of synchronization (3). Let us assume that A is reducible (or equivalently the associated graph G is not strongly connected). Then, by definition, there exists a permutation matrix P and an integer r > 0 such that (14) holds. But if α i j = 0 in A, then for reciprocity P i j = P ji = 0 and then α ji = 0 (i / = j). This property is sometimes referred to as bidirectionality of the graph (i.e., α i j = 0 if and only if α ji = 0 but α i j and α ji need not to be equal [14]). Therefore, the (14) has all zero entries. Since the permuted matrix P T AP is nonnegative and stochastic, so are submatrices B and D. By applying the permutation function π(k) = P k [1 · · · K] T , where P k is the kth row of matrix P, to the nodes' labels, we can write the system (12) as where τ(n) = Pτ(n). Therefore, the set of r nodes {π(1), . . . , π(r)} evolves independently from the remaining nodes {π(r + 1), . . . , π(K)}. Now, if either B or D are reducible, the reasoning above can be iterated bringing to the formation of multiple independent set of nodes evolving separately. At the end of this procedure, the system matrix can be O. Simeone and U. Spagnolini Figure 2: The rectangular topology considered in the example in Section 3.3.
written as a block matrix with irreducible stochastic blocks on the diagonal. Without loss of generality, let us then assume that B and D are irreducible. From the first part of the proof (see also Appendix A), it follows the two cluster of r and (K − r) nodes synchronize among themselves according to (3). Moreover, the steady state values of the timing vectors depend on the left eigenvectors of B and D according to (16): where is the r × 1 vector collecting the first r entries of τ(n) and The convergence of the dynamic system at hand could be also studied in terms of the subdominant eigenvalue of matrix A, similarly to approach commonly adopted in the context of the analysis of Markov chains [21]. In particular, the following results can be proved relating convergence to the multiplicity of eigenvalue 1.

Theorem 2.
The distributed synchronization (10) converges to a unique cluster of synchronized nodes as in (2) if and only if the subdominant eigenvalue λ 2 / = 1.
Proof. By recalling Theorem 1, it is enough to prove that: (i) if λ 2 = 1, then the graph is not strongly connected; (ii) if the graph is not strongly connected then, λ 2 = 1. Part (i) can be proved similarly to [13]; however, in Appendix A we give an alternative proof based on the measure σ in (15) of irreducibilty of A. Part (ii) does not hold in general for problems with directed graphs but it is easily shown under the reciprocity assumption similarly to Theorem 1.

Numerical results
Here, we present a numerical example to corroborate the analysis discussed above. A network of K = 4 nodes is considered where the nodes are divided into two groups,  (11), the performance is not affected by the value of C in (4) and it only depends on relative distances. Figure 3 shows the timing vector

τ(n) versus n.
After a transient where the nodes tend to synchronize in pairs within the two groups, the system reaches the steady state to the average value τ * /T = 0.475, as stated in Theorem 1, since the system matrix is easily shown to be doubly stochastic for this specific example.
In order to quantify the rate of convergence, from Theorem 2, we notice that the convergence of the synchronization protocol (10) depends on the subdominant eigenvalue λ 2 . In particular, as it is well known from the theory of linear difference equations, the rate of convergence is ruled by a term proportional to |λ 2 | n (see, e.g., [22]). If we define a threshold λ o , we could say that the protocol reaches the steady state condition at the time instant n o for which |λ 2 | no = λ o : n o = log λ o / log |λ 2 |. Therefore, we can take v = − log λ 2 (19) as a measure of the rate of convergence of the algorithm. Figure 4 shows the rate of convergence v versus the normalized distance D/d for ε = 0.3, 0.7. As expected the rate v decreases with increasing D/d and decreasing ε. Along with v, Figure 4 shows the measure of irreducibility (or strong connectivity) σ (15) as dashed lines. It is interesting to note that the rate of convergence v and the measure of irreducibility σ have the same behavior as a function of D/d and ε. This confirms that convergence is strictly related to the connectivity properties of the associated graph, as proved in Theorem 1.

Effect of fading: an example
In this section, the effect of fading on the rate of convergence v is investigated via simulation for linear, ring, and star topologies (see Figure 5). Rayleigh fading is assumed, that is, the fading amplitude G ki in (4) accounts for Rayleigh fading with unit average power. Fading is assumed to be constant for any n during the evolution of the algorithm.  three networks versus the number of nodes K (ε = 0.3). Notice that for K = 2 the three networks coincide and recall that only relative distances are of concern for the behavior of the system (10). As it is expected, the star topology has the largest rate of convergences whereas the linear network yields the slowest convergence.

TIME-VARYING FREQUENCY-SYNCHRONOUS NETWORK
Here, we reconsider the performance of the synchronization procedure (5) by removing the assumption of timeinvariance underlying the analysis of the previous section. However, we still assume a frequency-synchronous network. Overall, the system (5) can be written in vector form in terms of phases as (recall (12)) with the definition of the system matrix A(n) in Section 2.
In this case, the sensor network can be described by a sequence of directed graphs G(n) = (V, E (n), A(n)) of order K, defined similarly to Section 3. In particular, A(n) is the adjacency matrix [A(n)] i j = α i j (n) and the edge connecting the ith and the jth nodes, i / = j, belongs to E (n) if and only if α i j (n) > 0. At each time instant n, the dynamic system describing the synchronization process evolves as where A(n) = I−εL(n) with L(n) being the graph Laplacian at time n (see Section 3). Study of convergence of a family of algorithms encompassing (10) has been recently attempted in a few works (see [13,14] and references therein). In particular, adapting a result first presented in [14] to our case, we are able to relate the convergence of dynamic system (20) to the connectivity properties of the associated sequence of graphs G(n). We need the following definition.

Theorem 3. The distributed synchronization
Proof. Theorem 3 can be proved by specializing the proof of [14, Theorem 3] to our scenario. The basic idea is to exploit the convexity, or the contractive property, of transformation (10). An interesting remark is that reciprocity of fading plays here a key role as it did in the proof of Theorem 1 for the case of fixed topology. Reciprocity of fading translates into bidirectionality of the associated graphs (see Section 3). As proved in [14], in presence of unidirectional communication among nodes (i.e., nonreciprocal fading in our scenario), convergence of synchronization is not necessarily guaranteed if each sensor communicates to every other sensor (either directly or via intermediate nodes) in an interval [n 0 , ∞). On the contrary, in order to guarantee convergence in a unidirectional graph, a limit should be imposed on the time it takes for the graph to become strongly connected, that is, the interval in Theorem 3 should be modified as [n 0 , n 0 + T] where T ≥ 0 finite.

FREQUENCY-ASYNCHRONOUS NETWORK
In the previous sections, it was assumed that all the nodes have the same clock period T (frequency synchronous network). However, in practice, different nodes might have different frequencies {1/T k } K k=1 , and the question arises of whether or not the physical layer-based scheme (5) is still able to achieve synchronization on a strictly connected graph. For a time-invariant scenario (i.e., P ki (n) = P ki for any n and k / = i), it will be shown below that, in presence of a frequency mismatch, the scheme (10) is able to synchronize the clock periods of the nodes (recall (1)), but not their timing phases, so that the full synchronization condition (2) is not achieved.
In this regard, it should be noted that, while perfect synchronization (2) is necessary for many applications, in other scenarios having nodes with synchronized frequency is the only requirement (i.e., to ensure equal sensor duty cycles).
For a frequency-asynchronous time-invariant network, the considered synchronization scheme (5) reads or, in vector form (see (7)): Let us now denote a possible common value for the clock period of all nodes as T (to be determined) as in (1). It follows that the clock of the kth sensor can then be written for sufficiently large n as or equivalently, in vector form, as t(n) = nT · 1 + τ(n) with τ(n) = [τ 1 (n) · · · τ K (n)] T . We are interested in determining if such common frequency 1/T exists and, if so, whether eventually the phases τ(n) converge to the same value for n → ∞. The main conclusion is summarized in the theorem below, whose proof is inspired by the analysis of the convergence of coupled analog oscillators in [23].

Theorem 4.
With reference to (23), under the assumption that the graph G is strictly connected, the system (21) synchronizes the clocks of the K nodes to the common period where v is the normalized left eigenvector of matrix A corresponding to eigenvalue 1: However, the timing phases τ(n) remain generally mismatched and given for n → ∞ by with (·) † denoting the pseudoinverse and the definitions The theorem above states that, in presence of a frequency mismatch, the algorithm (21) is able to synchronize the frequencies of different nodes to the common clock period T in (24). However, the system does not lead to phasesynchronous clocks, and the phase error is determined by the frequency (period) mismatch ΔT according to (25). Notice that, if the network is such that the system matrix A is doubly stochastic (as in the example of Section 3.3), the eigenvector v reads 1/K and the common period T is in this case the average T = 1/K K k=1 T k . Moreover, with doubly stochastic matrix A, condition (25) simplifies as η = 1/K K k=1 τ k (0) since 1 T L † = 0 (see proof below for further details). Finally, we remark that if the frequency mismatch is ΔT = 0 (or equivalently T k = T), Theorem 4 follows from Theorem 1.
Proof. Under the assumption of a connected graph (or irreducible matrix A), according to Theorem 2 or [13], the Laplacian L is easily shown to have rank K − 1, where the one-dimensional null subspaces are defined by the relationships Using the latter equality, recalling (13) and the definition of common clock period T and phases τ(n) in (23), the vector difference equation (22) can be written as An equilibrium state τ(n + 1) = τ(n) = τ * for the difference equation (29) satisfies τ(n + 1) − τ(n) = 0, which yields the condition From (28), it follows that (i) in order for (30) to be feasible (i.e., for an equilibrium point to exist), the common clock period T must satisfy v T ΔT = 0 or equivalently (24); (ii) an equilibrium phase vector τ * must read τ * = (L † /ε)ΔT+η1 where η is an arbitrary constant. It remains to show that the system actually converges for n → ∞ to the equilibrium point τ * determined above, and to evaluate the constant η. Toward the goal of studying convergence, let us define τ (n) = τ(n) − (L † /ε)ΔT. With this change of variables, the difference equation (29) boils down to

EURASIP Journal on Wireless Communications and Networking
Timing error detection Loop filter where we used the relationship LL † ΔT = ΔT, which easily follows from the definition of pseudoinverse and (24). As a consequence of (31), as per Theorem 1, we have τ (n)→ v T τ (0). This expression is equivalent to (25), thus proving the theorem. Notice that from (31) the rate of convergence is the same as in the case of no frequency mismatch.
As a final remark, we notice that the study of timevarying frequency-asynchronous networks is a challenging task and is left for future work.

DISTRIBUTED TIME SYNCHRONIZATION AS COUPLED DISCRETE-TIME PLLs
The purpose of this section is to discuss the distributed synchronization algorithm investigated throughout the paper by casting it into the framework of discrete-time phase locked loops (PLLs) [24]. The discussion is not only be beneficial for a better understanding of the system, but it also allows us to generalize the system design. In order to appreciate the similarity with a discrete-time PLL, Figure 7 depicts the synchronization procedure (5) carried out at each sensor as a linear dynamic feedback system. Similarly to a discrete-time PLL, the adder at the input evaluates the timing error Δt k (n), that is then multiplied by ε and then fed to a voltage controlled clock (VCC) that updates the clock according to (5a). The constant ε plays the role of the loop filter in a discrete-time PLL. The procedure (5) can then be interpreted as a firstorder discrete-time PLL since the loop filter is a trivial pure gain [25]. From the discussion above, it is clear that the second or third order discrete PLLs 5 can be obtained by introducing a loop filter ε(z), with one or two poles respectively, instead of the constant ε in the synchronization system of Figure 7. For instance, in the case of a second-order loop, we can introduce a pole μ in the loop by setting ε(z) = ε/(1 − μz −1 ) with 5 As discussed in [25], loops are never built with order larger than three. 0 < μ < 1. The corresponding update rule (5a) modifies as t k (n + 1) = t k (n) + ε · Δt k (n + 1) The updating rule (32) essentially corrects the local period T k by the estimate of the common clock period t k (n) − t k (n − 1). The convergence analysis of the second-order loop (32) can be carried out similarly to the previous section where a firstorder loop was considered. In particular, the following results hold.

Theorem 5.
If the network of distributed PLL is strictly connected and the system (32) converges, then it synchronizes the clocks of the K nodes to the common period (24). However, under the same conditions, the timing phases τ(n) remain generally mismatched and given for n → ∞ by with (·) † denoting the pseudoinverse, and with definitions and (27).
Proof. The proof is along the lines of the proof of Theorem 4 (see Appendix B for details).
Comparing the statement of the previous theorem with the results derived for a first-order loop (Theorem 4), it can be seen that introducing a pole μ in the loop causes a reduction in the steady state phase error by a factor 1−μ. However, this reduction comes at the expense of decreased margins of stability. In fact, convergence cannot be guaranteed for all values of 0 < ε < 1 and 0 < μ < 1. We refer to Appendix C for further analysis on this point. Here, we illustrate this issue by means of an example. Consider a network with two nodes. In this case, we have α 12 = α 21 = 1 and the graph is connected. Figure 8 shows the four eigenvalues of the system matrix associated with (32) (see Appendices B and C for further details): for different values of the pole μ and ε = 0.9. Notice that the system matrix (35) is 4×4 since (32) is a system of two second order difference equations [22]. Moreover, one eigenvalue of A is 1 irrespective of the value of μ. The absolute value of the remaining eigenvalues tends to one for μ → 1, showing that increasing the value of the pole leads to lack of stability of the equilibrium point (33). Moreover, the value of μ at which a couple of eigenvalues acquire a nonzero imaginary part can be calculated exactly as a function of the spectrum of matrix A, as shown in Appendix C (see (C.5)). In order to corroborate the conclusions above, Figure 9 shows the standard deviation ξ(n) of the timing vector t(n) O. Simeone and U. Spagnolini To conclude, it is interesting to revisit the results of Theorems 1, 4, and 5 in the light of the analogy with conventional PLLs drawn above. It has been shown that in a strictly connected network: (i) a phase error is perfectly recovered for n → ∞ by the distributed synchronization algorithm (5) (Theorem 1); (ii) a frequency error is perfectly recovered at the expense of a phase mismatch for n → ∞ (Theorem 4); (iii) the residual phase mismatch caused by a frequency error can be reduced by introducing a pole in the control loop (Theorem 5). All these results can be read as the counterpart of known facts in the analysis of linearized PLLs, which assert that a first-order loop is indeed able to (i) recover phase errors and (ii) to achieve a constant phase error (referred to as static phase error) in case of a frequency mismatch [25]. Moreover, in this second case, it is interesting to notice how the phase errors (25) and (33) depend on the frequency mismatch ΔT exactly as the static phase error of a PLL [25]. Further results on large-scale randomly deployed networks can be found in [26].

IMPLEMENTATION OF DISTRIBUTED COUPLED DISCRETE-TIME OSCILLATORS
In the previous sections, it was assumed that each node is able to measure time differences and powers of other nodes so as to calculate the phase update Δt k (n). Here, we remove this assumption by presenting a practical scheme to implement the phase detector over a bandlimited noisy channel. Since the algorithm is based only on instantaneous power measurements by different nodes, it applies to both a time-invariant and time-variant scenario. The scheme is inspired by the proposal in [12]. A carrier frequency is dedicated to the synchronization channel, where each node, say the kth, transmits a bandlimited waveform g(t) (say a square-root raised cosine pulse) centered at times t k (n) with symbol period 1/F s (i.e., the time between peak and first zero). The symbol period 1/F s defines the timing resolution of the system. Each node works in an half-duplex mode and measures the received signal on a interval of duration T k around the current timing instant t k (n). Due to the half duplex constraint and the finite switching time between transmitting and receiving mode, each sensor is not able to measure the received signal in an interval of (unilateral) size θ around the firing instant t k (n). It follows that the observation win- Figure 10(a) illustrates a block diagram of the operations performed at the receiver side by each sensor. The receiver performs baseband filtering matched to the transmitted waveform and than samples the received signal at some multiple L of the symbol frequency F s , that is, LF s with L ≥ 1. Based on the N = LF s T samples received in the nth observation window, the kth node computes the update (21a) in case a first-order loop is employed, or (32) if a second-order loop is considered. Not knowing the exact timings and powers of other nodes, t i (n) and P ki with i / = k, the kth sensor cannot directly calculate the updating term Δt k (n) from (5b) and (6). Instead, it estimates these quantities from the received samples, as explained below. After matched filtering and sampling, the discrete-time baseband signal received by the kth node in the nth time period reads (sampling index m ranges within −N/2 < m ≤ N/2 with m = 0 corresponding to the firing time t k (n) of the kth node): where the average energy per symbol reads E ki = C/(d γ ki · F s ) (recall (4)); β ki denotes the Rayleigh fading coefficient, that is a zero-mean and unit-power complex (circularly symmetric) Gaussian random variable with |β ki | 2 = G ki ; and w(n, m) is the additive Gaussian noise with zero mean and power N 0 . Notice that the sample in the interval −ΔLF s ≤ m ≤ ΔLF s is not measured (i.e., zero) due to the half-duplex constraint and the switching time between receive and transmit mode of node k. A sketch of a possible realization of the received signal (36) is provided in Figure 10(b) using an arbitrary waveform g(t).
A simple estimate of ΔT k (n) can then be obtained as where J is the subset of time instants m ∈ (−N/2, −θLF s ) (θLF s , N/2] for which the received signal |y k (n, m)| 2 is above a given threshold as in [12]. The threshold is a system parameter that can be optimized. Being based solely on instantaneous power measurements (i.e., on samples |y k (n, m)| 2 ), the practical scheme (37) has the advantage that it does not require any a priori knowledge at the nodes about network topology or average received powers. For the example of Section 3.3 with no fading (β ki = 1 for every i and k), Figure 11 shows the standard deviation ξ(n) of the timing vector t(n) versus n, where ξ 2 (n)

and expectation E[·]
is taken with respect to noise. We are considering equal clock periods T k = T = 1 for k = 1, . . . , K. Moreover, it is assumed that all nodes transmit the same power and the signal-tonoise ratio for transmission to the closest node (e.g., from 2 to 1) is set to SNR = E 12 /N 0 = 15 dB. Other parameters are as follows: ε = 0.9; the threshold is set for simplicity to zero (see discussion below); distances satisfy D/d = 2; the normalized timing resolution is 1/F s = 0.01; the waveform g(t) is a raised cosine with roll-off δ = 0.2; the switching time is set to θ = 1/F s . 6 We first consider the first-order loop (21a) (or equivalently (32) with pole μ = 0) with different oversampling factors L = 1, 2, 5, 10, 15. It can be seen that the finite resolution of the system produces a performance floor for increasing n, that can be lowered by increasing the oversampling factor L. In any case, an upper bound on the accuracy of synchronization is set by the finite switching time θ = 0.01. This bound is reached for n and L sufficiently large. 7 Adding a pole in the loop as in (32) can increase the convergence speed as shown in Figure 11 for μ = 0.2, 0.4, 0.6. Notice that convergence speed could also be improved by setting an appropriately chosen threshold in (37) (not shown). Finally, Figure 11 shows that an upper bound on the performance of the practical implementation discussed here is set by the performance of the dynamic system (10), where the performance gap is due to noisy observations, finite resolution and finite switching time.

CONCLUSIONS
An increasing number of applications of sensor networks requires the availability of a common time reference to all the nodes. Aiming at scalability and complexity reduction, physical-layer based synchronization qualifies as a valid alternative to conventional packet-based synchronization. In this work, the convergence properties of such an approach have been investigated by modeling the network as a set of discrete-time coupled oscillators and relying on the analytical framework of algebraic graph theory. It has been shown that the system can be equivalently studied as a set of distributed discrete-time PLLs. This observation allowed to generalize the design of the synchronization process, and prove that well-established results on conventional PLLs extend naturally to distributed systems. Finally, practical details of the implementation of the synchronization algorithm have been discussed on a bandlimited Gaussian channel.

Lemma 1. Every eigenvalue of A satisfies
Moreover, if G is connected, only one eigenvalue λ = 1.
From (C.3) and (C.4), it can be seen that the eigenvalue λ of A (with the positive sign in (C.3)) corresponding to λ = 1 equals one. Moreover, by substituting λ = 1 in (C.3), this eigenvalue is unique. In order to prove stability of the system at hand, we then need to show that all the other 2K −1 eigenvalues λ(λ, μ) different from λ = 1 are such that | λ(λ, μ)| < 1. In general, it is not easy to find general conditions on ε and μ that guarantee the condition | λ(λ, μ)| < 1 for any λ in (C.4) (except λ = 1). However, a numerical study can be easily carried out from (C.3) and (C.4). For illustration, consider an undirected graph G. In this case, the eigenvalues of matrix A are real and from (C.4) satisfy −1 < λ < 1 (recalling that 0 < ε < 1). The absolute value | λ(λ, μ)| (C.3) (with positive sign, similar result holds for the negative sign) is plotted versus μ and λ in Figure 12. It can be seen that unless μ → 1, the stability condition is satisfied. Moreover, by studying the sign of the polynomial in (C.3), it is also possible to evaluate the value of the pole μ, say μ 0 (λ), such that for μ ≤ μ 0 (λ) the eigenvalues λ corresponding to λ are real and for μ > μ 0 (λ) are complex conjugate. This value is easily found to be μ 0 (λ) = 2 − λ − 2 1 − λ, ( C . 5 ) which is shown in Figure 12.
As a final remark, we notice that, for μ = 1, the absolute value | λ(λ, μ)| equals one irrespective of λ, leading to an unstable system.