# Maximum Likelihood Estimation of Clock Skew in IEEE 1588 with Fractional Gaussian Noise

Authors: Chagai Levy; Monika Pinchas
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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 (H) in the range of 0.5
<
H
<
1. 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 (0.5
<
H
<
1). 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.


---

## Body

## 1. Introduction

Recently, the topic of synchronization in network nodes has attracted significant attention for a number of applications [1–10]. 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 [25–30] as LRD characteristics such as fGn with Hurst exponent in the range of 0.5
<
H
<
1.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 H) is a widely used traditional traffic model and refer the reader to [33–35]. They [32] also add that traditional methods to synthesize traffic are based on fGN with a single parameter and refer the reader to [36–40]. 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 [42–48] that model the PDV as a Gaussian or exponential delay. These papers [42–48] do not need the symmetric assumption on the up- and down link unlike [20, 21, 23]. However, they ([42–48]) 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 [42–48]) 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.Table 1
Comparison between the various estimation methods.
Model
Reference
Up- and Downlink assumption
PDV model existence in the algorithm
PDV model
Exposure to PDV (2WD/OWD)
Basic concept of IEEE 1588
[20, 21, 23]
Symmetry
—
—
Up- and downlink
Dual slave clock in a slave
[18]
Asymmetry/symmetry
—
—
Uplink only
Gaussian delay model
[42, 46]
Asymmetry/symmetry
✓
Gaussian
Up- and downlink
Exponential delay model
[42–45, 47, 48]
Asymmetry/symmetry
✓
Exponential
Up- and downlink
MLE of clock skew in IEEE 1588 with fGn
(New model)
Asymmetry/symmetry
✓
fGn
Uplink onlyIn 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 0.5
<
H
<
1. 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 Section2. 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 T
m
.
1 when the first Sync message is sent. When a slave clock receives the Sync message, it measures and stores the time of reception T
s
.
1. The master then sends a Follow_Up message that contains the actual value of the timestamp T
m
.
1. At this point the slave calculates the master-to-slave delay D
m
2
s [18]:(1)
D
m
2
s
=
T
s
.
1
-
T
m
.
1
.On the other hand, the master-to-slave delay D
m
2
s 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 Message
D
m
2
s [22]:(2)
D
m
2
s
=
θ
+
Message
D
m
2
s
.Figure 1
PTP messaging timing diagram.In order to estimate the clock offsetθ the slave periodically sends the Delay_Req message and stores the transmission time with a timestamp T
s
.
1
′. When the master receives the Delay_Req message, it sends a Delay_Resp message, which conveys the timestamp T
m
.
1
′ of the reception time of the request message. The slave calculates the slave-to-master delay D
s
2
m from these timestamps as [18](3)
D
s
2
m
=
T
m
.
1
′
-
T
s
.
1
′
.On the other hand, the slave-to-master delay D
s
2
m 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 Message
D
s
2
m [22]:(4)
D
s
2
m
=
-
θ
+
Message
D
s
2
m
.Assuming that the master-to-slave message delay Message
D
m
2
s is equal to the slave-to-master message delay Message
D
s
2
m we may write [22](5)
D
w
=
Message
D
m
2
s
=
Message
D
s
2
m
,where D
w is the estimated message delay of the uplink or downlink. By using (1) to (5), the slave estimates D
w and the clock offset (θ) as(6)
D
w
=
D
m
2
s
+
D
s
2
m
2
,
(7)
θ
=
D
m
2
s
-
D
w
=
D
m
2
s
-
D
s
2
m
2
.The 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 frequencyf
0 in addition to a normalized frequency difference of x 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 (f
0
+
x
f
0) and slave's frequency is (f
0
-
x
f
0); 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 T
m
.
i is the time when the ith Sync message is sent. T
s
1
.
i and T
s
2
.
i are the times that slave clocks 1 and 2 measure at the receive of the ith Sync message, respectively. Notably, the dual slave clocks method requires only one-way time transfer. Clock 2 runs at a normalized clock period of (
1
+
ε
) 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 (
1
+
ε
), which is half of the slave clock 1 running at a normalized clock period of 2
(
1
+
ε
). Both slave clocks have the same offset, which is assumed to be θ. Thus, according to [18] we may express the times T
s
1
.
i, T
s
2
.
i as(8)
T
s
1
.
i
=
2
1
+
ε
T
m
.
i
+
D
m
2
s
+
θ
+
ϕ
1
,
i
,
T
s
2
.
i
=
1
+
ε
T
m
.
i
+
D
m
2
s
+
θ
+
ϕ
2
,
i
,where(i)
ϕ
1
,
i and ϕ
2
,
i are random jitters of slave clock 1 and 2 within the ith period, respectively.
(ii)
The random jitterϕ
1
,
i is assumed to be an independent Gaussian random variable (RV) with zero mean and variance σ
2.
(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ϕ
2
,
i is 2
σ
2 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].Figure 2
Synchronization of dual slave clocks.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 D
m
2
s we saw in [18], and X
i will be denoted by the PDV of the ith cycle. Therefore, (8) must be written in the following form:(9)
T
s
1
.
i
=
2
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
+
θ
+
ϕ
1
,
i
,
(10)
T
s
2
.
i
=
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
+
θ
+
ϕ
2
,
i
,where X
i is assumed to be fractional Gaussian noise (fGn) with zero mean with Hurst exponent H
x and variance σ
W
x
2.


