Abstract

To support system-wide synchronization accuracy and precision in the sub-microsecond range without using GPS technique, the precise time protocol (PTP) standard IEEE-1588 v2 is chosen. Recently, a new clock skew estimation technique was proposed for the slave based on a dual slave clock method that assumes that the packet delay variation (PDV) in the Ethernet network is a constant delay. However, papers dealing with the Ethernet network have shown that this PDV is a long range dependency (LRD) process which may be modeled as a fractional Gaussian noise (fGn) with Hurst exponent () in the range of . In this paper, we propose a new clock skew estimator based on the maximum likelihood (ML) technique and derive an approximated expression for the Cramer-Rao lower bound (CRLB) both valid for the case where the PDV is modeled as fGn (). Simulation results indicate that our new clock skew method outperforms the dual slave clock approach and that the simulated mean square error (MSE) obtained by our new proposed clock skew estimator approaches asymptotically the developed CRLB.

1. Introduction

Recently, the topic of synchronization in network nodes has attracted significant attention for a number of applications [110]. In many cases the global positioning system (GPS) is used to synchronize frequency and time accuracy in the sub-microsecond range. On the other hand, although GPS can offer reliable clock synchronization accuracy, network operators are seeking to minimize the use of GPS within their networks (especially for the indoor places) due to a weak indoor GPS signal [11]. Technologies like synchronous Ethernet (SyncE) [12], network time protocol (NTP) [13], and IEEE 1588v2 (PTP) [14] have been employed to deliver frequency synchronization (GPS, SyncE, NTP, and PTP) and time synchronization (GPS, NTP, and PTP) in Ethernet networks. This provides significant cost savings in network equipment as well as in ongoing installation and maintenance [15].

The SyncE is defined by the ITU [12] as a means of using Ethernet to transfer frequency via the Ethernet physical layer. It recovers frequency from one or more Ethernet signals and delivers timing to outgoing Ethernet signals. However, in order to use SyncE effectively, hardware changes in the switches and routers across the network are needed, which are not yet done.

The NTP was proposed for the first time in 1985 and has been revised several times [16]. The time filtering rules and the servoregulation algorithms implemented by the NTP clients have been designed for large networks; hence they are able to trace time down to a few milliseconds on geographically large-scale network [16]. However, the NTP protocol is not designed for local area networks (LANs) and, due to slow response and software clock implementations (timestamping at application layer), cannot compete with PTP performance even when special implementation are used [17].

IEEE 1588 v2 is also named as PTP, used for time synchronization in networked measurement and control systems. IEEE 1588 is a master-slave time synchronization protocol that targets clock synchronization accuracy of less than 1 μs [18] due to the precise timestamping at the physical layer. In the last decade, the PTP has experienced a large diffusion in systems that require high performance synchronization over Ethernet, since its simplicity and autoconfiguration are suitable for embedded and low cost applications [19].

There are two main functions established by the standard [20]: (i) establishing a master-slave hierarchy of clocks in which each slave synchronizes to its master and (ii) making the necessary information available for the slaves to perform this synchronization. Most of the basic approaches [20, 21, 23] involve clock synchronization that is based on the knowledge of the link propagation delay that is the master-to-slave message delay (uplink) and slave to master message delay (downlink). By assuming symmetric communication links (uplink delay = downlink delay) [20, 21, 23] estimate the clock offset, that is, the time difference between the master and the slave. However, in a realistic network the up- and downlink delays may be very different. Thus, this approach may lead to errors in clock synchronization in a realistic network. Recently, a clock synchronization method (frequency synchronization) was proposed [18] based on dual slave clocks in a slave, suitable for symmetric and asymmetric communication links. This estimation method [18] is based on the knowledge of the one-way delay (OWD) (on the uplink delay) instead on the knowledge of the two-way delay (2WD) (namely, on both directions up- and downlink delays) as introduced in the basic concept of the IEEE 1588 [22]. It does not assume having symmetric communication links (uplink delay = downlink delay) unlike [20, 21, 23]. Thus, its estimation method accuracy [18] does not depend on the up- and downlink symmetry unlike in [20, 21, 23]. Therefore, if the up- and downlink are asymmetric, the estimation method accuracy of [18] will not be affected by this asymmetry unlike how it is in the case of [20, 21, 23]. However, [18] does not take into account the PDV existing in the network. According to [24], PDV is a random delay due to the behavior of switches or routers in the network. The primary source of this random delay is the output queuing delay, caused when a PTP message arrives at a switch or router when the exit port is blocked by other traffic, and the PTP Sync message has to wait in a queue. This random delay (PDV) is characterized according to [2530] as LRD characteristics such as fGn with Hurst exponent in the range of .

