A Low-Computational Burden Closed-Form Approximated Expression for MSE Applicable for PTP with gfGn Environment

Abstract

The Precision Time Protocol (PTP) plays a pivotal role in achieving precise frequency and time synchronization in computer networks. However, network delays and jitter in real systems introduce uncertainties that can compromise synchronization accuracy. Three clock skew estimators designed for the PTP scenario were obtained in our earlier work, complemented by closed-form approximations for the Mean Squared Error (MSE) under the generalized fractional Gaussian noise (gfGn) model, incorporating the Hurst exponent parameter (H) and the a parameter. These expressions offer crucial insights for network designers, aiding in the strategic selection and implementation of clock skew estimators. However, substantial computational resources are required to fit each expression to the gfGn model parameters (H and a) from the MSE perspective requirement. This paper introduces new closed-form estimates that approximate the MSE tailored to match gfGn scenarios that have a lower computational burden compared to the literature-known expressions and that are easily adaptable from the computational burden point of view to different pairs of H and a parameters. Thus, the system requires less substantial computational resources and might be more cost-effective.

1. Introduction

In modern computer networked systems, achieving precise offset (time) and clock skew (frequency) synchronization is essential for ensuring optimal performances and efficiency [1,2,3,4]. The IEEE’s PTP (IEEE1588v2) [5] standard is built on a hardware time-stamped Master–Slave clock packet exchange design. It is a two-way delay (TWD) system, requiring time stamps from both Forward (the path from the Master to the Slave) and Reverse (the path from the Slave to the Master) paths for time and frequency synchronization [5]. The PTP is a popular solution in this domain, facilitating the synchronization of clocks across distributed systems [6,7]. However, PTP systems have challenges in the synchronization task, particularly in environments characterized by complex network dynamics and varying levels of noise and jitter [3]. The packets undergo fixed and random delays that impact synchronization accuracy [1,2,3,4,8]. The fixed delay is determined by the number of components the packet passes until reaching the destination clock [3]. Please note that the Master-to-Slave (Forward) and Slave-to-Master (Reverse) paths might not have the same constant delay [9] since, in each path, the number of components may be different. Packet delay variation (PDV), which is the term for the random delay, is influenced by network traffic [10], as packets may be subject to queuing delays while awaiting transmission through an unblocked output port [3]. Based on [11,12], the PDV significantly impacts synchronization accuracy. Namely, as PDV increases, synchronization accuracy decreases.

Based on the latest research [13,14,15,16,17], gfGn offers a sophisticated framework for modeling complex time series with long-range dependence (LRD) on network jitter and delay fluctuations. The gfGn model in this work is characterized by H and a, where $0.5\le H<1$ and $0<a<1$, respectively. Thus, the gfGn model addresses the challenges of the PDV in real PTP systems.

In our previous work [18,19], we derived three clock skew estimators: (1) the TWD clock skew estimator, which is dependent on both the Master-to-Slave and Slave-to-Master pathways; (2) the one-way delay (OWD) clock skew estimator that depends on the Master-to-Slave path; and (3) the OWD clock skew estimator that depends on the Slave-to-Master path. These estimators also apply to asymmetric paths (where the Forward and Reverse fixed delays are different) and are applicable to the gfGn environment. In addition, we presented in [18,19] three closed-form estimates that approximate the MSE for the gfGn environment related to those clock skew estimators.

The MSE expressions are significant tools for network designers, enabling them to determine the number of Sync messages required to meet system requirements from the MSE perspective. However, the closed-form estimates that approximate the MSE from [18,19] require substantial computational resources to fit the MSE expressions to specific pairs of the Hurst exponent and a parameters.

In this paper, we present new derived low-computational burden, closed-form approximated expressions for the MSE for the gfGn environment for $a\le 0.1$, related to the three clock skew estimators from [18,19]. The newly derived MSE expressions are based on a single expression that can be adapted to each Hurst exponent and a parameter through simple multiplications that require relatively low computational resources. Moreover, simulation results revealed that the low-computational burden MSE expressions are also suitable for the case when $a\le 0.6$.

The newly derived MSE expressions offer several advantages over our previous work: (1) Reduced Computational Burden: the new expressions require significantly less computational power, making them more efficient to implement in software as well as in hardware; (2) Efficiency: by reducing the computational load, the new expressions save time in running simulations; and (3) Cost-Effectiveness: The reduced need for computational resources translates into cost savings in hardware and software, making the implementation of PTP systems more efficient and economical. These benefits underscore the practicality and efficiency of the new low-computational burden, closed-form estimates that approximate the MSE, providing significant tools for network designers to optimize clock skew estimation in PTP systems.

This study’s contributions are as follows:

