WITH) STATES APARTMENT OF COMMERCE

PUBLICATION NBS TECHNICAL NOTE 636

Modeling of

Atomic Clock Performance

and Detection of

Abnormal Clock Behavior

U.S. APARTMENT

OF

COMMERCE

National

Rureau

loo Js U51S3

i<?9 3
NATIONAL BUREAU OF STANDARDS

The National Bureau of Standards1 was established by an act of Congress March 3,

1 901 . The Bureau's overall goal is to strengthen and advance the Nation's science and

technology and facilitate their effective application for public benefit. To this end, the Bureau conducts research and provides: (1) a basis for the Nation's physical measure- ment system, (2) scientific and technological services for industry and government, (3) a technical basis for equity in trade, and (4) technical services to promote public safety. The Bureau consists of the Institute for Basic Standards, the Institute for Materials Research, the Institute for Applied Technology, the Center for Computer Sciences and

Technology, and the Office for Information Programs. THE INSTITUTE FOR BASIC STANDARDS provides the central basis within the United States of a complete and consistent system of physical measurement; coordinates

that system with measurement systems of other nations; and furnishes essential services leading to accurate and uniform physical measurements throughout the Nation's scientific community, industry, and commerce. The Institute consists of a Center for Radiation Research, an Office of Measurement Services and the following divisions: Applied Mathematics — Electricity — Mechanics — Heat — Optical Physics —

Linac Radiation - — Nuclear Radiation 2 — Applied Radiation - — Quantum

Electronics :1 — Electromagnetics : ' — Time and Frequency 3 — Laboratory

Astrophysics ;; — Cryogenics 3

. THE INSTITUTE FOR MATERIALS RESEARCH conducts materials research lead- ing to improved methods of measurement, standards, and data on the properties of

well-characterized materials needed by industry, commerce, educational institutions, and Government; provides advisory and research services to other Government agencies; and develops, produces, and distributes standard reference materials. The Institute consists of the Office of Standard Reference Materials and the following divisions: Analytical Chemistry—Polymers—Metallurgy—Inorganic Materials—Reactor

Radiation—Physical Chemistry.

THE INSTITUTE FOR APPLIED TECHNOLOGY provides technical services to pro- mote the use of available technology and to facilitate technological innovation in industry and Government; cooperates with public and private organizations leading to the development of technological standards (including mandatory safety standards), codes and methods of test; and provides technical advice and services to Government agencies upon request. The Institute also monitors NBS engineering standards activities and

provides liaison between NBS and national and international engineering standards

bodies. The Institute consists of a Center for Building Technology and the following

divisions and offices: Engineering and Product Standards—Weights and Measures—Invention and

Innovation—Product Evaluation Technology—Electronic Technology—Technical Analysis—Measurement Engineering—Building Standards and Code Services4—Housing Technology 1 —Federal Building Technology4—Structures, Materials and Life Safety 4—Building Environment4—Technical Evaluation and

Application4—Fire Technology.

THE INSTITUTE FOR COMPUTER SCIENCES AND TECHNOLOGY conducts re- search and provides technical services designed to aid Government agencies in improving cost effectiveness in the conduct of their programs through the selection, acquisition, and effective utilization of automatic data processing equipment; and serves as the principal focus within the executive branch for the development of Federal standards for automatic data processing equipment, techniques, and computer languages. The Center

consists of the following offices and divisions: Information Processing Standards—Computer Information—Computer Services —Systems Development—Information Processing Technology.

THE OFFICE FOR INFORMATION PROGRAMS promotes optimum dissemination and accessibility of scientific information generated within NBS and other agencies of

the Federal Government; promotes the development of the National Standard Reference Data System and a system of information analysis centers dealing with the broader

aspects of the National Measurement System; provides appropriate services to ensure

that the NBS staff has optimum accessibility to the scientific information of the world, and directs the public information activities of the Bureau. The Office consists of the following organizational units: Office of Standard Reference Data—Office of Technical Information and

Publications—Library—Office of International Relations.

1 Headquarters and Laboratories at Gaithersburg, Maryland, unless otherwise noted ; mailing address Washington, D.C. 20234.

2 Part of the Center (or Radiation Research.

3 Located at Boulder, Colorado 80302.

' Part o£ the Center for Building Technology.
pNAL U

IS I ^,

Modeling of Atomic Clock Performance 753

• t>?6> and Detection of Abnormal Clock Behavior

73

C. 2-

William A. Gahter

Department of Engineering Design and Economic Evaluation

University of Colorado

Final Report on Research Sponsored by Time and Frequency Division

Institute ibr Basic Standards

"T4>i/> National Bureau of Standards

Boulder, Colorado 80302

NBS Technical notes are designed to supplement the

Bureau's regular publications program. They provide

a means for making available scientific data that are of transient or limited interest. Technical Notes may

be listed or referred to in the open literature.

**'3^7

U.S. DEPARTMENT OF COMMERCE, Frederick B. Dent, Secretary

NATIONAL BUREAU OF STANDARDS, Richard W. Roberts, Director

Issued March 1973

n^i
National Bureau of Standards Technical Note 636

Nat. Bur. Stand. (I.S.). Tech. Note 636, 39 pages (March 1973) C0DE1N: NBTINAE

For gale by the Superintendent of Documents, U.S. Government Printing Office, Washington, D. C. 20402

(Order l>\ SI) Catalog No. < 13.16:636)
Contents

Page

Introduction 1 Noise Model 4

Structure of the Detection Procedures 9 Experimental Development and Evaluation 13

4. 1. Development of the Detection Procedures 13

4.2. Evaluation of the Detection Procedures 14

Conclusions 29

References 30
List of Figures

Figure Page

1 A nominal run with (a, b) = (. 7, . 3) and no 17

detection in 512 days.

2 A nominal run with (a, b) = (. 7, . 3) and no 18

detection in 512 days.