According to [31], fGn is a commonly used model of network traffic with LRD. In [32], the authors state that fGN with a single parameter (with ) is a widely used traditional traffic model and refer the reader to [3335]. They [32] also add that traditional methods to synthesize traffic are based on fGN with a single parameter and refer the reader to [3640]. Please note that, in [30, 41], fGN was used to model the network traffic where the work in [26] is a supplement of [41] that supplies experimental investigations to show that fGN can yet be used for modeling autocorrelation functions of various types of network traffic consistently in the sense that the modeling accuracy (expressed by mean square error) is in the order of magnitude of 10−3. The conclusion of [30] was that the analyses of generating fGN data have shown that this generator should be recommended for practical simulation studies of high-speed telecommunication networks (Ethernet, ATM, VBR video traffic, Web traffic, Telnet, and FTP), since it is very accurate.

It should be pointed out that there are also other works [4248] that model the PDV as a Gaussian or exponential delay. These papers [4248] do not need the symmetric assumption on the up- and down link unlike [20, 21, 23]. However, they ([4248]) work with the 2WD messaging technique where both sides (uplink and downlink) are affected by the PDV unlike with the OWD messaging technique. Therefore, in cases where the uplink PDV is lower than the downlink PDV we may get with the OWD messaging method which is affected only by the uplink PDV, better performance compared with the 2WD messaging technique. As already explained earlier in this paper, [18] uses the OWD messaging technique (unlike [4248]) and does not assume a symmetric communication link unlike [20, 21, 23]. But, [18] does not take into account the PDV. Please refer to Table 1 for a brief comparison between the mentioned methods. Our motivation in this paper is to develop a new normalized frequency estimator for the clock skew that works with the OWD messaging technique, does not assume a symmetric communication link, and takes into account the PDV modeled with the fGN.

In this paper we propose a new normalized frequency estimator for the clock skew based on the dual slave clocks technique proposed in [18], where the PDV is modeled as fGn unlike how it was done in [18]. The clock skew is defined as the normalized difference in the clock frequencies between the slave and the master. This new method is based on the ML estimation technique and outperforms the existing method [18] from the MSE point of view in a realistic communication link. In addition we derive in this paper an approximated expression for the Cramer-Rao lower bound (CRLB) depending on the Hurst exponent in the range of . Simulation results indicate that the MSE obtained from our new proposed normalized frequency estimator declines and approaches asymptotically our new derived CRLB.

The paper is organized as follows. The system description is first introduced in Section 2. The proposed ML estimator and the CRLB are presented in Section 3. In Section 4 simulation results are presented and the conclusion is given in Section 5.

2. System Description

2.1. IEEE 1588: Basic Concept

In this subsection we describe the basic concept of IEEE 1588 based on [18]. Figure 1 (recalled from [18] and slightly modified) describes the PTP messaging. The master sends Sync messages at a default rate of once every 2 s. The master clock measures the time when the first Sync message is sent. When a slave clock receives the Sync message, it measures and stores the time of reception . The master then sends a Follow_Up message that contains the actual value of the timestamp . At this point the slave calculates the master-to-slave delay [18]:On the other hand, the master-to-slave delay is the sum of the time difference between the master and the slave clocks (i.e., clock offset) and the master-to-slave message delay [22]:

In order to estimate the clock offset the slave periodically sends the Delay_Req message and stores the transmission time with a timestamp . When the master receives the Delay_Req message, it sends a Delay_Resp message, which conveys the timestamp of the reception time of the request message. The slave calculates the slave-to-master delay from these timestamps as [18]On the other hand, the slave-to-master delay is the sum of the negative time difference between the slave and the master clocks (i.e., clock offset) and the slave-to-master message delay [22]:Assuming that the master-to-slave message delay is equal to the slave-to-master message delay we may write [22]where is the estimated message delay of the uplink or downlink. By using (1) to (5), the slave estimates and the clock offset () asThe slave adjusts its time to minimize the clock offset (), thereby synchronizing with the master clock. Please note that (7) assumes that there is no frequency offset between the master and the slave clocks.

Every oscillator has its unique clock frequency consisting of a nominal frequency in addition to a normalized frequency difference of parts per million (PPM) from the nominal frequency multiplied by this nominal frequency. In addition, every oscillator has a random jitter . Therefore the master and the slave clocks may have different frequencies and thus are not synchronized. For example, suppose that mater's frequency is () and slave's frequency is (); thus their frequencies are not synchronized. In that case, the clock offset () generally keeps increasing and may cause the communication link to fail. Thus, the importance of estimating the normalized frequency difference between the slave and the master that is clock skew ().

2.2. Dual Slave Clocks Method

In this subsection we present briefly the dual slave clock method in a slave as presented in [18]. The basic PTP concept makes the assumption that the up- and downlink delays are equal. However, in a realistic network the up- and downlink delays may be very different. Thus, this approach (where the up- and downlink delays are assumed to be equal) may lead to errors in the clock synchronization task in a realistic network. To eliminate the assumption of symmetric links, [18] proposed two slave clocks, triggering the clock counters in the slaves. Figure 2 (recalled from [18] and slightly modified) describes the PTP messaging required by the dual slave clocks method, where is the time when the th Sync message is sent. and are the times that slave clocks 1 and 2 measure at the receive of the th Sync message, respectively. Notably, the dual slave clocks method requires only one-way time transfer. Clock 2 runs at a normalized clock period of which is normalized to that of the master clock, where is the normalized clock frequency difference between slave clock 2 and the master clock (i.e., normalized clock skew). A D-type flip-flop is employed to divide the frequency by two. Therefore, slave clock 2 runs at a normalized clock period of , which is half of the slave clock 1 running at a normalized clock period of . Both slave clocks have the same offset, which is assumed to be . Thus, according to [18] we may express the times , aswhere(i) and are random jitters of slave clock 1 and 2 within the th period, respectively.(ii)The random jitter is assumed to be an independent Gaussian random variable (RV) with zero mean and variance .(iii)Since slave clock 2 is derived from slave clock 1 and the period of slave clock 2 is twice that of slave clock 1, the variance of is due to accumulation of random jitters at two positive edges of slave clock 1.Here ends the brief presentation of the dual slave clock method in a slave as presented in [18].

The uplink delay can be represented as a constant value with an additive zero mean random noise called PDV. The PDV as already was mentioned results from the behaviour of switches or routers in the packet network [24]. In the following, the constant value is the we saw in [18], and will be denoted by the PDV of the th cycle. Therefore, (8) must be written in the following form:where is assumed to be fractional Gaussian noise (fGn) with zero mean with Hurst exponent and variance .

3. ML Clock Skew Estimation for Fractional Gaussian Noise PDV

ML estimators have a number of attractive characteristics. For example, they are asymptotically unbiased and efficient.

In this section we derive the ML estimation and the CRLB for the clock skew based on the dual slaves clock model presented in [18] where the PDV is modeled with fGn unlike how it was done in [18].

In Section 3.1 we derive the ML estimation for the clock skew under the assumption of fGn PDV with Hurst exponent in the range of .

In Section 3.2 we derive an expression for the Cramer-Rao (CR) lower bound on the variance of the clock skew.

3.1. ML Estimation of Clock Skew

Theorem 1. For the following assumptions(i)the Hurst exponent is in the range of ,(ii),
is assumed to be fGn with zero mean with Hurst exponent and variance .