• The derivation of three low-computational burdens, closed-form approximated expressions for the MSE under the gfGn environment for $a\le 0.1$, related to the three clock skew estimators from our previous work [18,19]. The derived MSE expressions are more computationally efficient and easily adaptable to various combinations of H and a parameters. Unlike before, the Hurst exponent and a parameters are no longer in a nested loop in the MSE calculations. Thus, the MSE result can be obtained faster when changing the Hurst exponent and a parameters from one set to another. The newly derived expression for the MSE has an approximate execution time of a few microseconds in Matlab code when changing the Hurst exponent and a parameters from one set to another. Comparatively, when varying the Hurst exponent and a parameters from one set to another, the Matlab code execution time for the previously obtained expression for the MSE is around minutes or more.
• The derivation of three closed-form estimates that approximate the MSE for the scenario where the a parameter tends toward zero, related to the three clock skew estimators from our previous work [18,19]. These MSE expressions demonstrate improved performances compared to the performances of the MSE expressions for the Gaussian case. These expressions are given as a function of the MSE expressions for the Gaussian case, where it can be easily seen that the MSE expressions for the Gaussian case are multiplied by a factor lower than one.
• Providing a practical tool for system designers to estimate the clock skew estimator’s performance under the gfGn case for $a\le 0.6$, with reduced computational burden, offering a cost-effective alternative to the closed-form estimates that approximate the MSE from our previous work [18,19].

2. System Description

The PTP is a synchronization protocol that relies on a system of clocks (Master and Slave) exchanging packets with time stamps (for a PTP diagram, please refer to Figure 1, recalled from [18]).

These time stamps delineate the temporal relationship between the Master and Slave clocks [1,2,3]. Based on [1,2,3], we can write

$$\begin{array}{c}\hfill t_1\left[j\right]+d_{ms}+{\omega }_1\left[j\right]=t_2\left[j\right](1+\alpha )+Q;t_4\left[j\right]-d_{sm}-{\omega }_2\left[j\right]=t_3\left[j\right](1+\alpha )+Q\end{array}$$

(1)

The PTP time stamps are denoted as $t_1\left[j\right]$, $t_2\left[j\right]$, $t_3\left[j\right]$, and $t_4\left[j\right]$, where j represents the j-th Sync period for $j=1,2,3,\dots ,J$, with J representing the total number of Sync periods. $T_{\mathrm{syn}}$ refers to the Sync message period. Based on Figure 1, the Master initiates the Sync period by sending a Sync message to the Slave at time stamp $t_1\left[j\right]$. The Slave receives the Sync message at time stamp $t_2\left[j\right]$ and sends a Delay request message back to the Master at time stamp $t_3\left[j\right]$, which is received by the Master at time stamp $t_4\left[j\right]$. The fixed delays in the Forward and Reverse paths are respectively designated as $d_{\mathrm{ms}}$ and $d_{\mathrm{sm}}$, while the PDV in the Master-to-Slave and Slave-to-Master paths are named as ${\omega }_1\left[j\right]$ and ${\omega }_2\left[j\right]$, respectively. The offset is denoted as Q, while $\alpha$ is the clock skew. In this paper, the gfGn model is applied to describe the PDV based on [13,14]. According to [13,14], the PDV variance for $n=1,2$ is

$$\begin{array}{cc}\hfill & E\left[{\omega }_n\left[j\right]{\omega }_n\left[i\right]\right]={\sigma }_{{\omega }_n}^2\mathrm{for}i=j\hfill \\ \hfill & \\ \hfill & \mathrm{and}\hfill \\ \hfill & \\ \hfill & \mathrm{for}i\ne j\mathrm{and}p=F,R(\mathrm{where}0.5\le H_F,H_R<1,\mathrm{and}0<a_F,a_R\le 1),\hfill \\ \hfill & E\left[{\omega }_n\left[j\right]{\omega }_n\left[i\right]\right]={\displaystyle \frac{{\sigma }_{{\omega }_n}^2}{2}}\left({\left|\right|(j-i)}^{a_p}{|-1|}^{2H_p}{-2|(j-i)}^{a_p}{|}^{2H_p}+{\left(\right|(j-i)}^{a_p}{|+1)}^{2H_p}\right)\hfill \end{array}$$

(2)

The expectation operator on (.) is denoted as $E[.]$.

Based on (2) and [13,14,17,20,21], we notice the following cases:

• The fractional Gaussian noise (fGn) case: $a_p=1$.
• The Gaussian case: $a_p=1$ and $H_p=0.5$.

Therefore, the gfGn model answers on a wider range of cases. The newly derived low-computational burden approximations for the MSE presented in Section 3 are for the clock skew estimators [18,19], where the TWD clock skew estimator from [18] is

$$\hat{\alpha }={\displaystyle \frac{1}{J(J-1)}}\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\left({\displaystyle \frac{T_{1,j}\left(i\right)}{T_{2,j}\left(i\right)}}+{\displaystyle \frac{T_{4,j}\left(i\right)}{T_{3,j}\left(i\right)}}\right)-1$$

(3)

where the anticipated clock skew is $\hat{\alpha }$. $T_{l,j}\left(i\right)$ is given for $l=1,2,3,4$ by