3 A nominal run with (a, b) = (. 8, . 2) and no 18

detection in 512 days.

4 A nominal run with (a, b) = (.8,-2) and no 19

detection in 512 days.

5 A nominal run with (a, b) = (. 8, . 2) and no 1 9 detection in 512 days.

6 Detection of a J = -3.4 jump by eq (23) when 20

(a, b) = (. 8, .2) and t = 35.

c 7 Detection of a J = -3. 2 jump by eq (23) when 20

(a, b) = (. 8, . 2) and t = 34.

c

8 Detection of a J = 3. 5 jump by eq (23) when 21

(a, b) = (.7, . 3) and t = 25.

c

9 Detection of a J = 3. jump by eqs (24), (25), and 21

(26) when (a, b) - (. 7, . 3) and t = 35.

c 10 Detection of a D = 1.0 time drift by eq (31) when 22

(a, b) = (. 8, . 2) and t = 35.

c 11 Detection of a D = 1 . 5 time drift by eq (31) when 22

(a, b) - (.7, . 3) and t = 31.

c 12 Detection of a D = -. 1, Q = -. 02 time and rate drift 23

by eq (32) when (a, b) - (. 7, . 3) and t = 28.

13 Detection of a D = . 2, Q = . 05 time and rate drift 2 3 by eq (31) when (a, b) = (.7, .3) and t =320.

14 Detection at day 126 of noise change jS = . 6 when 24

(a, b) = (. 7, . 3) and t = 28.

c 15 Detection at day 301 of noise change /3 = . 1 when 24

(a, b) = (. 7, . 3) and t = 26.