The ML estimation for the clock skew may be defined as where is defined by denotes Gamma function and is the number of observation.

Comments(i) denotes the difference between two consecutive samples.(ii)The system equations are (9) and (10) instead of (8) in order to describe the system more closely to a real system with PDV.

Proof. Subtracting (10) from (9) yieldsSubtracting (13) of the th Sync message from that of the th Sync message indicates thatwhich can be also written asThus, may be expressed aswhere according to [42]According to [49, 50], the probability density function (PDF) of vector (of length ) may be given bywhere denotes the transpose of and denotes the determinant of the matrix .
is the covariance matrix of the vector :The elements of are given by [50, 51]where, according to [51, 52],Note that so .
can be decomposed as was done in [49, 50]whereAccording to [50] we can writewhereNext, we recall from [50] the expression for the log-likelihood function:The derivative of the log-likelihood function (26) with respect to yieldsNote that [50] has taken the derivative of (26) with respect to . Thus, our final expression is very different.
From (24) and (25) we can writeBy substituting (28) into (27) we obtainAccording to [50], the expression can be written aswhereThen, we substitute (30) into (29) and obtainLet us define for simplicity the next two variables:Inserting (34) and (35) into (33) givesNext, the derivative of the log-likelihood function (36) with respect to is carried out and set to zeroAppendix A shows in detail the derivation of the log-likelihood function (36) and the calculation of (37). In Appendix C we give the sum limitation of (36).
According to Appendices A and C we obtain the ML estimation of the clock skew as written in (11), where is obtained in Appendix B as given in (12).
This completes our proof.

3.2. CRLB for the ML Estimation of Clock Skew

Theorem 2. For the following assumption
(i)   , where is the fixed interval between consecutive PTP Sync messages.
The CRLB of the clock skew ML estimation may be defined aswhere denotes the variance of and is defined byand is defined by

Proof. The inequality of the Cramer-Rao bound is given byFirst, we will calculate the second derivative of the log-likelihood function.
From (36) and Appendix C we can write the derivative of the log-likelihood function with respect to asInserting (16) into (42) givesTherefore, according to Appendix D, the second derivative of the log-likelihood function with respect to may be written aswhere the full proof is provided in Appendix D.
Next, we will derive the mean of (44) and according to Appendix D we may writeFinally, by substituting (45) into (41) and using some properties obtained in Appendix D, we come back to the expression for the CRLB as given in (38).
This completes our proof.

4. Simulation

In this section we show the MSE performance comparison between our proposed expression for the clock skew (11) compared with the conventional dual slaves clock method proposed in [18]. We also test the performance of our new proposed expression for the clock skew (11) for two different values of and . In our simulation the Sync interval and the clock skew have the same values as in [18]: the Sync interval  sec and the clock skew  ppm.

Figures 3 and 4 show the MSE performance comparison between the proposed approach (11) and the conventional dual slaves clock method given in equation in [18] for and , respectively. Namely, the MSE performance comparisons (Figures 3 and 4) are carried out for the case where the PDV is modeled with fGn. The standard deviation is the same value as was taken in [18] for the random jitter ( ps). Figures 3 and 4 show that our proposed expression (11) for the clock skew outperforms (in the MSE point of view) the dual slaves clock method given in equation in [18] for both cases and . Thus we may conclude that the proposed algorithm for the clock skew [18] may not fit for a realistic Ethernet network where PDV exists. But it should be pointed out that the proposed algorithm for the clock skew [18] is very simple compared to our proposed approach (11). As a matter of fact the computational complexity of (11) is while the computational complexity of [18] is only (where is defined as and “” is a constant).

Figures 58 show the MSE performance obtained from our proposed clock skew estimator (11) as a function of for two different values of and .