$$T_{l,j}\left(i\right)=t_l[j+i]-t_l\left[j\right]$$

(4)

For both the Forward and Reverse pathways, the OWD clock skew estimators from [19] are

$${\hat{\alpha }}_F={\displaystyle \frac{2}{J(J-1)}}\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\left({\displaystyle \frac{T_{1,j}\left(i\right)}{T_{2,j}\left(i\right)}}\right)-1$$

(5)

and

$${\hat{\alpha }}_R={\displaystyle \frac{2}{J(J-1)}}\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\left({\displaystyle \frac{T_{4,j}\left(i\right)}{T_{3,j}\left(i\right)}}\right)-1,$$

(6)

where ${\hat{\alpha }}_F$ and ${\hat{\alpha }}_R$ are, respectively, the expected clock skews for the Forward and Reverse pathways. For the TWD clock skew estimator (3) applicable in the gfGn environment, the closed-form estimates that approximates the MSE is, according to [18], as follows:

$$E\left[e_{{gfGn}_T}^2\right]\approx {\left({\displaystyle \frac{1}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2+{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}\right)\left(\left(1+{\displaystyle \frac{1}{P}}\right)C+D\right)$$

(7)

with P, C, and D given by (10), (13), and (14), respectively.

The MSE [19] associated with the OWD clock skew estimators (5) and (6) that apply in the gfGn environment have, according to [19], the following closed-form estimates.

For the Forward case,

$$E\left[e_{{gfGn}_F}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2}{T_{syn}^2}}\right)\left(\left(1+{\displaystyle \frac{1}{P_F}}\right)C+D\right)$$

(8)

For the Reverse case,

$$E\left[e_{{gfGn}_R}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}\right)\left(C+D\right),$$

(9)

where P and $P_F$ are defined by [18,19] as

$$P={\displaystyle \frac{A}{B}}\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2+{\sigma }_{{\omega }_2}^2}{{\sigma }_{{\omega }_1}^4}}\right)T_{syn}^2;P_F={\displaystyle \frac{A}{B}}\left({\displaystyle \frac{T_{syn}^2}{{\sigma }_{{\omega }_1}^2}}\right)$$

(10)

A, B, C, and D ([18]) are

$$\begin{array}{cc}\hfill & A=(2\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}+\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}-\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}})\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & B=(12\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^4}}+6\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{{\left(ik\right)}^2}}+4\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m\ne j\\ m\ne j+i-k\\ m\ne j+i\\ m\ne j-k\end{array}}^{J-k}{\displaystyle \frac{1}{{\left(ik\right)}^2}})\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & C=\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}\left(2-f{G}_H^{*}(i,H,a)\right)\hfill \\ \hfill & +\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(1+{\displaystyle \frac{1}{2}}\left(f{G}_{\begin{array}{c}H\end{array}}^{*}(i-k,H,a)-f{G}_H^{*}(i,H,a)-f{G}_H^{*}(k,H,a)\right)\right)\hfill \\ \hfill & -\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(1-{\displaystyle \frac{1}{2}}\left(f{G}_H^{*}(i,H,a)-f{G}_H^{*}(k,H,a)+f{G}_H^{*}(i+k,H,a)\right)\right)\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & D=\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m\ne j\\ m\ne j+i\\ m\ne j-k\\ m\ne j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{2ik}}(f{G}_H^{*}(j-m,H,a)\hfill \\ \hfill & -f{G}_H^{*}(j+i-m,H,a)-f{G}_H^{*}(j-m-k,H,a)+f{G}_H^{*}(j+i-m-k,H,a))\hfill \end{array}$$

(14)

and $f{G}_H^{*}(.)$ is, according to [18], as follows:

$$f{G}_H^{*}(y,H,a)=\left[\left|\right|y^a{|-1|}^{2H}-2\left(\right|y^a{\left|\right)}^{2H}+\left(\right|y^a{|+1)}^{2H}\right]$$

(15)

The reader may refer to [18,19] for the proof of (3), (5)–(9). Please note that for the Gaussian case ($H=0.5$ with $a=1$), the function $f{G}_H^{*}(y,0.5,1)$ is zero. Thus, we have $C=A$ and $D=0$ (please refer to [18,19]). The closed-form estimates that approximate the MSE ([18,19]) for the TWD and OWD for the Forward and Reverse paths for the white-Gaussian environment are denoted in this paper as $E\left[e_{G_T}^2\right]$, $E\left[e_{G_F}^2\right]$, and $E\left[e_{G_R}^2\right]$, respectively. These MSE expressions are related to the three clock skew estimators ((3), (5), and (6)) and are given by