## 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 T
m
.
1 when the first Sync message is sent. When a slave clock receives the Sync message, it measures and stores the time of reception T
s
.
1. The master then sends a Follow_Up message that contains the actual value of the timestamp T
m
.
1. At this point the slave calculates the master-to-slave delay D
m
2
s [18]:(1)
D
m
2
s
=
T
s
.
1
-
T
m
.
1
.On the other hand, the master-to-slave delay D
m
2
s 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 Message
D
m
2
s [22]:(2)
D
m
2
s
=
θ
+
Message
D
m
2
s
.Figure 1
PTP messaging timing diagram.In order to estimate the clock offsetθ the slave periodically sends the Delay_Req message and stores the transmission time with a timestamp T
s
.
1
′. When the master receives the Delay_Req message, it sends a Delay_Resp message, which conveys the timestamp T
m
.
1
′ of the reception time of the request message. The slave calculates the slave-to-master delay D
s
2
m from these timestamps as [18](3)
D
s
2
m
=
T
m
.
1
′
-
T
s
.
1
′
.On the other hand, the slave-to-master delay D
s
2
m 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 Message
D
s
2
m [22]:(4)
D
s
2
m
=
-
θ
+
Message
D
s
2
m
.Assuming that the master-to-slave message delay Message
D
m
2
s is equal to the slave-to-master message delay Message
D
s
2
m we may write [22](5)
D
w
=
Message
D
m
2
s
=
Message
D
s
2
m
,where D
w is the estimated message delay of the uplink or downlink. By using (1) to (5), the slave estimates D
w and the clock offset (θ) as(6)
D
w
=
D
m
2
s
+
D
s
2
m
2
,
(7)
θ
=
D
m
2
s
-
D
w
=
D
m
2
s
-
D
s
2
m
2
.The 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 frequencyf
0 in addition to a normalized frequency difference of x 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 (f
0
+
x
f
0) and slave's frequency is (f
0
-
x
f
0); 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 T
m
.
i is the time when the ith Sync message is sent. T
s
1
.
i and T
s
2
.
i are the times that slave clocks 1 and 2 measure at the receive of the ith Sync message, respectively. Notably, the dual slave clocks method requires only one-way time transfer. Clock 2 runs at a normalized clock period of (
1
+
ε
) 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 (
1
+
ε
), which is half of the slave clock 1 running at a normalized clock period of 2
(
1
+
ε
). Both slave clocks have the same offset, which is assumed to be θ. Thus, according to [18] we may express the times T
s
1
.
i, T
s
2
.
i as(8)
T
s
1
.
i
=
2
1
+
ε
T
m
.
i
+
D
m
2
s
+
θ
+
ϕ
1
,
i
,
T
s
2
.
i
=
1
+
ε
T
m
.
i
+
D
m
2
s
+
θ
+
ϕ
2
,
i
,where(i)
ϕ
1
,
i and ϕ
2
,
i are random jitters of slave clock 1 and 2 within the ith period, respectively.
(ii)
The random jitterϕ
1
,
i is assumed to be an independent Gaussian random variable (RV) with zero mean and variance σ
2.
(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ϕ
2
,
i is 2
σ
2 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].Figure 2
Synchronization of dual slave clocks.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 D
m
2
s we saw in [18], and X
i will be denoted by the PDV of the ith cycle. Therefore, (8) must be written in the following form:(9)
T
s
1
.
i
=
2
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
+
θ
+
ϕ
1
,
i
,
(10)
T
s
2
.
i
=
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
+
θ
+
ϕ
2
,
i
,where X
i is assumed to be fractional Gaussian noise (fGn) with zero mean with Hurst exponent H
x and variance σ
W
x
2.