c 16 Detection at day 271 of noise change (8 = . 5 when 2 5

(a, b) = (. 7, . 3) and t = 28.

c

IV
List of Figures (cont'd. )

Figure Page

17 Detection at day 1 1 1 of noise change 0! = . 9 when 2 5

(a, b) = (. 7, . 3) and t = 24.

c 18 Detection at day 216 of noise change Ct = . 5 when 26

(a, b) = (. 7, . 3) and t = 28.

c 19 False inference at day 111 that a > .7 when actually 26

j8 = . 5. (a, b) = (. 7, . 3) and t = 26.

c 20 False inference at day 81 that Ci > . 7 when actually 27

|8 = .1. (a, b) = (.7, . 3) and t = 31.

c 21 False detection of a jump by eq (23) at day 164 when 27

actually j8 = . 5. (a, b) = (. 7, . 3) and t = 26.

22 False detection of a drift by eq (31) at day 316 when 28

actually jS = . 4. (a, b) = (. 7, . 3) and t =28.

23 False detection of a drift by eq (31) at day 276 when 28

actually (8 = . 4. (a, b) = (. 8, . 2) and t = 29.

c

Table

List of Tables

Function sigma (a, b, T)

Modeling of Atomic Clock Performance

and Detection of Abnormal Clock Behavior

William A. Ganter

We have assumed that the nominal performance of an

atomic clock can be well characterized by (1) a noise amplitude and (2) a mix of white and flicker pure noise processes.

A number of specific kinds of changes are assumed that a clock might encounter. We assume that these changes can occur with either sign and with varying magnitudes. The

changes considered are a noise amplitude different from

nominal, a flicker component in the noise mix which is different from nominal, a step jump in the time counts for

a clock ind a linear (frequency jump) or quadratic (frequency

drift) trend in the time counts for a clock.

Detection of a change is accomplished with a multiple

sequential test having compound limits. The test is designed

to respond quickly to an actual change but to make few incorrect detections (identify the wrong change) or false detections

(when no change from nominal operation has occurred 1

. When a change is detected for a clock, the laboratory time

scale can be adjusted accordingly for this condition.

Key words: Atomic clock model; detection; flicker noise;

sequential test; time scale.

1. Introduction

A recent paper by Allan, Gray, and Machlan [l] gives a comprehensive view of the AT(NBS) time scale. These authors discuss the

influence of an individual commercial cesium beam clock on this time

scale. Of special interest for this report is a section of [l] on clock

modeling where the belief is expressed that a clock can experience
changes in its rate due to either internal or external perturbations. In

addition, the interpretation of "time" from a clock is subject to counting

noise. A prominent model of clock noise is a mix of two kinds of noise:

white FM and flicker FM noise processes. These noise processes are discussed in [2], [3], and [4]. In [l] the authors disclose that a weight w. is used to incorporate

measured counts from clock i into the AT(NBS) scale. They state that

this weight depends upon the quality of performance that clock i is expected to give. This means that the intensity and the kinds of noise

processes affect the weighting. The assumption of the quality of performance for clock i must be derived from a comparison of clock i with other clocks or from its direct calibration with a primary standard.

Under the assumption that the noise process of a clock might

change or that its rate might change, this research addresses the

problem of detecting changes in the performance of a clock. The noise

changes considered are intensity (amplitude) changes in either the 'white

FM or flicker FM component. The other changes are a jump in time, a

drift in time, or a drift in rate. The objective is to explore detection

schemes which have the following properties:

1. Make few false detections;

2. Ignore very minute changes;

3. Respond rapidly to large changes; and

4. Identify the type of change correctly.

Unfortunately these desired characteristics are in conflict in a number of ways. Very rapid detection of a change causes more false

detections. Trying to detect very small changes causes more false

detections. Trying to distinguish between many changes causes more

detections for the wrong reason. Therefore, a good detection scheme

is one which maintains a suitable balance between the goals of the

detection procedure. It is also one which responds quickly to actual
changes but makes only an acceptable number of false detections. A

design goal for the procedures developed in this project was to make a

false detection on the average only once per year with daily testing for

changes.

Detection schemes are used in many other applications. Two

very prominent examples are quality control and submarine warfare.

Other examples are early failure detection for operating machinery and

aircraft engines, medical diagnosis, and environmental control. Graphs

comparing the false detection frequency and the speed of detection of

various magnitudes of changes are often called operating characteristic curves. Another related notion is that machines are conceived to operate

in a particular manner by their designers. Often a failure of some type occurs causing abnormal performance from the machine. A detection

scheme can be devised to recognize non-nominal behavior when it is important to be aware of such conditions.

The statistical techniques used in this research are moving

averages, geometric smoothing, estimation, prediction, likelihood

functions, and comparisons of quantities against limits. The procedures

(which are really hypothesis tests) are multiple sequential tests using

compound limits to better distinguish between possible changes. Geometric smoothing is often employed to require a persistence of evidence

before a detection is made.

In section 2, we discuss the noise model of a clock that we have

used in this study. Section 3 shows the structure of the multiple detection

schemes which resulted from the experimental development using simulation of a clock. The bulk of the effort in this study was devoted to the

development of a multiple test which was consistent with the design

goals. Section 4 discusses this experimental development and then

discusses the results of an experimental evaluation of the procedures

of section 3. Section 5 gives conclusions and recommendations for
additional study. No attempt was made to specify any adjustments to be

made to a time scale after a change was detected. A change might

suggest a revision of the w. weights or a one time increment to be

applied to the scale.

2. Noise Model

The noise process about the rate of an atomic clock considered

in this study has a white noise FM component and a flicker noise FM

component. The white FM process has constant spectral density of the

frequency fluctuation while the flicker process has a spectral density

of the frequency fluctuation inversely proportional to the Fourier frequency. Let us consider other properties of these stochastic processes,

the autocovariance function and the Allan variance. The Allan variance

is specified as a sample statistic (time averaging of a process realization)

in [2], [3], and [4]. It will be convenient to view this measure of frequency

stability in its expectation form in this section.

Let W for t = 0, 1, . . . , • be a white noise stochastic process.

Let t be a day index; however, other counting intervals are equally valid.

Let us also denote a realization of the W process by w . We also assume

that random variable W is normally distributed with mean and unit

variance. The process outcomes at various days are independent, thus

the autocovariance function

AW

(T) = Cov(W

fc

, Wfc+T) = 0, for T > 0, and (1) A

w

(0) = Var(W

fc

) = 1.

Let <J (T) denote the Allan variance for the W noise process,

define

We

E[(W

t -Wt+1 ) 2

] <(1)
This expectation expands to E[W2

] - 2E[W

t Wfc+1

] + E[W

t+1

2

] 2

E[W

fc

2

] - 2E[W

t ]E[W

fc+1

] + E[W

t+i

2

] __

E[W

fc

3

] + E[W

fc+i

3

] =

2

'

since W and W are independent random variables. By definition we

have that W is distributed chi-square with 1 degree of freedom; the

expected value of a chi-square random variable with 1 degree of freedom

is 1. Thus,

Next consider

wt + wm

- wfc+2 + w

t+3 ^

((

v

S(2) =

_ 1 + 1 _ 2(0] 1 + 1 _ 1_ 4 " 4 4 ~ 2 '

In similar fashion it is easily shown that

<V3(T) = T • (3)

Let F be a flicker noise stochastic process. It is shown in [4] that

CT 2 (T) = c, a constant. (4) F

However, the autocovariance function

A (T) 4 for any T . (5) F

The flicker process is thought to be positively autocorrelated for all T;

this result should be demonstrable using eq (4) and the definition of

A (T) since many similar terms are present. It is necessary to divide
F by c so that

<T*(T) =1 (6)

in this model.

Let us define noise stochastic process

Z

t = a Wt + b F

t

, (7)

where a + b = 1. The Allan variance for this process is

o*(T) = 7 + b

2

. (8) Z T

For T= 1 in eq (8) we see that

2 , , E[(aW

t -aWm) 3

] E[ (bF^ - bF^ ) 2

]

cr z

a) = 5 + 1

= a

2 aw

2 (l) + b

2 aF2 (l).

Henceforth, let us denote the Allan variance of the Z process by simply

2 A 2 a (T). We will let a (T) be an estimate of the Allan variance from a particular process realization. Let us also define the notation

sigma (a, b, T) =yjj- + b

2

. (9)

Table 1 gives values of this function for a and b equal to 0. , .1, .2, . . .

,

1. and T + 1, 2, 4 and 8. Let us also define a particular estimate of the Allan variance for

the Z process with outcomes z by

T

/— — \2 £ ( vvi':

CTJ(T) = 1=1

, (10)

2(7- 1) where z". is the arithmetic average of the i-th group of T Z outcomes.

For comparison, a common estimator for the autocovariance function is
for T= 0, 1,

of Z..

T-T

AT

(r) = ^

E(zfc - S) ( Zfc+T

- i) , di:

, T-l and where z is the average over all T outcomes

TABLE 1. FUNCTION SIGMA (a, b, T).

. 1 .2

. 3 .4

1

. 5 .6

.7

.8

.9

1.

0. . 1 . 3

sigma (a, b, 1) b

.4 . 5 .7 1.0

0.000 . 100 .200 . 300 .400 . 500 .600 .700 . 800 . 900 1.000

. 100 . 141 .224 .316 .412 . 510 .608 .707 . 806 . 906 1.005

.200 .224 .283 .361 .447 . 539 .632 .728 . 825 .922 1.020

. 300 .316 .361 .424 . 500 .583 .671 .762 . 854 .949 1. 044

.400 .412 .447 . 500 . 566 .640 .721 . 806 . 894 .985 1.077

. 500 . 510 . 539 . 583 .640 .707 .781 .86 . 943 1.030 1. 118

.600 .608 .632 .671 .721 .781 .849 . 922 1.000 1.082 1.166

.700 .707 .728 .762 . 806 .860 . 922 .990 1. 063 1. 140 1.221

. 800 .806 . 82 5 . 854 .894 . 943 1. 000 1.063 1. 131 1.204 1.281

. 900 . 906 . 922 .949 . 985 1. 030 1.082 1. 140 1.204 1.273 1.345

1. 000 1.005 1.020 1. 044 1. 077 1. 118 1.166 1.221 1-281 1.345 1.414

0.

. 1 .2

.3

.4

1 . 5 .6

.7

. 8 .9

1.0

sigma (a, b, 2)

b

.4 . 5 .7 1.0

0. 000 . 100 .200 . 300 .400 .500 .600 .700 . 800 . 900 1.000

. 071 . 122 .212 . 308 .406 .505 .604 .704 . 803 .903 1.002

. 141 . 173 .245 . 332 .424 .520 .616 .714 . 812 .911 1.010

.212 .235 .292 . 367 .453 .543 .636 .731 . 828 .92 5 1.022

.283 . 300 .346 .412 .490 .574 .663 .755 . 849 .943 1.039

.354 .367 .406 .464 . 534 .612 .696 .7 84 . 875 .967 1.061

.424 .436 .469 . 52 . 583 .6 56 .735 .819 . 906 • 995 1.086

.495 .505 . 534 . 579 .636 .704 .77 8 . 857 . 941 1.027 1. 116

. 566 . 574 .600 .640 .693 .755 . 82 5 . 900 . 980 1.06 3 1. 149

.636 .644 .667 .704 .752 .809 . 87 5 . 946 1. 022 1. 102 1. 185

.707 .714 .735 .768 . 812 . 866 . 927 .995 1.068 1. 145 1.225
. 1

Table 1. (continued)

sigma (a, b, 4)

b

.3 .4 .5 .6 .9 1.

0. 0. 000 . 100 .200 . 300 .400 . 500 .600 .700 . 800 . 900 1. 000

. 1 . 050 . 112 .2 06 . 304 .403 . 502 .602 .702 . 802 . 901 1. 001

.2 . 100 .141 .224 .316 .412 . 510 .608 .707 . 806 . 906 1.005

. 3 . 150 . 180 .250 . 335 .427 . 522 .618 .716 . 814 . 912 1.011

.4 .200 .224 .283 .361 .447 . 539 .632 .728 . 825 . 922 1. 020

. 5 .250 .269 . 320 . 391 .472 . 559 .650 .743 . 838 . 934 1. 031

.6 .300 .316 .361 .424 . 500 . 583 .671 .762 . 854 .949 1. 044

.7 . 350 . 364 .403 .461 . 532 .610 .695 .783 . 873 . 966 1. 059

. 8 .400 .412 .447 . 500 . 566 .640 .721 . 806 . 894 . 985 1.077

.9 .450 .461 .492 . 541 .602 .673 .750 . 832 . 918 1. 006 1.097

1.0 . 500 . 510 .539 . 583 .640 .707 .781 . 860 . 943 1. 030 1. 118

sigma (a, b, 8)

0. . 1 .2 . 3 .4 . 5 .6 .7 . 8 .9 1.0

0. 0.000 . 100 .200 . 300 .400 . 500 .600 .700 . 800 . 900 1.000

. 1 . 035 . 106 .203 . 302 .402 . 501 .601 .701 . 801 . 901 1.001

.2 . 071 . 122 .212 . 308 .406 .505 .604 .704 . 803 . 903 1.002

. 3 . 106 . 146 .226 . 318 .414 .511 .609 .708 .807 . 906 1.006

.4 . 141 . 173 .245 . 332 .424 . 520 .616 .714 . 812 . 911 1.010

. 5 . 177 .203 .267 . 348 .437 .530 .625 .722 . 819 . 917 1.016

.6 .212 .235 .2 92 . 367 .453 . 543 .636 .731 . 828 . 92 5 1. 022

.7 .247 .267 . 318 . 389 .470 . 558 .649 .742 . 837 .933 1.030

.8 .283 . 300 .346 .412 .490 . 574 .663 .755 . 849 . 943 1.039

.9 . 318 . 334 .376 .437 . 511 . 593 .679 .769 . 861 .955 1.049

1.0 . 3 54 . 367 .406 .464 . 534 .612 .696 .784 . 87 5 . 967 1.061
3. Structure of the Detection Procedures

In this section we merely present the detection procedures developed in this research. The rationale behind them is defended in the next

section. Let T denote the time from a clock at day t; we update time by

Wi+ v (12)

Let us first consider the tests for a noise intensity change from

eq (7) of either of two types:

Z

t = a W + 8 F (13)

or Z

fc = aWt + bFt

, (14)

where Or and B are changed values of a and b, respectively.

An estimate of CT(T) is obtained by letting T = 16 in eq (10). Every

5 days thereafter let our estimate of 0"(T) be replaced by

.95 CT(T) + . 05 CT(T) . (15)

At day 16 we start the estimator with 0"(t) = O" , (T). The noise test is 16

performed every 5 days.

It proved necessary to adjust the 0" , (T) values to make them more

closely correspond to the theoretical values of the sigma (a, b, T) function.

The multiplicative factors at T = 1, 2, 4, 8 are . 97, 1.02, 1.05, 1.11,

respectively. The need for this correction stems from a minor imperfection in the flicker generator [5] and from a skewed sampling distribution for & (T). The sampling skew is suggested by analogy to the

1

following example. Let random variable X be distributed chi-square

with 1 degree of freedom. The E[X] = 1 and the Pr{X<l]is approximately

equal to .7. So we infer by analogy that in a small sample of X outcomes,

a number of low values is very probable. Since 0"

, (8) is estimated

16

from a single sample it often assumes a value lower than sigma (a, b, 8). Different correction factors would likely be better for T ^ 16.
A likelihood test is used to detect a noise change. Assume that

CT(1) has a greater value than nominal in a realization. This suggests

that a noise intensity increase may have occurred such that either

tt > a or 3 > b; we need to determine which one or neither (the test does

not explicitly look for a change in both).

If CT (1 ) > min{sigma (a + . 1, b, 1); sigma (a, b + . 1, 1)1, (16) a threshold value, let us start the following procedure. The first day

that CT(1) exceeds the threshold value let x = x = and x = 1. At other

times we update the x. values. We interpolate twice on (7(1) to obtain the

best fit to both sigma (a, b, 1) and sigma (a , b, 1). Then we compute

the absolute differences

| CT(T) - sigma (a, b' T) |

| CT(T) - sigma (a' b, T) |

| <7(T) - sigma (a, b, T ) | .

We update the x. by replacement signified as

*i

= - 75

*i + - 25

(

i 1 i )

' (17)

Vd! d2 d 3 /

In similar manner, if

CT(1 ) < max {sigma (a - . 1, b, 1); sigma (a,b-.l,l)} (18)

we start y = y = and y = 1 the first time the threshold is exceeded

and update on subsequent consecutive times. We interpolate on (T(l

again for both possible changes and compute d. as above and update the

y. by replacement signified as

V T= E

1, 2, 4, 8

d = 2 T= E

1, 2, 4, 8 V T= E

1, 2, 4, 8

t _l_

d.

i <

+ t

+

1 d

3

yi =.75y. + .25l-7--7--

T-J . (19)

10
3 3 Thus, the V x. = £ y. = 1 always.

i=l i=l

The test using compound limits is when x > . 55 and x < .20 detect P > b

when y > . 55 and y < . 20 detect & > a

A (20) when x > . 55 and x < .20 detect P<b

when y > . 55 and y < .20 detect tt < a . This test has 7 parameters. For review they are

1. T = 16 day averaging interval;

2. Smooth every 5 days into <T(T);

3. . 95 smoothing factor in eq (15);

4. The . 1 threshold in eqs (16) and (18);

5. .75 smoothing factor in eqs (17) and (19);

6. . 55 limit in eq (20); and

7. . 20 limit in eq (20) . Two procedures are used for detection of a time jump change.

This change is signified by replacing z with z plus or minus a step as shown as follows:

z

fc = zfc ± J . (21)

The first procedure has 3 parameters and either detects or fails

to detect jump J on day t. Consider the filter

pt =.90p

t _ 1 + .10ztl , (22)

where p = p = z . We let p predict the value z and if

|p

fc - zj > 3.4 + .6b (23)

we say that a jump must have been present at day t. The second procedure using 5 parameters is a backup to the first

since if zJ > 2. 8 + .6b (24)

and

11
(25)

and

|zk |<.8|z| for all k = t-10, ..., t-1, t+1, ..., t+10 (26)

we detect that a jump occurred at day t. A jump at day t can only be

detected at day t or at day t+10 using these procedures.

The drift change conditions can be represented by letting

z = z. ± D ± Q(t - t ), (27)

t t c where t is the day that the drift started,

c A predictor

r = .95 r + .05z (28)

is compared to a 5-day moving average

t£ z k

z c = k=t-4 (29)

5

_

in order to update one of four test quantities P , P , N , or N . The

initial values at day 1 are P^ = N = . 75 and P = N = . 50. The update

ri ri J—

<

J—

i

is signified in replacement form as follows:

if r > and z" > 0, then PTT = . 8PTT + . 2z"_

t 5 H H 5

or if r <; and z r < 0, then N = . 8N - . 2z"

5 ti ri b

C*,C\\ or if r > and z" ^ 0, then NT = . 8NT - . 2"z_

t 5 Li-ib

or if r £ and z\. > 0, then P_ = . 8P T + . 2z"

t 5 J_i J_/ b The drift test is if

(PH

- N

R

> . 5 and P

L

> N

L

)

or (PR- N

H

> .7) (3D

or (P_ > . 9NU and P__ > . 8)

Li ri ri

12
detect a drift in the positive direction, or if (P - N > . 5 and P > N ) or (P - N > .7) H H L L H H

(32) or (NT > . 9P„ and N > . 8) L ri H

detect a drift in negative direction. This test using 9 parameters is made every 5 days.

4. Experimental Development and Evaluation

4. 1. Development of the Detection Procedures

Considerable effort was required to verify and normalize the Z

generator so that cr J 1

) and 0" (1) averaged out to 1. The experimental

evaluation of the noise model was by simulation. A sequence of z values

was generated for t = 1,2,..., 512 days. At day t , usually around day

c

30, a change could be introduced of the type shown in eqs (13), (14),

(21), or (27). If no change were made at day t the run was said to be

nominal.

The CT(T) estimator was studied. The parameter choices made

were considered to be suitable. Longer averaging times have both

advantages and disadvantages. A major drawback is that the estimator

would require more than 16 days to start and would be less responsive

to changes in the a and b variables. In the noise detection scheme the

. 1 threshold parameter was chosen because it was near the outer edges

of the variation experienced in (7(1) under nominal conditions. It is possible to detect noise changes of less than ± . 1 variation around either

a or b.

In the case of the jump change the limit values in eqs (23) and (24)

were chosen to balance falsely detecting a jump and actually detecting ones just inside the noise level. The .6b term in these limits reflects

the apparently greater variation in the F process than in the W process.

13
The drift condition proved to be complex and the most difficult to

detect. The motivation behind the P , N , P , and N quantities

(meaning positive high, negative high, positive low, and negative low,

respectively) is as follows: The Z process is autocorrelated due to bF . The direction of the process is generally guided by the low frequency "energy" in the bF process, while the higher frequencies

generally come from the aW process. The drifts D or Q of eq (27) can be in a positive or negative direction. (If the signs are different

we must wait for Q(t - t ) to dominate. ) Let us assume that the r

c t predictor in eq (28) indicates a positive direction and that a positive

drift was introduced at day t . We contend that z" will generally be

c 5 greater in magnitude when it is in the same direction as r and generally

smaller in magnitude when it is running counter to r . The reason for

this condition is that D and Q(t - t ) are additive in one case and subc tractive in the counter case. Thus in this example explanation P

should be getting larger and N should be getting smaller. By reversing

J—; signs in the above example N would get larger and P smaller.

ri L

Several other schemes were tried using this same notion. They

proved less effective until a very large number of days beyond day t . For D only not zero in eq (27) this scheme is complicated by the predictor

gradually incorporating the D term into the prediction of z" . It should

5 also be noted that D does not affect <t (r); the Q term gradually will,

1

however.

4.2. Evaluation of the Detection Procedures

We estimate that the time between false detections is on the average

about 400 days when b = . 3. (The design center for these schemes was

a = . 7, b = . 3. ) Recall that 1 year was desired overall. Very accurate

estimates could not be made as many nominal runs went to 512 days

without a false detection. Figures 1 to 5 display cumulative noise for

5 nominal runs that went to 512 days without detection.

14
An evaluation of the jump change was made in 27 trials where | J | ranged from 3. to 4. 6. We experienced 21 jumps detected and 6 jumps

missed. Figures 6 to 9 display 4 of these 21 correctly detected jumps.

Of the 6 trials where the jump was missed, 4 went to 512 days nominally

(as they should have) while 1 drift and 1 Oi < a condition were mistakenly

detected. As |j| increases, the chance of successfully detecting the

jump also increases toward certainty.

In an evaluation of the time drift change consisting of 12 trials

with |D| equal to 1. or 1.5 and Q = 0, we obtained 9 correct detections,

2 missed detections which reported nominal to 512 days, and 1 P < b

mistake. The average number of days to make a |D| = 1.0 detection

was 83 while a |d| = 1.5 detection was reduced to 27 days. Figures

10 and 11 show a correct D = 1. and D = 1. 5 detection, respectively.

A |d| value around 3 might easily be mistaken for a jump. Values less

than 1 require on the average more than 83 days to detect while values

very near zero may not be detected much faster than the interval between

false detections.

The rate drift change combined with a small time drift was evaluated in 1 trials; this resulted in 8 correct detections, 1 jump by mistake

on the 6th day of drift and 1 (X < a condition on the 35th day of drift.

Had these 2 improper detections not been made so quickly, it is likely

that all 10 rate drifts would have been properly detected. In the trials

where |d| = . 1 and |Q| = .02 the average number of days until detection

was 76 while with |d| = .2 and |Q| = .05 the average was reduced to only 37 days. In practice, it is likely that Q values would be near zero

requiring a longer period of drift before proper detection is made.

Figures 12 and 13 display correct detections of combined time and rate

drifts.

15
An evaluation of the noise change consisted of 45 trials where

0! = a ± .2 and j3 = b ± .2. The overall results were that 19 correct

detections were made; in addition we had 9 cases where the direction

of the intensity was correctly sensed (a fi change was reported instead

of an a or vice versa), while 17 mistakes or missed detections occurred.

Of the 17 incorrect cases 11 were drifts, 4 were jumps, and 2 "were nominal to 512 days. In the correct detection cases the average number

of days observed to detect P = b + . 2, j3 = b - . 2, Of = a + . 2, 0! = a - . 2 were about 200, 175, 150, and 110 days, respectively. Figures 14 to 18

display 5 correct noise change detections.

The most successful detection by percentage was the (X = a - .2

followed by j3 = b - . 2, Ct= a + . 2, with /3 = b + .2 being the hardest

change to detect. It appears that the reason for |S = b - .2 being detected

slowly (175 days) was the difficulty in distinguishing this change from

the 0! = a - .2 change rather than from nominal. Most of the mistaken

drift detections occurred on the j3 = b + . 2 and a. = a + . 2 trials. It is

easily seen why these mistakes occurred in trials where the Z noise

process has larger than nominal outcomes as opposed to the cases where Z is lower than nominal. In the 9 cases where we detected the

wrong reason for a noise intensity increase or decrease but the correct

direction, it was because the likelihood of the change indicated was

greater than the other non-nominal alternative. In the cases where

Ci = a - .2 was detected instead of the proper /3 = b - .2 condition, the

detection interval averaged around 100 days. When more days were

taken to make the detection, more chance of getting the proper reason was shown; logically, this represents an effectively larger sample size

and thus should do better. Figures 19 to 23 display 5 noise changes

that were not correctly detected by procedure in eq (20).

16
An alternative way to check for a noise change might be to compute

<7 __(T) for T = 1, 2, 4, 8, 16, 32, 64, and perform an interpolation on

12 8

CT 00 (1) and ^ioo(^) - Then we could compute a likelihood of any change.

12 8 12 8 However, in this method one may go a number of 12 8 day intervals

before detecting a change since the cr (T) estimator experiences a

12 8 fairly substantial amount of variation about the theoretical sigma

(a, b, t) values.

+20

CD

.1-1 o

uou

-100

o o o O o

o o o O o

.—

(M ro * in Day

Figure 1. A nominal run with (a, b) = (. 7, . 3) and no detection in 512 days.

17
o

u

o

—IU

+ 150

+ 100

+ 50

Day

Figure 2. A nominal run with (a, b) = (. 7, . 3)

and no detection in 512 days.

CO

u

o—i U

+ 65

+ 50

+ 25

o o

f\i Day m

Figure 3. A nominal run with (a, b) = (. and no detection in 512 days.

.21

18
+45

+30

03

•1-4 o

53o

o

I—

U

+ 10

+2. 5

Figure 4. A nominal run with (a, b) = (. 8, . 2)

and no detection in 512 days.

+10

ID

•|Ho

53^-50

oU

-100

DayFigure 5. A nominal run with (a, b) = (. 8, . 2)

and no detection in 512 days.

19
+ 85

ox

u

o

»—

U

to

•iHo

u

o

r—

U

+ 50

Figure 6. Detection of a J = 3.4 jump by eq (23)

when (a, b) = (. 8, . 2) and t = 35.

c

Figure 7. Detection of a J = -3. 2 jump by eq (23)

when (a, b) = (. 8, . 2) and t = 34.

c 20
+ 55

+25

o

o

o

I—

I

U

Figure 8.

DayDetection of a J = when (a, b) = (. 7,

3. 5 jump by eq (23)

3) and t 25.

+20

v +10 eno

o

o

I—

I

u

Day

Figure 9. Detection of a J = 3.0 jump by eqs (24-26)

when (a, b) - (. 7, . 3) and t =35.

21
+ 55

+25

to

• Ho

oo

.—I U

DayFigure 10. Detection of a D = 1.0 time drift by

eq (31) when (a, b) = (. 8, . 2) and t =35.

Day

Figure 11. Detection of a D = 1.5 time drift by

eq (31) when (a, b) - (.7, . 3) and t = 31

22
Day

Figure 1 2 . Detection of a D = - . 1 , Q = -. 02 time and

rate drift by eq (32) when (a, b) - (. 7, . 3) and t c = 28.

+ 25

ox

u

oU

-25

Day

Figure 13. Detection of a D = . 2, Q = . 05 time and

rate drift by eq (31) when (a, b) = (. 7, . 3) and t = 320.

23
01

1-1O

uo

I—

U

-50

Day -1 -1 Figure 14. Detection at day 126 of noise change 8 = . 6 when (a, b) (.7, . 3) and t = 2! c

Figure 15. Detection at day 301 of noise change 8 = . 1 when (a, b) = (. 7, 24

3) and t = 26.

c
Figure 16. Detection at day 271 of noise change (8 when (a, b) = (. 7, . 3) and t = 28.

+ 10

o

u

oU

DayFigure 17. Detection at day 1 1 1 of noise change 0/ when (a, b) = (. 7, . 3) and t = 24.

25
+ 80

tn

• r-l o

u

o

I—

I

U

+ 50

+2 5

^ DayFigure 18. Detection at day 216 of noise change Oi = . 5 when (a, b) = (. 7, 3) and t = 28.

c

en

•r-l OX

uoU

Day

Figure 19. False inference at day 111 that (X> .7 when

actually jS = . 5. (a, b) = (. 7, 3) and t 26.

c 26
+ 30

o

u

oU

+20

+ 10

DayFigure 20. False inference at day 81 that a > .7 when

3) and t c actually j3 = . l . (a, b) = (. 7, = 31,

too

ooU

o

in Figure 21,

Day^

False detection of a jump by eq (23) at day 164

when actually j3 = . 5. (a, b) = (. 7, and t = Z< 27
a;

CO

•r-l OXuo5

DayFigure 22. False detection of ? drift by eq (31) at day 316

when actually = . 4. (a, b) = (. 7, . 3) and t = 21

+ 140

+ 100

0)

to H

o

u

o

I—

I

U +50

Day

Figure 23. False detection of a drift by eq (31) at day 276

when actually /S = . 4. (a, b) = 2) and t = 29.

c 28
5. Conclusions

We can conclude that the false detection interval came out as

desired. Also, nearly all of the runs where some change is introduced

a detection is obtained; however, sometimes it is an incorrect detection.

Rationale as to why mistakes are being made is available and enhances

our insight into the difficulty of establishing the cause of any apparently

non-nominal behavior of an individual clock in an ensemble. The

experimental evaluation reported in section 4 was not intended to be

comprehensive, but rather to give a general idea of successes and

failures of these detection procedures.

The 24 parameters used is an unfortunately large number, however,

the results indicate that even more sophisticated schemes might be

necessary to obtain better discrimination; this could require even more

parameters.

The problems of B > b in the noise process looking like drift and

3 < b appearing as though Of < a since 8 cannot be allowed to go below

zero when a lower than average (7(1 ) estimate is obtained can be resolved

by specializing the test to the b intensity to a greater degree. We can

also be comforted by the contention that jumps within the noise level

should not degrade a time scale very much.

It could be very useful to know more about the conditional distribution of Z given knowledge of Z, for all k ^ t as would other information

t+T & & k

about the Z noise process in terms of classical mathematical statistics.

It would also be useful to know more about the sampling distribution for

CTrpfT). In short, there are many more things to be learned about

detection of changes in atomic clock performance. However, it is hoped

that some of the concepts discussed in this report will aid the ongoing

effort of improving atomic time scales.

29
6. References

[l] Allan, D. W. , Gray, J. E., and MachIan, H. E. , "The National

Bureau of Standards atomic time scale system: Generation,

dissemination, precision, and accuracy, " IEEE Trans. Instrum.

andMeas., IM-21 , No. 4, pp. 388-391 (November 1972).

[2] Allan, D. W. , "Statistics of atomic frequency standards," Proc.

IEEE, 54, No. 2, pp. 221-230 (February 1966).

[3] Barnes, J. A., et al. , "Characterization of frequency stability,"

IEEE Trans. Instrum. andMeas., IM-20, No. 2, pp. 105-120

(May 1971).

[4] Barnes, J. A., and Allan, D. W. , "A statistical model of flicker

noise," Proc. IEEE, 54, No. 2, pp. 199-207 (February 1966).

[5] Barnes, J. A., and Jarvis, Stephen, Jr., "Efficient numerical

and analog modeling of flicker noise processes, " NBS Technical

Note 604 (June 1971).

30
FORM NBS-114A (1-71) 2. Gov't Accession No.

IT.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. NBS-TN-636

3. Recipient's Accession No.

4. TITLE AND SUBTITLE Modeling of Atomic Clock Performance and Detection

of Abnormal Clock Behavior

5. Publication Date March 1973

6. Performing Organization Code

7. AUTHOR(S)

William A. Ganter

8. Performing Organization

10. Project/Task/Work Unit No. 2730101

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS, Boulder Labs.

DEPARTMENT OF COMMERCE

Boulder, Colorado 80302

11. Contract/Grant No.

12. Sponsoring Organization Name and Address Same as 9

13. Type of Report & Period Covered Final

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant bibliography or literature sutvey, mention it here.) We have assumed that the nominal performance of an atomic clock can be well

characterized by (1) a noise amplitude and (2) a mix of white and flicker pure noise

processes. A number of specific kinds of changes are assumed that a clock might

encounter. We assume that these changes can occur with either sign and with

varying magnitudes. The changes considered are a noise amplitude different from

nominal, a flicker component in the noise mix which is different from nominal, a

step jump in the time counts for a clock and a linear (frequency jump) or quadratic

(frequency drift) trend in the time counts for a clock.

Detection of a change is accomplished with a multiple sequential test having

compound limits. The test is designed to respond quickly to an actual change but

to make few incorrect detections (identify the wrong change) or false detections

(when no change from nominal operation has occurred). When a change is detected

for a clock, the laboratory time scale can be adjusted accordingly for this condition.

17. KEY WORDS (Alphabetical order, separated by semicolons) Atomic clock model; detection; flicker noise; sequential test; time scale,

18. AVAILABILITY STATEMENT

HEl UNLIMITED.

I I FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE

TO NTIS.

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF PAGES

39

22. Price 50

USCOMM-DC 662<4-P71

NBS TECHNICAL PUBLICATIONS

PERIODICALS JOURNAL OF RESEARCH reports National Bureau of Standards research and development in physics, mathematics, and chemistry. Comprehensive

scientific papers give complete details of the work,

including laboratorv data, experimental procedures, and theoretical and mathematical analyses. Illustrated with photographs, drawings, and charts. Includes

listings of other NBS papers as issued.

Published in two sections, available separately:

• Physics and Chemistry (Section A)

Papers of interest primarily to scientists working in these fields. This section covers a broad range of physical and chemical research, with major emphasis on standards of physical measurement, fundamental

constants, and properties of matter. Issued six times a year. Annual subscription: Domestic, $17.00; Foreign, $21.25.

• Mathematical Sciences (Section B)

Studies and compilations designed mainly for the mathematician and theoretical physicist. Topics in mathematical statistics, theory of experiment design, numerical analysis, theoretical physics and chemistry,

logical design and programming of computers and

computer systems. Short numerical tables. Issued quarterly. Annual subscription: Domestic, $9.00; Foreign,

$11.25.

TECHNICAL NEWS BULLETIN

The best single source of information concerning the Bureau's measurement, research, developmental, co- operative, and publication activities, this monthly

publication is designed for the industry-oriented

individual whose daily work involves intimate contact with science and technology

for engineers, chemists,

physicists, research managers, product-development

managers, and company executives. Includes listing of

all NBS papers as issued. Annual subscription: Domestic, $6.50; Foreign, $8.25.

NONPERIODICALS Applied Mathematics Series. Mathematical tables, manuals, and studies. Building Science Series. Research results, test methods, and performance criteria of building materials, components, systems, and structures. Handbooks. Recommended codes of engineering and industrial practice (including safety codes) de- veloped in cooperation with interested industries, professional organizations, and regulatory bodies. Special Publications. Proceedings of NBS conferences, bibliographies, annual reports, wall charts, pamphlets, etc. Monographs. Major contributions to the technical

literature on various subjects related to the Bureau's

scientific and technical activities. National Standard Reference Data Series. NSRDS provides quantitative data on the physical and chemical properties of materials, compiled from

the world's literature and critically evaluated. Product Standards. Provide requirements for sizes, types, quality, and methods for testing various industrial products. , These standards are developed co- operatively with interested Government and industry groups and provide the basis for common understanding of product characteristics for both buyers and

sellers. Their use is voluntary.

Technical Notes. This series consists of communications and reports (covering both other-agency and

NBS-sponsored work) of limited or transitory interest. Federal Information Processing Standards

Publications. This series is the official publication within the Federal Government for information on

standards adopted and promulgated under the Public Law 89—306. and Bureau of the Budget Circular A—86

entitled, Standardization of Data Elements and Codes

in Data Systems. Consumer Information Series. Practical information, based on NBS research and experience, cover- ing areas of interest to the consumer. Easily understandable language and illustrations provide useful background knowledge for shopping in today's tech- nological marketplace.

BIBLIOGRAPHIC SUBSCRIPTION SERVICES

The following current-awareness and literature-survey bibliographies are issued periodically by the Bureau:

Cryogenic Data Center Current Awareness Service (Publications and Reports of Interest in Cryogenics). A literature survey issued weekly. Annual subscription: Domestic, $20.00; foreign, $25.00.

Liquefied Natural Gas. A literature survey issued quarterly. Annual subscription: $20.00. Superconducting Devices and Materials. A literature survey issued quarterly. Annual subscription: $20.00. Send subscription orders and remittances for the preceding bibliographic services to the U.S. Department

of Commerce, National Technical Information Service, Springfield, Va. 22151. Electromagnetic Metrology Current Awareness Service (Abstracts of Selected Articles on Measurement

Techniques and Standards of Electromagnetic Quantities from D-C to Millimeter-Wave Frequencies). Issued monthly. Annual subscription: $100.00 (Special rates for multi-subscriptions). Send subscription order and

remittance to the Electromagnetic Metrology Information Center, Electromagnetics Division, National Bureau

of Standards, Boulder, Colo. 80302. Order NBS publications (except Bibliographic Subscription Services) from : Superintendent of Documents. Government Printing Office, Washington, D.C. 20402.
U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington. D.C. 20234

OFFICIAL BUSINESS

Penalty for Private Use. $300

POSTAGE AND FEES PAID

U.S. DEPARTMENT OF COMMERCE