$$\begin{array}{c}E\left[e_{G_T}^2\right]\approx {\left({\displaystyle \frac{1}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2+{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}A\right)\left(1+{\displaystyle \frac{1}{P}}\right)\\ \\ E\left[e_{G_F}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left[{\displaystyle \frac{{\sigma }_{{\omega }_1}^2}{T_{syn}^2}}A(1+{\displaystyle \frac{1}{P_F}})\right]\\ \\ E\left[e_{G_R}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left[{\displaystyle \frac{{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}A\right]\end{array}$$

(16)

3. The Low-Computational Burden Approximated Expressions for the MSE

The closed-form estimates that approximate the MSE under the gfGn model in (7)–(9) necessitate significant computational burden resources for each specific set of Hurst exponent and a parameters. Therefore, a need was seen in having closed-form approximated expressions for the MSE that have a lower computational burden for each set of Hurst exponent and a parameters values.

Theorem 1. 

We apply the following assumptions:

1. 

The total number of Sync messages is set to 500. This is based on our previous work [18,19], where 500 Sync messages result in an accuracy of ${10}^{-12}$ from the MSE perspective. Thus, for consistency and comparability, we adopt the same number in this paper.

2. 

The a parameter is in the range of $a\le 0.1$. With the help of this assumption, (15) can be simplified so that the set of Hurst exponent and a parameters are no longer in a nested loop in the MSE calculations. Thus, the MSE result can be obtained faster when changing the Hurst exponent and a parameters from one set to another.

In the following, we denote the closed-form estimates that approximate the MSE for the TWD, OWD for the Forward path, and OWD for the Reverse path under the gfGn environment for $a\le 0.1$, as $E\left[e_{T,(a\le 0.1)}^2\right]$, $E\left[e_{F,(a\le 0.1)}^2\right]$, and $E\left[e_{R,(a\le 0.1)}^2\right]$, respectively. These closed-form estimates that approximate the MSE associated with the three clock skew estimators ((3),(5), and (6)) are given by

$$E\left[e_{T,(a\le 0.1)}^2\right]\approx {\left({\displaystyle \frac{1}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2+{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}\right)\left(\left(1+{\displaystyle \frac{1}{P}}\right)\tilde{C}+\tilde{D}\right)$$

(17)

$$E\left[e_{F,(a\le 0.1)}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_1}^2}{T_{syn}^2}}\right)\left(\left(1+{\displaystyle \frac{1}{P_F}}\right)\tilde{C}+\tilde{D}\right)$$

(18)

$$E\left[e_{R,(a\le 0.1)}^2\right]\approx {\left({\displaystyle \frac{2}{\left(J\right(J-1\left)\right)}}\right)}^2\left({\displaystyle \frac{{\sigma }_{{\omega }_2}^2}{T_{syn}^2}}\right)\left(\tilde{C}+\tilde{D}\right)$$

(19)

where $\tilde{C}$ is

$$\begin{array}{cc}\hfill & \tilde{C}\approx A\left(2-0.5\left(2^{2H}\right)\right)+H\left(a^{2H}V-4a+2^{2H}a\right)\left(-C_1+{\displaystyle \frac{C_2}{2}}-{\displaystyle \frac{C_3}{2}}\right)\hfill \end{array}$$

(20)

$C_1$, $C_2$, and $C_3$ are

$$\begin{array}{c}\hfill C_1=\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}\left(ln\left(i\right)\right)\end{array}$$

(21)

$$\begin{array}{c}\hfill C_2=\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(ln\left(\right|i-k\left|\right)-ln\left(i\right)-ln\left(k\right)\right)\end{array}$$

(22)

$$\begin{array}{c}\hfill C_3=\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(ln\left(i\right)-ln\left(k\right)+ln(i+k)\right)\end{array}$$

(23)

where $\tilde{D}$ is defined as

$$\begin{array}{cc}\hfill & \tilde{D}\approx H\left(a^{2H}V-4a+2^{2H}a\right)D_1\hfill \end{array}$$

(24)

and $D_1$ is

$$\begin{array}{cc}\hfill D_1=& \sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m\ne j\\ m\ne j+i\\ m\ne j-k\\ m\ne j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{2ik}}\left(ln\right(|j-m\left|\right)-ln\left(\right|j+i-m\left|\right)\hfill \\ \hfill & -ln\left(\right|j-m-k\left|\right)+ln\left(\right|j+i-m-k\left|\right))\hfill \end{array}$$

(25)

Comments: Although our derivation was based on the assumption that $a\le 0.1$, we noticed, via simulation results, that they are also valid for a wider range of a, namely, for $a\le 0.6$. Additionally, please note that according to (20) and (24), it may be easy to calculate the MSE for each Hurst exponent and a parameters, since the expressions (A, $C_1$, $C_2$, $C_3$, and $D_1$) that require computational resources are not dependent on the parameters (H and a). Therefore, the expressions in (20) and (24) can be easily adjusted to each Hurst exponent and a parameters through simple multiplications.

Proof of Theorem 1. 

We use, in the following, parameter $T_D$ to signify the difference between the numbers of $n-th$ and $m-th$ Sync messages ($T_D=m-n$). For example, when $n=1$ and $m=500$, indicating the difference between the first and last Sync messages, $T_D=499$. Similarly, for the difference between two contiguous Sync message numbers, with $m=n+1$, we have $T_D=1$. Before we start our proof, we use some series expansion from [22]. The first series expansion [22] is

$${\left(T_D\right)}^a=e^{ln{\left(T_D\right)}^a}=e^{a\left(ln\left(T_D\right)\right)}\approx 1+a\left(ln\left(T_D\right)\right)+{\displaystyle \frac{{\left(a\left(ln\left(T_D\right)\right)\right)}^2}{2!}}+{\displaystyle \frac{{\left(a\left(ln\left(T_D\right)\right)\right)}^3}{3!}}+\dots$$

(26)

where $ln$ is the $lan$ function, and the ! is the factorial function. Please note that the series expansion is conducted for fixed $T_D$.

The second is the binomial series [22]:

$${(1+\overline{x})}^{2H}\approx 1+2H\left(\overline{x}\right)+{\displaystyle \frac{2H(2H-1)}{2!}}\left({\overline{x}}^2\right)+{\displaystyle \frac{2H(2H-1)(2H-2)}{3!}}\left({\overline{x}}^3\right)+\dots$$

(27)

where $|\overline{x}|<1$.

In the first step, we simplify (15) based on assumptions 1–2. We approximate ${\left(T_D\right)}^a$ only using the two first parts of the series expansion (26); therefore, (15) can be written as

$$f{G}_H^{*}(T_D,H,a)\approx \left[{\left(a\left(ln\left(T_D\right)\right)\right)}^{2H}-2{\left(1+a\left(ln\left(T_D\right)\right)\right)}^{2H}+{\left(2+a\left(ln\left(T_D\right)\right)\right)}^{2H}\right]$$

(28)

In the second step, we simplify the three parts of the right expression in (28), with the first two parts of the binomial series (27). Thus, ${(1+a\left(ln\left(T_D\right)\right))}^{2H}$ from (28) can be defined as

$${(1+a\left(ln\left(T_D\right)\right))}^{2H}\approx 1+2H\left(a\left(ln\left(T_D\right)\right)\right)$$

(29)

Based on (29), the expression ${(2+a\left(ln\left(T_D\right)\right))}^{2H}$ in (28), can be written as

$${(2+a\left(ln\left(T_D\right)\right))}^{2H}=2^{2H}{\left(1+{\displaystyle \frac{a\left(ln\left(T_D\right)\right)}{2}}\right)}^{2H}\approx 2^{2H}(1+H\left(a\left(ln\left(T_D\right)\right)\right))$$

(30)

The expression ${\left(a\left(ln\left(T_D\right)\right)\right)}^{2H}$ from (28) is more difficult to approximate. In this paper, we use the following approximation:

$${\left(a\left(ln\left(T_D\right)\right)\right)}^{2H}=a^{2H}{\left(ln\left(T_D\right)\right)}^{2H}\approx a^{2H}H\left(ln\left(T_D\right)\right)V$$

(31)

where V is a constant correction factor to the power of $ln\left(T_D\right)$, as given in

Table 1. V values for $1<T_D<500$.

			Hurst	Exponent	Parameter
		0.5	0.6	0.7	0.8	0.9
a parameter	0.6	2	1.8	1.45	0.85	0.3
	0.5	2	1.85	1.5	1.05	0.3
	0.1	2	2.2	2.1	1.2	1
	0.05	2	1.9	1.8	1.2	0.25
	0.01	2	2.1	2.1	2	2
	1 $\times {10}^{-6}$			Irrelevant

. Based on (29)–(31), we can write (28) as

$$\begin{array}{cc}\hfill f{G}_H^{*}(T_D,H,a\le 0.1)& \approx \left[\left(a^{2H}H\left(ln\left(T_D\right)\right)V\right)-2\left(1+2Ha\left(ln\left(T_D\right)\right)\right)+2^{2H}\left(1+Ha\left(ln\left(T_D\right)\right)\right)\right]\hfill \\ \hfill & \approx (-2+2^{2H})+ln\left(T_D\right)\left(a^{2H}HV-4Ha+2^{2H}Ha\right)\hfill \end{array}$$

(32)

Please note that determining the optimal value for V is challenging due to its dependence on $T_D$. It is essential to ascertain the best value for V relative to the varying $T_D$ values. Thus, we conducted simulations to identify the approximated V value that would produce the lowest difference between the results obtained from the approximated expressions for the MSE and those obtained by the clock skew estimators ((3), (5), and (6)) performances from the MSE perspective for $1<T_D<500$. For each combination of a and H parameters, a corresponding V value was determined across different $T_D$ values. The V values for $1<T_D<500$ are given in Table 1. The last row of the Table 1 reads ’Irrelevant’, since the magnitude of V does not impact the outcome when the a parameter approaches zero. In the simulation, as the a parameter tended toward zero, V was maintained at a fixed value.

Based on (32), we can calculate $\tilde{C}$ and $\tilde{D}$, where $\tilde{C}=C(J,H,a\le 0.1)$ and $\tilde{D}=D(J,H,a\le 0.1)$. First, by putting (32) into (13), we have

$$\begin{array}{cc}\hfill & \tilde{C}\approx [\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}\left(2-(-2+2^{2H})-ln\left(i\right)\left(a^{2H}HV-4Ha+2^{2H}Ha\right)\right)\hfill \\ \hfill & +\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}(1+{\displaystyle \frac{1}{2}}((2-2^{2H})+\left(a^{2H}HV-4Ha+2^{2H}Ha\right)\hfill \\ \hfill & \left(ln\right(i-k\left)\right)-\left(ln\right(i\left)\right)-\left(ln\right(k\left)\right)\left)\right)\hfill \\ \hfill & -\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}(1-{\displaystyle \frac{1}{2}}((-2+2^{2H})+\left(a^{2H}HV-4Ha+2^{2H}Ha\right)\hfill \\ \hfill & \left(ln\right(i\left)\right)-\left(ln\right(k\left)\right)+\left(ln\right(i+k\left)\right)\left)\right)]\hfill \end{array}$$