## 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 Section3.1 we derive the ML estimation for the clock skew under the assumption of fGn PDV with Hurst exponent in the range of 0.5
<
H
x
<
1.In Section3.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 of0.5
<
H
x
<
1,
(ii)
ξ
i
=
X
i
-
X
i
-
1,
ξ
i is assumed to be fGn with zero mean with Hurst exponent H
ξ and variance σ
W
ξ
2.The ML estimation for the clock skew may be defined as(11)
ε
^
=
2
∑
i
=
1
N
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
(
i
,
j
)
-
1
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
A
N
,
i
,
j
,
k
,
H
ξ
∑
i
=
1
N
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
‍
·
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
A
(
N
,
i
,
j
,
k
,
H
ξ
)
·
Δ
T
m
.
j
(
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
)
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
-
1
-
1
,where A
(
N
,
i
,
j
,
k
,
H
ξ
) is defined by(12)
A
N
,
i
,
j
,
k
,
H
ξ
=
Γ
N
Γ
-
H
ξ
+
0.5
Γ
N
-
H
ξ
+
0.5
2
·
Γ
i
-
k
-
H
ξ
+
0.5
Γ
N
-
H
ξ
+
0.5
-
i
+
k
·
Γ
(
j
-
k
-
H
ξ
+
0.5
)
Γ
(
N
-
H
ξ
+
0.5
-
j
+
k
)
·
Γ
N
-
i
+
k
Γ
N
-
j
+
k
Γ
i
-
k
Γ
j
-
k
-
1
·
1
i
-
k
1
j
-
k
-
1
N
-
i
+
k
1
N
-
j
+
k
.
Γ
(
·
) denotes Gamma function and N 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) yields(13)
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
=
T
s
1
.
i
-
T
s
2
.
i
+
ϕ
2
,
i
-
ϕ
1
,
i
.Subtracting (13) of the (
i
-
1
)th Sync message from that of the ith Sync message indicates that(14)
(
1
+
ε
)
(
T
m
.
i
-
T
m
.
(
i
-
1
)
+
X
i
-
X
i
-
1
)
=
(
T
s
1
.
i
-
T
s
1
.
(
i
-
1
)
)
-
(
T
s
2
.
i
-
T
s
2
.
(
i
-
1
)
)
+
(
ϕ
2
,
i
-
ϕ
2
,
i
-
1
)
-
(
ϕ
1
,
i
-
ϕ
1
,
i
-
1
)which can be also written as(15)
1
+
ε
Δ
T
m
.
i
+
ξ
i
=
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
.Thus, ξ
i may be expressed as(16)
ξ
i
=
ε
′
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
-
Δ
T
m
.
i
,where according to [42](17)
ε
′
≜
1
1
+
ε
;
ε
≠
0
.According to [49, 50], the probability density function (PDF) of vector ξ
_ (of length N) may be given by(18)
p
ξ
_
;
H
ξ
,
σ
W
ξ
2
=
1
2
π
N
/
2
R
1
/
2
exp
⁡
-
1
2
ξ
_
T
R
-
1
ξ
_
,where ξ
_
T denotes the transpose of ξ
_ and R denotes the determinant of the matrix R.
R is the covariance matrix of the vector ξ
_:(19)
R
H
ξ
,
σ
W
ξ
2
=
E
ξ
_
ξ
_
T
.The elements of R are given by [50, 51](20)
R
i
j
=
r
i
-
j
=
r
k
=
σ
W
ξ
2
(
-
2
d
ξ
)
!
-
d
ξ
!
2
Γ
(
1
-
d
ξ
)
Γ
(
k
+
d
ξ
)
Γ
(
d
ξ
)
Γ
(
k
+
1
-
d
ξ
)
=
σ
W
ξ
2
Γ
-
2
H
ξ
+
2
Γ
k
+
H
ξ
-
0.5
Γ
-
H
ξ
+
1.5
Γ
H
ξ
-
0.5
Γ
k
+
1.5
-
H
ξ
,where, according to [51, 52],(21)
d
ξ
=
H
ξ
-
0.5
,
d
ξ
∈
(
-
0.5,0.5
)
;
H
ξ
∈
(
0,1
)
.Note that H
x
≠
0.5 so H
ξ
≠
0.5.
R can be decomposed as was done in [49, 50](22)
R
=
σ
W
ξ
2
R
1
=
σ
W
ξ
2
r
1
0
r
1
1
⋯
r
1
N
-
1
r
1
1
r
1
0
r
1
1
·
r
1
2
·
⋯
·
·
·
⋯
·
r
1
N
-
1
r
1
N
-
2
⋯
r
1
0
,where(23)
r
1
i
-
j
=
r
[
i
-
j
]
σ
W
ξ
2
.According to [50] we can write(24)
R
2
=
r
2
[
0
]
r
2
[
1
]
⋯
r
2
[
N
-
1
]
r
2
[
1
]
r
2
[
0
]
r
2
[
1
]
·
r
2
[
2
]
·
⋯
·
·
·
⋯
·
r
2
[
N
-
1
]
r
2
[
N
-
2
]
⋯
r
2
[
0
]
,where(25)
r
2
i
-
j
=
r
1
i
-
j
r
1
0
.Next, we recall from [50] the expression for the log-likelihood function:(26)
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
=
log
⁡
p
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
=
-
N
2
log
⁡
⁡
(
2
π
)
-
N
2
log
⁡
⁡
(
σ
W
ξ
2
)
-
1
2
log
⁡
R
1
-
1
2
σ
W
ξ
2
ξ
_
T
R
1
-
1
ξ
_
.The derivative of the log-likelihood function (26) with respect to ε
′ yields(27)
∂
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
∂
∂
ε
′
-
1
2
σ
W
ξ
2
ξ
_
T
R
1
-
1
ξ
_
.Note that [50] has taken the derivative of (26) with respect to d. Thus, our final expression is very different.
From (24) and (25) we can write(28)
R
1
=
r
1
0
R
2
.By substituting (28) into (27) we obtain(29)
∂
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
-
1
2
σ
W
ξ
2
r
1
0
∂
∂
ε
′
ξ
_
T
R
2
-
1
ξ
_
.According to [50], the expression ξ
_
T
R
2
-
1
ξ
_ can be written as(30)
ξ
_
T
R
2
-
1
ξ
_
=
1
P
N
-
1
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
i
-
k
a
N
-
1
j
-
k
-
a
N
-
1
N
-
i
+
k
a
N
-
1
N
-
j
+
k
ξ
i
ξ
j
,where(31)
P
N
-
1
=
P
N
-
2
1
-
H
ξ
-
0.5
2
N
-
0.5
-
H
ξ
2
;
P
0
=
1
,
(32)
a
k
i
=
k
i
i
-
d
ξ
-
1
!
k
-
d
ξ
-
i
!
-
d
ξ
-
1
!
k
-
d
ξ
!
for
1
≤
i
≤
k
.Then, we substitute (30) into (29) and obtain(33)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
-
1
2
σ
W
ξ
2
r
1
[
0
]
1
P
N
-
1
×
∂
∂
ε
′
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
(
i
-
k
)
a
N
-
1
(
j
-
k
)
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
-
a
N
-
1
(
N
-
i
+
k
)
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
N
-
j
+
k
ξ
i
ξ
j
.Let us define for simplicity the next two variables:(34)
A
N
,
i
,
j
,
k
,
H
ξ
≜
a
N
-
1
(
i
-
k
)
a
N
-
1
(
j
-
k
)
-
a
N
-
1
N
-
i
+
k
a
N
-
1
N
-
j
+
k
,
(35)
K
≜
-
1
2
σ
W
ξ
2
r
1
0
1
P
N
-
1
.Inserting (34) and (35) into (33) gives(36)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
K
∂
∂
ε
′
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
A
N
,
i
,
j
,
k
,
H
ξ
ξ
i
ξ
j
.Next, the derivative of the log-likelihood function (36) with respect to ε
′ is carried out and set to zero(37)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
0
.Appendix 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 AppendicesA and C we obtain the ML estimation of the clock skew ε
^ as written in (11), where A
(
N
,
i
,
j
,
k
,
H
ξ
) 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)Δ
T
m
.
i
=
Δ
T
m
.
j
=
β, where β is the fixed interval between consecutive PTP Sync messages.
The CRLB of the clock skew ML estimation may be defined as(38)
Var
ε
^
-
ε
≥
1
+
ε
2
P
N
-
1
∑
i
=
1
N
∑
j
=
1
N
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
B
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
,where Var
[
·
] denotes the variance of [
·
] and B
(
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
) is defined by(39)
B
(
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
)
=
A
N
,
i
,
j
,
k
,
H
ξ
Γ
1.5
-
H
ξ
Γ
j
-
i
+
H
ξ
-
0.5
Γ
H
ξ
-
0.5
Γ
j
-
i
+
1.5
-
H
ξ
·
3
+
β
2
C
i
,
j
,
H
x
,
σ
W
x
2and C
(
i
,
j
,
H
x
,
σ
W
x
2
) is defined by(40)
C
(
i
,
j
,
H
x
,
σ
W
x
2
)
=
σ
W
x
2
Γ
(
-
2
H
x
+
2
)
Γ
(
-
H
x
+
1.5
)
1
Γ
(
H
x
-
0.5
)
·
2
Γ
j
-
i
+
H
x
-
0.5
Γ
j
-
i
+
1.5
-
H
x
-
Γ
j
-
i
-
1.5
+
H
x
Γ
j
-
i
-
H
x
+
0.5
-
Γ
i
-
j
-
1.5
+
H
x
Γ
i
-
j
-
H
x
+
0.5
.Proof.
The inequality of the Cramer-Rao bound is given by(41)
Var
ε
^
-
ε
≥
-
E
∂
2
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
2
-
1
.First, 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 ε as(42)
∂
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
=
K
∂
∂
ε
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
ξ
i
ξ
j
.Inserting (16) into (42) gives(43)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
=
K
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
∑
d
2
d
2
·
∂
∂
ε
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
1
+
ε
∑
d
2
d
2
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
1
+
ε
-
Δ
T
m
.
i
·
(
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
+
Δ
ϕ
2
,
j
-
Δ
ϕ
1
,
j
)
(
1
+
ε
)


