The NTP ClockCombining
Algorithm
Introduction
A common problem in synchronization subnets is
systematic timeoffset errors resulting from asymmetric transmission paths, where the networks or transmission media in one direction
are substantially different from the other. The errors can range from microseconds on highspeed ring networks to large fractions of a
second on satellite/landline paths. It has been found experimentally that these errors can be considerably reduced by combining the
apparent offsets of a number of time servers to produce a more accurate working offset. Following is a description of the combining
method used in the NTP implementation for the Fuzzball [MIL88b]. The method is similar to that used by national standards laboratories
to determine a synthetic laboratory time scale from an ensemble of cesium clocks [ALL74b]. These procedures are optional and not
required in a conforming NTP implementation.
In the following description the stability of a
clock is how well it can maintain a constant frequency, the accuracy is how well its frequency and time compare with national
standards and the precision is how precisely these quantities can be maintained within a particular timekeeping system. Unless
indicated otherwise, The offset of two clocks is the time difference between them, while the skew is the frequency difference (first
derivative of offset with time) between them. Real clocks exhibit some variation in skew (second derivative of offset with time),
which is called drift.
Determining Time and
Frequency
Figure 9 shows the overall organization of the NTP
timeserver model. Timestamps exchanged with possibly many other subnet peers are used to determine individual roundtrip delays and
clock offsets relative to each peer as described in the NTP specification. As shown in the figure, the computed delays and offsets are
processed by the clock filter to reduce incidental timing noise and the most accurate and reliable subset determined by the
clockselection algorithm. The resulting offsets of this subset are first combined as described below and then processed by the
phaselocked loop (PLL). In the PLL the combined effects of the filtering, selection and combining operations is to produce a
phasecorrection term. This is processed by the loop filter to control the local clock, which functions as a voltagecontrolled
oscillator (VCO). The VCO furnishes the timing (phase) reference to produce the timestamps used in all calculations.
Clock Modelling
The International Standard (SI) definition of time
interval is in terms of the standard second: the duration of 9,192,631,770 periods of the radiation corresponding to the transition
between the two hyperfine levels of the ground state of the cesium133 atom. Let u represent the standard unit of time interval so
defined and v = 1 over u be the standard unit of frequency. The epoch, denoted by t, is defined as the reading of a counter that runs
at frequency v and began counting at some agreed initial epoch t0, which defines the standard or absolute time scale. For the purposes
of the following analysis, the epoch of the standard time scale, as well as the time indicated by a clock will be considered
continuous. In practice, time is determined relative to a clock constructed from an atomic oscillator and system of counter/dividers,
which defines a time scale associated with that particular oscillator. Standard time and frequency are then determined from an
ensemble of such timescales and algorithms designed to combine them to produce a composite time scale approximating the standard time
scale.
Let T(t) be the time displayed by a clock at epoch t
relative to the standard time scale:
T(t) = 1/2 D(t sub 0 )[t  t sub 0 ] sup 2 + R(t sub
0 )[t  t sub 0 ] + T(t sub 0 ) + x(t),
where D(t sub 0 ) is the fractional frequency drift
per unit time, R(t sub 0 ) the frequency and T(t sub 0 ) the time at some previous epoch t0. In the usual stationary model these
quantities can be assumed constant or changing slowly with epoch. The random nature of the clock is characterized by x(t), which
represents the random noise (jitter) relative to the standard time scale. In the usual analysis the secondorder term D(t sub 0 ) is
ignored and the noise term x(t) modelled as a normal distribution with predictable spectral density or autocorrelation
function.
The probability density function of time offset p (t
 T(t)) usually appears as a bellshaped curve centered somewhere near zero. The width and general shape of the curve are determined