(33)

After rearranging (33), we have the same expression as (20). For $\tilde{D}$, we substitute (32) into (14). Thus, we have the same expression as (24), and this completes our proof. □

Please note that for the Gaussian case, the derived closed-form estimates that approximate MSEs (17), (18), and (19) are almost the same as the closed-form estimates that approximate the MSE given in (16) for the Gaussian case. In order to show this, please refer to Table 1 and (20). Based on Table 1, where the correction factor V is set to 2, and (20), $\tilde{C}$ can be written as

$$\begin{array}{cc}\hfill & \tilde{C}\approx A\left(2-0.5\left(2\right)\right)+H\left(2a-4a+2a\right)\left(-C_1+{\displaystyle \frac{C_2}{2}}-{\displaystyle \frac{C_3}{2}}\right)\approx A\hfill \end{array}$$

(34)

Next, $\tilde{D}$ may be written based on (24) as

$$\begin{array}{cc}\hfill & \tilde{D}\approx H\left(2a-4a+2a\right)D_1\approx 0\hfill \end{array}$$

(35)

that is, we approximated those equations for the MSE valid for the Gaussian case:

$$\begin{array}{c}E{\left[e_{T,(a\le 0.1)}^2\right]}_{(H_F=H_R=0.5)}\approx E\left[e_{G_T}^2\right]\\ \\ E{\left[e_{F,(a\le 0.1)}^2\right]}_{(H_F=0.5)}\approx E\left[e_{G_F}^2\right]\\ \\ E{\left[e_{R,(a\le 0.1)}^2\right]}_{(H_R=0.5)}\approx E\left[e_{G_R}^2\right]\end{array}$$