∑
d
2
d
2
(
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
+
Δ
ϕ
2
,
j
-
Δ
ϕ
1
,
j
)
(
1
+
ε
)
-
Δ
T
m
.
j
.Therefore, according to Appendix D, the second derivative of the log-likelihood function with respect to ε may be written as(44)
∂
2
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
2
=
2
1
+
ε
2
K
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
β
+
ξ
i
+
β
+
ξ
j
ξ
i
·
β
+
ξ
i
ξ
j
+
β
+
ξ
j
β
+
ξ
i
+
β
+
ξ
j
ξ
i
·
β
+
ξ
i
ξ
j
+
β
+
ξ
j
β
+
ξ
i
+
β
+
ξ
j
ξ
i
,where the full proof is provided in Appendix D.
Next, we will derive the mean of (44) and according to Appendix D we may write(45)
E
∂
2
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
2
=
2
1
+
ε
2
·
K
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
·
E
3
ξ
i
ξ
j
+
2
β
ξ
j
+
2
β
ξ
i
+
β
2
.Finally, 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.


## 3.1. ML Estimation of Clock Skew

Theorem 1.
For the following assumptions(i)
the Hurst exponent is in the range of0.5
<
H
x
<
1,
(ii)
ξ
i
=
X
i
-
X
i
-
1,
ξ
i is assumed to be fGn with zero mean with Hurst exponent H
ξ and variance σ
W
ξ
2.The ML estimation for the clock skew may be defined as(11)
ε
^
=
2
∑
i
=
1
N
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
(
i
,
j
)
-
1
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
A
N
,
i
,
j
,
k
,
H
ξ
∑
i
=
1
N
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
‍
·
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
A
(
N
,
i
,
j
,
k
,
H
ξ
)
·
Δ
T
m
.
j
(
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
)
d
d
d
d
d
d
d
d
d
d
d
d
d
d
d
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
+
Δ
T
m
.
i
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
-
1
-
1
,where A
(
N
,
i
,
j
,
k
,
H
ξ
) is defined by(12)
A
N
,
i
,
j
,
k
,
H
ξ
=
Γ
N
Γ
-
H
ξ
+
0.5
Γ
N
-
H
ξ
+
0.5
2
·
Γ
i
-
k
-
H
ξ
+
0.5
Γ
N
-
H
ξ
+
0.5
-
i
+
k
·
Γ
(
j
-
k
-
H
ξ
+
0.5
)
Γ
(
N
-
H
ξ
+
0.5
-
j
+
k
)
·
Γ
N
-
i
+
k
Γ
N
-
j
+
k
Γ
i
-
k
Γ
j
-
k
-
1
·
1
i
-
k
1
j
-
k
-
1
N
-
i
+
k
1
N
-
j
+
k
.
Γ
(
·
) denotes Gamma function and N 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) yields(13)
1
+
ε
T
m
.
i
+
D
m
2
s
+
X
i
=
T
s
1
.
i
-
T
s
2
.
i
+
ϕ
2
,
i
-
ϕ
1
,
i
.Subtracting (13) of the (
i
-
1
)th Sync message from that of the ith Sync message indicates that(14)
(
1
+
ε
)
(
T
m
.
i
-
T
m
.
(
i
-
1
)
+
X
i
-
X
i
-
1
)
=
(
T
s
1
.
i
-
T
s
1
.
(
i
-
1
)
)
-
(
T
s
2
.
i
-
T
s
2
.
(
i
-
1
)
)
+
(
ϕ
2
,
i
-
ϕ
2
,
i
-
1
)
-
(
ϕ
1
,
i
-
ϕ
1
,
i
-
1
)which can be also written as(15)
1
+
ε
Δ
T
m
.
i
+
ξ
i
=
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
.Thus, ξ
i may be expressed as(16)
ξ
i
=
ε
′
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
-
Δ
T
m
.
i
,where according to [42](17)
ε
′
≜
1
1
+
ε
;
ε
≠
0
.According to [49, 50], the probability density function (PDF) of vector ξ
_ (of length N) may be given by(18)
p
ξ
_
;
H
ξ
,
σ
W
ξ
2
=
1
2
π
N
/
2
R
1
/
2
exp
⁡
-
1
2
ξ
_
T
R
-
1
ξ
_
,where ξ
_
T denotes the transpose of ξ
_ and R denotes the determinant of the matrix R.
R is the covariance matrix of the vector ξ
_:(19)
R
H
ξ
,
σ
W
ξ
2
=
E
ξ
_
ξ
_
T
.The elements of R are given by [50, 51](20)
R
i
j
=
r
i
-
j
=
r
k
=
σ
W
ξ
2
(
-
2
d
ξ
)
!
-
d
ξ
!
2
Γ
(
1
-
d
ξ
)
Γ
(
k
+
d
ξ
)
Γ
(
d
ξ
)
Γ
(
k
+
1
-
d
ξ
)
=
σ
W
ξ
2
Γ
-
2
H
ξ
+
2
Γ
k
+
H
ξ
-
0.5
Γ
-
H
ξ
+
1.5
Γ
H
ξ
-
0.5
Γ
k
+
1.5
-
H
ξ
,where, according to [51, 52],(21)
d
ξ
=
H
ξ
-
0.5
,
d
ξ
∈
(
-
0.5,0.5
)
;
H
ξ
∈
(
0,1
)
.Note that H
x
≠
0.5 so H
ξ
≠
0.5.
R can be decomposed as was done in [49, 50](22)
R
=
σ
W
ξ
2
R
1
=
σ
W
ξ
2
r
1
0
r
1
1
⋯
r
1
N
-
1
r
1
1
r
1
0
r
1
1
·
r
1
2
·
⋯
·
·
·
⋯
·
r
1
N
-
1
r
1
N
-
2
⋯
r
1
0
,where(23)
r
1
i
-
j
=
r
[
i
-
j
]
σ
W
ξ
2
.According to [50] we can write(24)
R
2
=
r
2
[
0
]
r
2
[
1
]
⋯
r
2
[
N
-
1
]
r
2
[
1
]
r
2
[
0
]
r
2
[
1
]
·
r
2
[
2
]
·
⋯
·
·
·
⋯
·
r
2
[
N
-
1
]
r
2
[
N
-
2
]
⋯
r
2
[
0
]
,where(25)
r
2
i
-
j
=
r
1
i
-
j
r
1
0
.Next, we recall from [50] the expression for the log-likelihood function:(26)
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
=
log
⁡
p
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
=
-
N
2
log
⁡
⁡
(
2
π
)
-
N
2
log
⁡
⁡
(
σ
W
ξ
2
)
-
1
2
log
⁡
R
1
-
1
2
σ
W
ξ
2
ξ
_
T
R
1
-
1
ξ
_
.The derivative of the log-likelihood function (26) with respect to ε
′ yields(27)
∂
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
∂
∂
ε
′
-
1
2
σ
W
ξ
2
ξ
_
T
R
1
-
1
ξ
_
.Note that [50] has taken the derivative of (26) with respect to d. Thus, our final expression is very different.
From (24) and (25) we can write(28)
R
1
=
r
1
0
R
2
.By substituting (28) into (27) we obtain(29)
∂
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
-
1
2
σ
W
ξ
2
r
1
0
∂
∂
ε
′
ξ
_
T
R
2
-
1
ξ
_
.According to [50], the expression ξ
_
T
R
2
-
1
ξ
_ can be written as(30)
ξ
_
T
R
2
-
1
ξ
_
=
1
P
N
-
1
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
i
-
k
a
N
-
1
j
-
k
-
a
N
-
1
N
-
i
+
k
a
N
-
1
N
-
j
+
k
ξ
i
ξ
j
,where(31)
P
N
-
1
=
P
N
-
2
1
-
H
ξ
-
0.5
2
N
-
0.5
-
H
ξ
2
;
P
0
=
1
,
(32)
a
k
i
=
k
i
i
-
d
ξ
-
1
!
k
-
d
ξ
-
i
!
-
d
ξ
-
1
!
k
-
d
ξ
!
for
1
≤
i
≤
k
.Then, we substitute (30) into (29) and obtain(33)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
-
1
2
σ
W
ξ
2
r
1
[
0
]
1
P
N
-
1
×
∂
∂
ε
′
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
(
i
-
k
)
a
N
-
1
(
j
-
k
)
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
-
a
N
-
1
(
N
-
i
+
k
)
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
a
N
-
1
N
-
j
+
k
ξ
i
ξ
j
.Let us define for simplicity the next two variables:(34)
A
N
,
i
,
j
,
k
,
H
ξ
≜
a
N
-
1
(
i
-
k
)
a
N
-
1
(
j
-
k
)
-
a
N
-
1
N
-
i
+
k
a
N
-
1
N
-
j
+
k
,
(35)
K
≜
-
1
2
σ
W
ξ
2
r
1
0
1
P
N
-
1
.Inserting (34) and (35) into (33) gives(36)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
K
∂
∂
ε
′
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
N
‍
A
N
,
i
,
j
,
k
,
H
ξ
ξ
i
ξ
j
.Next, the derivative of the log-likelihood function (36) with respect to ε
′ is carried out and set to zero(37)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
′
=
0
.Appendix 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 AppendicesA and C we obtain the ML estimation of the clock skew ε
^ as written in (11), where A
(
N
,
i
,
j
,
k
,
H
ξ
) 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)Δ
T
m
.
i
=
Δ
T
m
.
j
=
β, where β is the fixed interval between consecutive PTP Sync messages.
The CRLB of the clock skew ML estimation may be defined as(38)
Var
ε
^
-
ε
≥
1
+
ε
2
P
N
-
1
∑
i
=
1
N
∑
j
=
1
N
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
B
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
,where Var
[
·
] denotes the variance of [
·
] and B
(
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
) is defined by(39)
B
(
N
,
i
,
j
,
k
,
H
,
β
,
σ
W
2
)
=
A
N
,
i
,
j
,
k
,
H
ξ
Γ
1.5
-
H
ξ
Γ
j
-
i
+
H
ξ
-
0.5
Γ
H
ξ
-
0.5
Γ
j
-
i
+
1.5
-
H
ξ
·
3
+
β
2
C
i
,
j
,
H
x
,
σ
W
x
2and C
(
i
,
j
,
H
x
,
σ
W
x
2
) is defined by(40)
C
(
i
,
j
,
H
x
,
σ
W
x
2
)
=
σ
W
x
2
Γ
(
-
2
H
x
+
2
)
Γ
(
-
H
x
+
1.5
)
1
Γ
(
H
x
-
0.5
)
·
2
Γ
j
-
i
+
H
x
-
0.5
Γ
j
-
i
+
1.5
-
H
x
-
Γ
j
-
i
-
1.5
+
H
x
Γ
j
-
i
-
H
x
+
0.5
-
Γ
i
-
j
-
1.5
+
H
x
Γ
i
-
j
-
H
x
+
0.5
.Proof.
The inequality of the Cramer-Rao bound is given by(41)
Var
ε
^
-
ε
≥
-
E
∂
2
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
2
-
1
.First, 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 ε as(42)
∂
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
=
K
∂
∂
ε
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
i
,
j
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
ξ
i
ξ
j
.Inserting (16) into (42) gives(43)
∂
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
=
K
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
∑
d
2
d
2
·
∂
∂
ε
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
1
+
ε
∑
d
2
d
2
Δ
T
s
1
.
i
-
Δ
T
s
2
.
i
+
Δ
ϕ
2
,
i
-
Δ
ϕ
1
,
i
1
+
ε
-
Δ
T
m
.
i
·
(
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
+
Δ
ϕ
2
,
j
-
Δ
ϕ
1
,
j
)
(
1
+
ε
)