by x(t), which depends on the oscillator precision and jitter characteristics, as well as the measurement system and its transmission
paths. Beginning at epoch t0 the offset is set to zero, following which the bell creeps either to the left or right, depending on the
value of R(t sub 0 ) and accelerates depending on the value of D(t sub 0 ).
Development of a Composite
Timescale
Now consider the time offsets of a number of real
clocks connected by real networks. A display of the offsets of all clocks relative to the standard time scale will appear as a system
of bellshaped curves slowly precessing relative to each other, but with some further away from nominal zero than others. The bells
will normally be scattered over the offset space, more or less close to each other, with some overlapping and some not. The problem is
to estimate the true offset relative to the standard time scale from a system of offsets collected routinely between the
clocks.
A composite time scale can be determined from a
sequence of offsets measured between the n clocks of an ensemble at nominal intervals tau. Let R sub i (t sub 0 ) be the frequency and
T sub i (t sub 0 ) the time of the ith clock at epoch t0 relative to the standard time scale and let "^" designate the associated
estimates. Then, an estimator for Ti computed at t0 for epoch t sub 0 + tau is
neglecting secondorder terms. Consider a set of n
independent timeoffset measurements made between the clocks at epoch t sub 0 + tau and let the offset between clock i and clock j at
that epoch be T sub ij (t sub 0 + tau ), defined as
Note that T sub ij =  T sub ji and T sub ii =
0. Let w sub i ( tau ) be a previously determined weight factor associated with the ith clock for the nominal interval tau. The basis
for new estimates at epoch t sub 0 + tau is
That is, the apparent time indicated by the jth
clock is a weighted average of the estimated time of each clock at epoch t sub 0 + tau plus the time offset measured between the jth
clock and that clock at epoch t sub 0 + tau.
An intuitive grasp of the behavior of this algorithm
can be gained with the aid of a few examples. For instance, if w sub i ( tau ) is unity for the ith clock and zero for all others, the
apparent time for each of the other clocks is simply the estimated time T hat sub i (t sub 0 + tau ). If w sub i ( tau ) is zero for
the ith clock, that clock can never affect any other clock and its apparent time is determined entirely from the other clocks.
If w sub i ( tau ) = 1 / n for all i, the apparent time of the ith clock is equal to the average of the time estimates computed
at t0 plus the average of the time offsets measured to all other clocks. Finally, in a system with two clocks and w sub i ( tau ) = 1
/ 2 for each, and if the estimated time at epoch t sub 0 + tau is fast by 1 s for one clock and slow by 1 s for the other, the
apparent time for both clocks will coincide with the standard time scale.
In order to establish a basis for the next interval
tau, it is necessary to update the frequency estimate R hat sub i (t sub 0 + tau ) and weight factor w sub i ( tau ). The average
frequency assumed for the ith clock during the previous interval tau is simply the difference between the times at the beginning and
end of the interval divided by tau. A good estimator for R sub i (t sub 0 + tau ) has been found to be the exponential average of
these differences, which is given by
R hat sub i (t sub 0 + tau ) = R hat sub i (t sub 0
) + alpha sub i [ R hat sub i (t sub 0 )  {T sub i (t sub 0 + tau )  T sub i (t sub 0 )} over tau ],
where alpha sub i is an experimentally
determined weight factor which depends on the estimated frequency error of the ith clock. In order to calculate the weight factor w
sub i ( tau ), it is necessary to determine the expected error epsilon sub i ( tau ) for each clock. In the following, braces ""
indicate absolute value and brackets "<>" indicate the infinite time average. In practice, the infinite averages are computed as
exponential time averages. An estimate of the magnitude of the unbiased error of the ith clock accumulated over the nominal interval
tau is
where epsilon sub i ( tau ) and epsilon sub e ( tau
) are the accumulated error of the ith clock and entire clock ensemble, respectively. The accumulated error of the entire ensemble
is
Finally, the weight factor for the ith clock is
calculated as
When all estimators and weight factors have been
updated, the origin of the estimation interval is shifted and the new value of t0 becomes the old value of t sub 0 + tau.
While not entering into the above calculations, it
is useful to estimate the frequency error, since the ensemble clocks can be located some distance from each other and become isolated
for some time due to network failures. The frequencyoffset error in Ri is equivalent to the fractional frequency yi,
y sub i
= { nu sub i~~nu sub I } over nu sub I
measured between the ith time scale and the standard
time scale I. Temporarily dropping the subscript i for clarity, consider a sequence of N independent frequencyoffset samples y(j) (j=
1, 2,... ,N) where the interval between samples is uniform and equal to T. Let tau be the nominal interval over which these samples
are averaged. The Allan variance sigma sub y sup 2 ( N,T,tau ) [ALL74a] is defined as
A particularly useful formulation is N = 2 and T =
tau:
so that
While the Allan variance has found application when
estimating errors in ensembles of cesium clocks, its application to NTP is limited due to the computation and storage burden. As
described in the next section, it is possible to estimate errors with some degree of confidence using normal byproducts of NTP
processing algorithms.
Application to NTP
The NTP clock model is somewhat less complex than
the general model described above. For instance, at the present level of development it is not necessary to separately estimate the
time and frequency of all peer clocks, only the time and frequency of the local clock. If the timekeeping reference is the local clock
itself, then the offsets available in the peer.offset peer variables can be used directly for the T sub ij quantities above. In
addition, the NTP localclock model incorporates a typeII phaselocked loop, which itself reliably estimates frequency errors and
corrects accordingly. Thus, the requirement for estimating frequency is entirely eliminated.
There remains the problem of how to determine a
robust and easily computable error estimate epsilon sub i. The method described above, although analytically justified, is most
difficult to implement. Happily, as a byproduct of the NTP clockfilter algorithm, a useful error estimate is available in the form of
the dispersion. As described in the NTP specification, the dispersion includes the absolute value of the weighted average of the
offsets between the chosen offset sample and the n1 other samples retained for selection. The effectiveness of this estimator was
compared with the above estimator by simulation using observed timekeeping data and found to give quite acceptable results.
The NTP clockcombining algorithm can be implemented
with only minor modifications to the algorithms as described in the NTP specification. Although elsewhere in the NTP specification the
use of generalpurpose multiply/divide routines has been successfully avoided, there seems to be no way to avoid them in the
clockcombining algorithm. However, for best performance the localclock algorithm described elsewhere in this document should be
implemented as well, since the combining algorithms result in a modest increase in phase noise which the revised localclock algorithm
is designed to suppress.
ClockCombining
Procedure
The result of the NTP clockselection procedure is a
set of survivors (there must be at least one) that represent truechimers, or correct clocks. When clock combining is not implemented,
one of these peers, chosen as the most likely candidate, becomes the synchronization source and its computed offset becomes the final
clock correction. Subsequently, the system variables are adjusted as described in the NTP clockupdate procedure. When clock combining
is implemented, these actions are unchanged, except that the final clock correction is computed by the clockcombining
procedure.
The clockcombining procedure is called from the
clockselect procedure. It constructs from the variables of all surviving peers the final clock correction THETA. The estimated error
required by the algorithms previously described is based on the synchronization distance LAMBDA computed by the distance procedure, as
defined in the NTP specification. The reciprocal of LAMBDA is the weight of each clockoffset contribution to the final clock
correction. The following pseudocode describes the procedure.
begin clockcombining procedure
temp1 <
0;
temp2 <
0;
for (each peer
remaining on the candidate list) /* scan all survivors */
LAMBDA < distance (peer);
temp < 1 over peer.stratum * NTP.MAXDISPERSE + LAMBDA;
temp1 < temp1 + temp; /* update weight and offset
*/
temp2 < temp2 + temp * peer.offset;
endif;
THETA <
temp2 over
temp1;
/* compute final correction */
end clockcombining procedure;
The value THETA is the final clock correction used
by the localclock procedure to adjust the clock.