(36)

Next, we deal with the case of a tending to zero.

The New Obtained Expressions for the MSE When a Is Tending to Zero

A notable scenario arises when the parameter a tends toward zero ($a\to 0$). In this context, as the Hurst exponent parameter increases and the a parameter tends toward zero, the closed-form estimates that approximate the MSE for the non-Gaussian cases exhibit superior performances compared to the performances obtained with the MSE expressions valid for the Gaussian case. In order to show this, we assume first that ${\left(T_D\right)}^a\approx 1$ for $a\to 0$. Therefore, (15) can be defined as

$$f{G}_H^{*}(x,H,0)\approx \left[-2+2^{2H}\right]$$

(37)

Now, by putting (37) into (13) and (14), we can write

$$\begin{array}{cc}\hfill & C(J,H,0)\approx \sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}\left(2-(-2+2^{2H})\right)\hfill \\ \hfill & +\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(1+{\displaystyle \frac{1}{2}}\left(-(-2+2^{2H})\right)\right)\hfill \\ \hfill & -\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}\left(1-{\displaystyle \frac{1}{2}}\left((-2+2^{2H})\right)\right)\hfill \end{array}$$

(38)

and

$$\begin{array}{cc}\hfill & D(J,H,0)\approx 0\hfill \end{array}$$

(39)

After rearranging (38), we have:

$$\begin{array}{cc}\hfill & C(J,H,0)\approx (2-0.5\left(2^{2H}\right))\hfill \\ \hfill & (\sum _{i=1}^{J-1}{\displaystyle \frac{J-i}{i^2}}2+\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\\ k\ne i\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j\\ m=j+i-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}}-\sum _{i=1}^{J-1}\sum _{j=1}^{J-i}\sum _{\begin{array}{c}k=1\end{array}}^{J-1}\sum _{\begin{array}{c}m=1\\ m=j+i\\ m=j-k\end{array}}^{J-k}{\displaystyle \frac{1}{ik}})\hfill \\ \hfill & \approx (2-0.5\left(2^{2H}\right))A\hfill \end{array}$$