Since the simulation tool (MATLAB) can not calculate the Gamma function for an argument higher than 175, we can not run the simulation for . Therefore, we have to reduce the value for in order to be able to show results also for higher values of the Hurst exponent for which the MSE approaches the CRLB. Nevertheless, in order to show it, Figures 5 and 6 were produced with a lower value for (than chosen in Figures 7-8) in order to show that the MSE approaches the CRLB also for higher values of the Hurst exponent as expected. According to Figures 5-6 we can see that the MSE declines and approaches the CRLB (38) as increases. However, according to Figures 7 and 8 we can see that even for the higher values for the MSE declines and goes in the right direction towards the CRLB, which means that for greater values for we will see a complete match between the MSE to the CRLB.

5. Conclusion

In this paper we developed a new clock skew method for the case where the PDV is modeled as fGn, as proposed in the literature for a realistic Ethernet network. Simulation results indicate that our new proposed method for the clock skew leads the system with a lower MSE compared with the literature proposed dual slave clocks in a slave. In this paper we derived also the CRLB ( dependent) for our proposed clock skew method and have shown that the simulated MSE obtained from our proposed method achieves asymptotically our new derived CRLB.

Appendices

A. Log-Likelihood Function Derivation

In this appendix we derive in detail the derivative of the log-likelihood function with respect to .

We also show the calculation of the log-likelihood function maximum value.

Equation (36) determinesSince is not dependent on , we can writeThen, we substitute (16) into (A.2) and obtainNow, we will derive the internal expression with respect to and getNext, we set (A.4) to zero in order to find the log-likelihood function maximum value:After some basic algebraic operations we getHence, by substituting (17) into (A.6) we obtainSince we can say thatNow, we will multiply the numerator and denominator in the first term of (A.7) by :Inserting (A.8) into (A.9) givesSince the random jitters and are assumed to be independent we will writeInserting (A.11) into (A.10) givesSince the random jitter is assumed to be an independent Gaussian RV with zero mean and variance of we can writeThereforewhere only when . because occurs times. because or occurs times.

We will use the same technique regarding :Thereforewhere only when . because occurs times. because or occurs times.

Inserting (A.14) and (A.16) into (A.12) givesFinally, according to Appendix C we may write (A.17) with the proper sum limitation:

B. Derivation

In this appendix we derive in detail the expression for .

Earlier we defined in (34):By setting (32) into (B.1) we obtainNow, by substituting (21) into (B.2) and performing some algebraic operation we obtainWe will write with the Gamma function instead of the factorial function because .

The main character of the Gamma function isThereforeFinally, we will arrange the expression for more elegantlyThis is the final form of mentioned in (12).

C. Log-Likelihood Function First Derivative Sum Limitation

In this Appendix we derive in detail the sum limitations: .

We can get the above limitation by two different ways:(i)examining the limitations on the Gamma function arguments that exist in ,(ii)examining the limitation on the expression mentioned in [50] that built up the expression.The first way:

(i) The Gamma function has simple poles at .

Now, under the above knowledge, we examine all the Gamma function arguments existing in (B.6):(i.A) always exists;(i.B)   exists under (1.J);(i.C) always exists;(i.D) exists under (1.K);(i.E) always exists;(i.F) always exists;(i.G) always exists;(i.H) always exists;(i.I) always exists;(i.J);(i.K);(i.L) exists under (1.J);(i.M) exists under (1.K);(i.N) always exists;(i.O) always exists.From (i.J) and (i.K) we get the limitation: .

Note that (i.A)–(i.I) and (i.L)–(i.O) are always true for the abovementioned limitations.

Therefore, the conclusion from the first way is

The second way:

(ii) From the limitation mentioned in [50, page 2981], we can learn that is only defined for and .

Therefore(ii.A)   and ,(ii.B) and ,(ii.C) and ,(ii.D) and .From (ii.A) and (ii.B) and from the fact that we get that .

From (ii.C) and (ii.D) we get the limitation: .

In conclusion we get the same limitations as in (C.1).

D. Log-Likelihood Function-Second Derivative

In this Appendix we derive in detail the Cramer-Rao lower bound (CRLB).

The inequality of the Cramer-Rao bound as given in (41) isIn order to calculate the second derivative of the log-likelihood function we will write the first derivative of the log-likelihood function with respect to .