∑
d
2
d
2
(
Δ
T
s
1
.
j
-
Δ
T
s
2
.
j
+
Δ
ϕ
2
,
j
-
Δ
ϕ
1
,
j
)
(
1
+
ε
)
-
Δ
T
m
.
j
.Therefore, according to Appendix D, the second derivative of the log-likelihood function with respect to ε may be written as(44)
∂
2
log
⁡
L
(
ξ
_
;
H
ξ
,
σ
W
ξ
2
)
∂
ε
2
=
2
1
+
ε
2
K
·
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
β
+
ξ
i
+
β
+
ξ
j
ξ
i
·
β
+
ξ
i
ξ
j
+
β
+
ξ
j
β
+
ξ
i
+
β
+
ξ
j
ξ
i
·
β
+
ξ
i
ξ
j
+
β
+
ξ
j
β
+
ξ
i
+
β
+
ξ
j
ξ
i
,where the full proof is provided in Appendix D.
Next, we will derive the mean of (44) and according to Appendix D we may write(45)
E
∂
2
log
⁡
L
ξ
_
;
H
ξ
,
σ
W
ξ
2
∂
ε
2
=
2
1
+
ε
2
·
K
∑
i
=
1
N
‍
∑
j
=
1
N
‍
∑
k
=
1
min
⁡
⁡
(
i
,
j
)
-
1
‍
A
N
,
i
,
j
,
k
,
H
ξ
·
E
3
ξ
i
ξ
j
+
2
β
ξ
j
+
2
β
ξ
i
+
β
2
.Finally, 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 σ
W
x
2 and H
x. In our simulation the Sync interval β and the clock skew ε have the same values as in [18]: the Sync interval β
=
2 sec and the clock skew ε
=
40 ppm.Figures3 and 4 show the MSE performance comparison between the proposed approach (11) and the conventional dual slaves clock method given in equation (
11
) in [18] for H
x
=
0.55 and H
x
=
0.8, 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 σ
W
x is the same value as was taken in [18] for the random jitter (±
40 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 (
11
) in [18] for both cases H
x
=
0.55 and H
x
=
0.8. 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 O
(
N
3
) while the computational complexity of [18] is only O
(
N
) (where O
(
Q
~
) is defined as lim
Q
~
→
0
(
O
Q
~
/
Q
~
)
=
r
const and “r
const” is a constant).Figure 3
MSE comparison between the proposed ML estimator for the clock skew (11) to the conventional dual slave method ((
11
) in [18]) as a function of the observation number N. The standard deviation σ
W
x is ±
40 ps.Figure 4
MSE comparison between the proposed ML estimator for the clock skew (11) to the conventional dual slave method ((
11
) in [18]) as a function of the observation number N. The standard deviation σ
W
x is ±
40 ps.Figures5–8 show the MSE performance obtained from our proposed clock skew estimator (11) as a function of N for two different values of σ
W
x
2 and H
x.Figure 5
MSE and CRLB of the clock skew as a function of the observation numberN. The standard deviation σ
W
x of the PDV is ±
0.04 ps.Figure 6
MSE and CRLB of the clock skew as a function of the observation numberN. The standard deviation σ
W
x of the PDV is ±
0.04 ps.Figure 7
MSE and CRLB of the clock skew as a function of the observation numberN. The standard deviation σ
W
x of the PDV is ±
40 ps.Figure 8
MSE and CRLB of the clock skew as a function of the observation numberN. The standard deviation σ
W
x of the PDV is ±
40 ps.Since the simulation tool (MATLAB) can not calculate the Gamma function for an argument higher than 175, we can not run the simulation forN
>
60. Therefore, we have to reduce the value for σ
W
x
2 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 σ
W
x
2 (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 N increases. However, according to Figures 7 and 8 we can see that even for the higher values for σ
W
x
2 the MSE declines and goes in the right direction towards the CRLB, which means that for greater values for N 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 (H 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.



---

*Source: 174289-2015-04-01.xml*