(40)

Please note that for $a\to 0$ and based on (20), we may write $\tilde{C}$ as

$$\begin{array}{cc}\hfill \tilde{C}& \approx {\left[A\left(2-0.5\left(2^{2H}\right)\right)+\left(-C_1+{\displaystyle \frac{C_2}{2}}-{\displaystyle \frac{C_3}{2}}\right)H\left(a^{2H-1}V-4+2^{2H}\right)a\right]}_{a\to 0}\hfill \\ \hfill & \approx A(2-0.5\left(2^{2H}\right))\hfill \end{array}$$

(41)

Based on (24) and for $a\to 0$, we may write $\tilde{D}$ as

$$\tilde{D}\approx H{D}_1\left(a^{2H-1}V-4+2^{2H}\right)a{|}_{a\to 0}\approx 0$$

(42)

Next, based on (37)–(42), the closed-form estimates that approximate the MSE under the gfGn environment for $a\to 0$ can be written as

$$E{\left[e_{gfG{n}_T}^2\right]}_{(a_F,a_R\to 0)}\approx E{\left[e_{T,(a\le 0.1)}^2\right]}_{(a_F,a_R\to 0)}\approx E\left[e_{G_T}^2\right](2-0.5\left(2^{2H}\right))$$

(43)

$$E{\left[e_{gfG{n}_F}^2\right]}_{(a_F\to 0)}\approx E{\left[e_{F,(a\le 0.1)}^2\right]}_{(a_F\to 0)}\approx E\left[e_{G_F}^2\right](2-0.5\left(2^{2H}\right))$$

(44)

$$E{\left[e_{gfG{n}_R}^2\right]}_{(a_R\to 0)}\approx E{\left[e_{R,(a\le 0.1)}^2\right]}_{(a_R\to 0)}\approx E\left[e_{G_R}^2\right](2-0.5\left(2^{2H}\right))$$

(45)

For example, let us set $H_F=0.9$ and $a_F\to 0$. Then, the closed-form approximated expressions for the MSE for the Forward path is $E\left[e_{G_F}^2\right](a_F\to 0,H_F=0.9)\approx E\left[e_{G_F}^2\right]\left(0.2589\right)$. Thus, it can be clearly seen in this example that the MSE expression for the case where the a parameter approaches zero is smaller (thus having improved performances) compared to the MSE expression valid for the Gaussian case.

4. Results

In this study, we tested our new closed-form approximated expressions for the MSE (17)–(19) for the gfGn case against the simulated results obtained from the clock skew estimators performances from the MSE perspective. Initially, we present simulation results for various values of the a parameter, specifically those less than or equal to $0.1$, to demonstrate the accuracy of the new MSE expressions for the gfGn case when $a\le 0.1$, including the special case where the a parameter tends to zero. Subsequently, we analyze simulation results for a parameters greater than $0.1$ ($a=0.5,0.6$) to illustrate the applicability of the new MSE expressions for cases where a is less than or equal to $0.6$.

Let us start with the simulation results for $a\le 0.1$. Figure 2, Figure 3 and Figure 4 show the performances of clock skew estimators (5) and (6) compared to the closed-form estimates that approximate the MSE (18) and (19) under the gfGn case for $a\le 0.1$. According to Figure 2, Figure 3 and Figure 4, it can be seen that the performances obtained with the MSE expressions are quite similar to the simulated performances of the clock skew estimators (MSE perspective).

Figure 5 illustrates the performances of the clock skew estimators (5) and (6) compared to the closed-form estimates that approximate the MSE when $a\to 0$: (44) and (45). According to Figure 5, for $a\to 0$, those MSE expressions closely align with the simulated clock skew estimators’ performances. As was noted in Section 3, when the a parameter tends to zero and the Hurst exponent parameter increases, the closed-form estimates that approximate the MSE demonstrate better performances compared to the performances obtained with the MSE expressions valid for the Gaussian case. Indeed, Figure 5 illustrates this phenomenon.

In Section 3, we presented the closed-form estimates that approximate the MSE for the case where the a parameter is less than $0.1$. However, our simulations revealed that it is feasible to identify values for V that yields a reasonable approximation to the MSE expressions, even for values of the a parameter that exceed $0.1$ within the range of $T_D$ values between 1 and 500. This observation underscores the existence of a V value capable of facilitating a robust approximation to (15), even when the a parameter is greater than $0.1$. It is noteworthy that when the a parameter exceeds $0.1$, the closed-form estimates that approximate the MSE become a little bit less accurate compared to the case where $a\le 0.1$. Consequently, the newly derived closed-form estimates that approximate MSEs (17)–(19) remain applicable even for a values below $0.6$. Figure 6, Figure 7 and Figure 8 show the performances of clock skew estimators (3), (5), and (6) compared to the closed-form estimates that approximate MSEs (17)–(19). According to Figure 6, Figure 7 and Figure 8, even for an a parameter greater than 0.1, those MSE expressions obtain quite similar results to the simulated performances of the clock skew estimators (MSE perspective). Thus, although our new proposed approximated expressions for the MSE were obtained for $a\le 0.1$, they can also be applied according to simulation results for $a\le 0.6$.