We can write it from (36)Then, we substitute (16) for the th sample and (16) for the th sample into (D.2) and obtainNext, we will take the derivative according to the product role and getNow, let us calculate the second derivative of the log-likelihood function by deriving (D.4) with respect to Inserting (16) into (D.5) we obtainThereforeNow, we will calculate the mean of the second derivative of the log-likelihood function:Because and are deterministic values we can saySince and are assumed to be with zero mean, we can writeNow, by substituting (D.10) into (D.9) we obtainAccording to Appendices E and F and by substituting (20) into (D.11) we obtainwhere

Then, by substituting (35) and (D.13) into (D.1) we obtain

Let us divide the numerator and denominator by

In order to find we will set into (23) and obtainSimilarly, in order to find we will set into (20) and obtainBy substituting (D.19) into (D.18) we obtainThen, we substitute (D.20) into (D.17) and obtain

Let us divide the numerator and denominator by

Finally, by substituting (D.14) into (D.22) we obtain

According to Appendices E and G and the above equation (D.23) the CRLB is defined bywhereand is obtained by

E. Calculation

In this Appendix we show in detail how to calculate .

According to equation given in [50], the covariance function of a fractionally differenced Gaussian noise process (fdGn) is given byTherefore, we can writeOn the other hand, we defined in this paper (assumption in subsection “ML Estimation of Clock Skew”) , asBy substituting (E.3) into we can say thatNext, by substituting (E.1) into (E.4) we obtainBy comparing (E.2) with (E.5) we getFinally, by substituting (D.14) into (E.6), we get

Therefore, the definition of is given bywhere

F. Calculation

In this Appendix we show in detail how to calculate .

We will calculate according to the Rescaled-Range () analysis.

According to [51, 52] in order to estimate the Hurst exponent of a time series ( denotes vector of data set), we should execute the following process.

First, the time series of full length is divided into a number of shorter time series of length .

The average Rescaled-Range is then calculated for each value of .

The Rescaled-Range is calculated as follows.

(i) Mean calculation:

The mean of the time series is calculated as follows:

(ii) Adjusted mean series calculation:

(iii) Cumulative deviate series calculation:

(iv) Range series calculation:

(v) Standard deviation series calculation:

(vi) Average over all the partial time series of length of the Rescaled-Range :And finally are (vii) plotting the logarithm as a function of and fitting a straight line according to the linear regression method. The slope of the line gives .

Figure 9 shows the plot of the Rescaled-Range analysis as a function of (where is defined as the log function with base of 2) for the case of .

Operating the above analysis on several data sets (corresponding to several values of ) and by taking the mean value gives the following results:The values for are denoted by and the values for are denoted by .

By a second-order polynominal fit preformed in Excel, we find the following strong () connection between and :where is the Pearson product-moment correlation coefficient and the value can be interpreted as the proportion of the variance in attributable to the variance in .

is a measure of the dependence between two variables, giving a value between and −1 inclusive, where 1 is total positive correlation. It is widely used in the sciences as a measure of the degree of dependence between two variables.

can be calculated by the following equation:where(i) and are the mean of the and data sets, respectively.(ii) and are the standard deviation of the and data sets, respectively.

G. Derivation

In this Appendix we derive in detail the definition of and its sum limitations: .

Inserting (E.8) into (D.23) gives

Therefore, is defined bywhere is defined asNext we will explain the sum limitations of .

First, the sum limitations contains at least the same limitations as in (shown in Appendix C) because is a product of with additional term. Therefore, the first limitation will beSecond, we will examine the limitations on the Gamma function arguments existing in the second term of :As already mentioned in this paper, the Gamma function has simple poles at ; therefore(A) always exists;(B) always exists under ;(C) always exists;(D) always exists under ;(E) always exists;(F) always exists;(G) always exists;(H) always exists under ;(I) always exists under ;(J) always exists under ;(K) always exists under ;(L) always exists under ;(M) always exists under ;Finally, we will examine the terms in the (G.5) denominator:(i) always exists; .(ii) always exists; .Therefore, the conclusion is that has the same limitations like .

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors would like to thank the anonymous reviewers for their helpful comments.