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.