In the newly derived expressions for MSEs (17)–(19), the Hurst exponent and a parameters are connected with $\tilde{C}$ (20) and $\tilde{D}$ (24). In the previously obtained expression for MSEs (7)–(9), the Hurst exponent and a parameters are connected with C (13) and D (14). Please note that the Hurst exponent and a parameters are in a nested loop in C (13) and D (14) unlike in $\tilde{C}$ (20) and $\tilde{D}$ (24). Thus, when varying the Hurst exponent and a parameters from one set to another, the Matlab code execution time for the previously obtained expression for the MSE is longer than for the newly obtained one. To show this, we tested the Matlab code execution time for $\tilde{C}$, $\tilde{D}$, C, and D when varying the Hurst exponent and a parameters from one set to another. Since $C1$ (21), $C2$ (22), $C3$ (23), and $D1$ (25) are not a function of H and a, they do not have to be calculated again when varying the Hurst exponent and a parameters from one set to another. We used MATLAB Version: 9.7.0.1296695 (R2019b) Update 4; Operating System: Microsoft Windows 10 Pro Version 10.0. In the following test, we changed the Hurst exponent and a parameters from one set to $H=0.8$ and $a=0.6$ with $V=0.85$. We obtained the following results when $C1$ (21), $C2$ (22), $C3$ (23), and $D1$ (25) were not calculated again:

• For the total number of Sync periods equal to $\mathbf{100}$ ($\mathbf{J}=\mathbf{100}$):
• The averaged Matlab code execution time for D is $63.8585$ s. The averaged results were obtained for 10 trials.
• The averaged Matlab code execution time for C is $1.8535$ s. The averaged results were obtained for 10 trials.
• The averaged Matlab code execution time for $\tilde{D}$ is $0.5015$ ms. The averaged results were obtained for 10,000 trials.
• The averaged Matlab code execution time for $\tilde{C}$ is $0.544$ ms. The averaged results were obtained for 10,000 trials.
• For the total number of Sync periods equal to $\mathbf{200}$ ($\mathbf{J}=\mathbf{200}$):
• The averaged Matlab code execution time for D is $1132.7$ s. The averaged results were obtained for 5 trials.
• The averaged Matlab code execution time for C is $17.24$ s. The averaged results were obtained for 5 trials.
• The averaged Matlab code execution time for $\tilde{D}$ is $0.7107$ ms. The averaged results were obtained for 10,000 trials.
• The averaged Matlab code execution time for $\tilde{C}$ is $0.9047$ ms. The averaged results were obtained for 10,000 trials.

Next, we changed the Hurst exponent and a parameters from one set to $H=0.8$ and $a=0.6$ with $V=0.85$. But now, $C1$ (21), $C2$ (22), $C3$ (23), and $D1$ (25) were calculated again. We received the following results:

• For the total number of Sync periods equal to $\mathbf{100}$ ($\mathbf{J}=\mathbf{100}$):
• The averaged Matlab code execution time for $\tilde{D}$ is $2.9462$ s. The averaged results were obtained for 10 trials.
• The averaged Matlab code execution time for $\tilde{C}$ is $0.3125$ s. The averaged results were obtained for 10 trials.
• For the total number of Sync periods equal to $\mathbf{200}$ ($\mathbf{J}=\mathbf{200}$):
• The averaged Matlab code execution time for $\tilde{D}$ is $36.1930$ s. The averaged results were obtained for 5 trials.
• The averaged Matlab code execution time for $\tilde{C}$ is $2.3896$ s. The averaged results were obtained for 5 trials.

According to the results, the newly derived MSE has, approximately, an averaged Matlab code execution time that is thirty times faster compared to the averaged Matlab code execution time of the previously obtained MSE expression. It should be pointed out that this advantage in the Matlab code execution time may be even higher for $J>200$.

5. Discussion

The newly proposed expressions for the MSE applicable to the gfGn case for $a\le 0.1$ offer a significant advantage over our previous MSE expressions due to their reduced computational burden. These MSE expressions can be easily adapted to any pair of Hurst exponent and a parameters with minimum computational resources. Additionally, we derived MSE expressions for the case where the a parameter tends to zero. In this scenario, these newly derived MSE expressions are functions of the obtained MSE expressions related to the Gaussian case and require minimum computational burden. Simulation results revealed that the new closed-form approximations for $a\le 0.1$ also apply to cases where $a\le 0.6$. Thus, these MSE expressions provide a computationally efficient solution for the gfGn case with $a\le 0.6$. In summary, the new expressions can significantly reduce the time spent running simulations, reduce computational load, and lower the requirements for software and hardware. This makes the implementation of PTP systems more efficient and cost-effective.