An Unbiased FIR Filter for TIE Model of a Local Clock in ... - IEEE Xplore

4 downloads 0 Views 523KB Size Report
K-degree polynomial model of a local clock in global posi- tioning system (GPS)-based timekeeping in the presence of noise that is not obligatory Gaussian.
862

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 53, no. 5, may 2006

An Unbiased FIR Filter for TIE Model of a Local Clock in Applications to GPS-Based Timekeeping Yuriy S. Shmaliy, Senior Member, IEEE Abstract—An unbiased finite impulse response (FIR) filter is proposed to estimate the time-interval error (TIE) -degree polynomial model of a local clock in global positioning system (GPS)-based timekeeping in the presence of noise that is not obligatory Gaussian. Generic coefficients for the unbiased FIRs are derived. The low-degree FIRs and noise power gains are given. An estimation algorithm is proposed and examined for the TIE measurements of a crystal clock in the presence of the uniformly distributed sawtooth noise induced by the multichannel GPS timing receiver. Based upon this algorithm, we show that the unbiased FIR estimates are consistent with the reference (rubidium) measurements and fit them better than the standard Kalman filter.

K

I. Introduction ast and accurate estimation and adjustment of a local clock performance, making possible for a variety of modern digital systems to operate in common time with minimum “slips”, is of importance for the global positioning system (GPS)-based timekeeping. Here, the time interval error (TIE) between the GPS time and the local clock time is measured in the presence of noise induced by the GPS timing receiver. The TIE then is estimated, and the correction signal adjusts the clock for the GPS time. The standard deviation of the noise using commercially available receivers is about 30 ns, can reach 10–20 ns [1], and may be improved by removal of systematic errors to no less than 3–5 ns [1], [2]. Having such a large noise, the measured data is usually not appropriate for clock correction, and a TIE tracking filter with a large time constant is applied. To obtain filtering in an optimum way, a TIE model of a local clock must be known for the filter memory. Such a model was proposed in [3] as the second degree Taylor polynomial and is now basic in timekeeping, being practically proven. In the discrete time, it may be written as:

F

x1 (n) = x1 (0) + x2 (0)τ n +

x3 (0) 2 2 τ n + w1 (n, τ ), 2

(1)

Manuscript received April 26, 2005; accepted May 3, 2005. This material was discussed December 7–9, 2004 at the Precise Time and Time Interval (PTTI) Meeting, Washington, DC. The work was supported by the CONACyT Project, J38818-A. Y. S. Shmaliy is with the Guanajuato University, Salamanca, Gto., Mexico (e-mail: [email protected]). He also is with the Kharkiv National University of Radio Electronics, Kharkiv, Ukraine.

Fig. 1. Typical TIE 1 PPS sawtooth noise induced by a GPS timing sensor SynPaQ III.

where n = 0, 1, . . . ; τ = tn − tn−1 is a time step multiple to the 1 s; tn is a discrete time; x1 (0) is an initial time error; x2 (0) is an initial fractional frequency offset of a local clock from the reference frequency; x3 (0) is an initial linear fractional frequency drift rate; and w1 (n, τ ) is a random component caused by the oscillator noise and environment. In GPS-based measurements, (1) is observed via the mixture: ξ1 (n) = x1 (n) + v1 (n),

(2)

in which v1 (n) is a noisy component induced at the receiver (noise of a measurement set is usually small). In modern receivers, such as the Motorola M12+ (see [4]), and SynPaQ III GPS Sensor (Synergy Systems, LLC, San Diego, CA), a random variable v1 (n) is uniformly distributed owing to the sawtooth noise caused by a principle of the 1 PPS (one pulse per second) signal formation. Fig. 1 shows a typical structure of v1 (n) in a short time; that is the modulo 2vmax Brownian TIE associated with a phase of a receiver local oscillator, where vmax is the noise bound (in SynPaQ III, it is about 50 ns. If the sawtooth correction is used, the noise is reduced by the factor of about 5 and becomes near Gaussian). In long-term measurements, v1 (n) exhibits nonstationary excursions due to uncertainties in the GPS time [5] with different satellite sets in a view. To estimate optimally the states of different clocks (atomic and crystals), several filtering algorithms have been examined during a couple of decades. For the state space equivalent of (1), the problem was formulated by Allan and Barnes in [6] to apply Kalman filtering. The solutions then were given in [7], [8]. Thereafter, the Kalman algorithm was used by many authors [9]–[13], and its applications in time scales were recently outlined in [14]. The problem with the Kalman filter arises when noise is not white and the model demonstrates temporary uncertainties; thus, the estimate may become biased and noisy. To overcome, advanced Kalman algorithms were proposed in [15], [16] being, however, not yet adopted for GPS applications.

c 2006 IEEE 0885–3010/$20.00 

shmaliy: proposed finite impulse response filter

863

An alternative approach is known as finite impulse response (FIR) filtering allowing for noise of arbitrary distribution. In contrast to the infinite impulse response (IIR) structures (including Kalman filters), FIR structures have inherent properties, such as a bounded input/bounded output (BIBO) stability and robustness against temporary model uncertainties and round-off errors [17]. They may be used independently or combined with Kalman filters [12], [18]. A general optimal FIR filtering algorithm with embedded unbiasedness for state space models was recently proposed in [19]. Especially for GPS-based timekeeping and a linear TIE polynomial model, an unbiased FIR filter was designed and studied in [20]. In this paper, we propose a new unbiased FIR filter and algorithm intended for real time estimating the TIE Kdegree polynomial model in the presence of a GPS noise of arbitrary distribution (with or without the sawtooth correction). The rest of the paper is organized as follows. In Section II, we formulate the problem. In Section III, a design of an unbiased FIR filter is given (the necessary coefficients are postponed in Appendix I). The low-degree FIRs and noise-power gains are derived in Section IV. An unbiased FIR-filtering algorithm for a single multichannel GPS timing receiver is described in Section V and applied for a local crystal clock without the sawtooth correction. The FIR estimates are compared here to the reference measurement (rubidium) and to those obtained using the 3-state Kalman filter (Appendix II). Concluding remarks are drawn in Section VI.

The most common situation that may be assumed in timekeeping is when all or several clock states are observable by (5) via M (independent or dependent) measurements in the presence of correlated or uncorrelated noises. Then the observation vector is ξ(n) = [ξ1 (n)ξ2 (n) . . . ξM (n)]T and: ⎡ c11 ⎢0 C=⎢ ⎣. . . 0

0 c22 ... 0

... 0 ... 0 ... ... . . . cMM

⎤ ... 0 . . . 0⎥ ⎥, . . . 0⎦ ... 0

(7)

is a measurement matrix [M × (K + 1)] with, typically, unit components cuu , u ∈ [1, M ]. The clock noise vector is w(n, τ ) = [w1 (n, τ )w2 (n, τ ) . . . wK+1 (n, τ )]T with the components caused by the oscillator noises. In view of Fig. 1, the noise vector v(n) = [v1 (n)v2 (n) . . . vM (n)]T contains correlated or uncorrelated components that are not obligatory Gaussian. It is important that the GPS noise v(n) dominates on a horizon N ; that is, typically, wu2 (n, τ )N  vl2 (n)N . Therefore, w(n, τ ) may be neglected in the averaging FIR procedure (it cannot be discarded in the Kalman filter). The problem now formulates as follows. Given a noisefree TIE model (4), w(n, τ ) = 0, that is observable via (5) in the presence of a mean zero noise, v(n) = 0, distributed arbitrary with a known covariance v(n)vT (n). We want to derive a real-time, unbiased FIR estimator of the clock states λ(n), supposing in this paper that M = K + 1, and C is an identity matrix.

II. Problem Formulation Most commonly, the TIE polynomial model projects ahead on a horizon of N points from the start point n = 0 with the K-degree Taylor polynomial: x1 (n) =

K 

p p

xp+1

p=0

= x1 + x2 τ n +

τ n + w1 (n, τ ) p! x3 2 2 x4 3 3 τ n + τ n · · · + w1 (n, τ ), 2 6

(3)

where xl+1 ≡ xl+1 (0), l ∈ [0, K], are initial states of the clock and w1 (n, τ ) is a noise with known properties. By extending the time derivatives of the TIE model to the Taylor series, (3) and (2) become, respectively: λ(n) = A(n)λ(0) + w(n, τ ), ξ(n) = Cλ(n) + v(n),

(4) (5)

where λ(n) = [x1 (n)x2 (n) . . . xK+1 (n)]T is a vector [(K + 1) × 1] of the clock states and a time-varying clock matrix [(K + 1) × (K + 1)] is: ⎡ ⎤ 1 τ n τ 2 n2 /2 . . . (τ n)K /K! ⎢0 1 τ n . . . (τ n)K−1 /(K − 1)!⎥ ⎢ ⎥ ⎢ . . . (τ n)K−2 /(K − 2)!⎥ A(n) = ⎢0 0 1 ⎥. ⎣. . . . . . . . . ⎦ (6) ... ... 0 0 0 ... 1

III. An Unbiased FIR Filter Using N points of the nearest past, the FIR estimate ˆ λ(n) = [ˆ x1 (n)ˆ x2 (n) . . . x ˆK+1 (n)]T at n-th point is given by, neglecting w(n, τ ): ˆ λ(n) =

N −1 

H(i)ξ(n − i)

i=0

(8a)

= q(n)Γ = d(n)Γ + r(n)Γ, where an auxiliary unit matrix (N ×1) is Γ = [1 1 . . . 1]T and the matrix [(K + 1) × (K + 1)] of unknown FIRs is: ⎡ ⎤ hK (i) 0 ... 0 ⎢0 ⎥ hK−1 (i) . . . 0 ⎥, H(i) = ⎢ ⎣. . . ... ... ... ⎦ 0 0 . . . h0 (i)

(9)

in  which the l-th FIR has inherent properties: hl (i) = N −1 hl (i), 0 ≤ i ≤ N − 1 and i=0 hl (i) = 1. The es0, otherwise

864

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 53, no. 5, may 2006

timates vector and its deterministic and random constituents, all of [(K + 1) × N ], are, respectively: q(n) = [H(0)ξ(n) H(1)ξ(n − 1) . . . H(N − 1)ξ(n − N + 1)],

(10)

d(n) = [H(0)Cλ(n) H(1)Cλ(n − 1) . . . H(N − 1)Cλ(n − N + 1)],

(11)

r(n) = [H(0)v(n) H(1)v(n − 1) . . . H(N − 1)v(n − N + 1)].

(12)

ˆ The mean square error (MSE) of λ(n) is written as: T ˆ ˆ [λ(n) − λ(n)] J(n) = [λ(n) − λ(n)]

= [λ(n) − d(n)Γ]T [λ(n) − d(n)Γ]

(13)

in which G = [1 0 . . . 0]T is of ⎡ 1 1 1 ... ⎢ 0 1 2 ... ⎢ 2 F=⎢ ⎢ 0 1 2 ... ⎣. . . . . . . . . . . . 0 1 2l . . .

[(l + 1) × 1] and ⎤ 1 N −1 ⎥ ⎥ (N − 1)2 ⎥ ⎥. ⎦ ... l (N − 1)

(22)

We notice that (21) is inherent for any other unbiased estimators, e.g., for the minimum variance unbiased (MVU) and best linear unbiased estimator (BLUE). It follows from Kalman filtering that an optimum estimate is achieved if the model degree is equal to a number of the filter states minus one. This means that the components in (18) also may be substituted by the l-degree polynomial, such as:

+ [r(n)Γ] [r(n)Γ], T

producing the estimate bias and variance, respectively:

hl (i) =

l 

ajl ij ,

(23)

j=0

ˆ bias[λ(n)] = λ(n) − d(n)Γ, ˆ var[λ(n)] = [r(n)Γ]T [r(n)Γ].

(14) (15)

where ajl are still unknown coefficients. Embedded (23), the constraint (21) becomes:

A. Generic Coefficients for the FIR of an Unbiased Filter

DB = G,

The necessary and sufficient condition for the unbiased estimate follows straightforwardly from (14); that is:

where the FIR coefficients matrix is B = [a0l a1l . . . all ]T , and an auxiliary quadratic matrix [(l + 1) × (l + 1)] is given by: ⎡ ⎤ d0 d1 d2 . . . dl ⎢ d1 d2 d3 . . . dl+1 ⎥ ⎥ ⎢ ⎥ D=⎢ (25) ⎢ d2 d3 d4 . . . dl+2 ⎥ , ⎣. . . . . . . . . . . . . . . ⎦ dl dl+1 dl+2 . . . d2l

λ(n) = d(n)Γ,

(16)

providing the rule to derive the FIRs for the clock states: ⎡ ⎤ ⎡ T ⎤ x1 (n) WK λ1 (n) ⎢x2 (n) ⎥ ⎢WT λ2 (n) ⎥ ⎢ ⎥ = ⎢ K−1 ⎥, (17) ⎣. . . ⎦ ⎣. . . ⎦ xK+1 (n) W0T λK+1 (n) where Wl = [hl (0)hl (1) . . . hl (N − 1)]T , ⎡ ⎤ xK+1−l (n) ⎢xK+1−l (n − 1) ⎥ ⎥. λK+1−l (n) = ⎢ ⎣. . . ⎦ xK+1−l (n − N + 1)

(18)

For the clock (K + 1 − l)-th state, (17) thus yields a relation: xK+1−l (n) =

WlT λK+1−l (n).

x ˆK+1−l (n) = =

(21)

M(j+1)1 , |D|

(26)

in which |D| and M(j+1)1 are the determinant and minor of (25), respectively. Determined ajl and hl (n), the unbiased FIR estimate of the clock (K + 1 − l)-th state becomes1 :

(20)

It is seen that (17) is valid for arbitrary n. Setting n = 0, expressing the components in (19) with the ldegree polynomial, by (3), and involving the property N −1 i=0 hl (i) = 1, we go from (17) to the unbiasedness (or deadbeat) constraint: FWl = G,

N −1 where the generic component dm = i=0 im , m ∈ [0, 2l], is determined by the Bernoulli polynomials (Appendix I). An analytic solution of (24) yields the generic coefficients for (23): ajl = (−1)j

(19)

(24)

=

N −1 

hl (i)ξK+1−l (n − i), i=0 WlT ΞK+1−l (n), GT (FFT )−1 FΞK+1−l (n),

(27a) (27b) (27c)

1 The second reviewer mentioned that (23) can be written as W = l FT B; therefore, (24) follows from (21) with D = FFT . We notice that (24) then may be solved by B = (FFT )−1 G and (27b) rewritten as (27c).

shmaliy: proposed finite impulse response filter

where:

⎡ ⎤ ξK+1−l (n) ⎢ξK+1−l (n − 1) ⎥ ⎥. ΞK+1−l (n) = ⎢ ⎣. . . ⎦ ξK+1−l (n − N + 1)

865

(28)

B. Estimate Noise The estimate noise variance (15) now may be rewritten as: ˆ var[λ(n)] = [r(n)Γ]T [r(n)Γ] =

K 

T WK−p Rp+1 (n)WK−p ,

(29)

p=0

where the autocorrelation matrix Rl (n) of (N × N ) has a generic component Rl (i, j) = vl (i)vl (j), i, j ∈ [n, n − N + 1]. Accordingly, the estimate variance associated with (27b) calculates: 2 σK+1−l (n) = WlT RK+1−l (n)Wl .

Fig. 2. Unique FIRs of the unbiased filter: h0 (i), by (32); h1 (i), by (33); h2 (i), by (34); and h3 (i), by (35). An example is given for N = 60.

h2 (i) =

(30) h3 (i) =

It is important that, by large N , the sawtooth noise becomes delta-correlated2. This degenerates Rl (n) to the 2 diagonal form with the components Rl (i, i) = σvl (i) and tends (30) toward: 2 (n) σK+1−l

=

2 σv(K+1−l) (n)WlT Wl ,

(31)

2 where σv(K+1−l) (n) is a variance of the noise perturbing the (K + 1 − l)-th clock state at nth point.

IV. Applications to the Clock TIE Polynomial Model In applications, K is identified for the filter memory on a horizon [0, N − 1]) by the clock precision. Typically, it is assumed that (3) fits atomic clocks, by K ≤ 1, and crystal clocks, by K ≤ 2. However, K = 3 may be required for low-precision crystal clocks. Below we derive the unbiased FIRs for all these cases.

Setting l = 0, 1, 2, 3 and using the coefficient dm (Appendix I), we first calculate (26). Then (23) leads to the unique FIRs, namely: 1 , N

2(2N − 1) − 6i h1 (i) = , N (N + 1)

(34)

8(2N 3 − 3N 2 + 7N − 3) − 20(6N 2 − 6N + 5)i N (N + 1)(N + 2)(N + 3) 120(2N − 1)i2 − 140i3 (35) + . N (N + 1)(N + 2)(N + 3)

Fig. 2 sketches (32)–(35) for N = 60. A uniform FIR (32) corresponds to simple averaging and is optimal in a sense of a minimum-produced noise. This FIR is practically proven to be reasonable in GPS-based common view measurements [5]. A linear FIR (33) was derived in [20] by using linear regression to compensate a bias of simple averaging. Its kernel starts with a maximum

2(2N −1)

∼ h2 (0) = N (N +1)

= N4 > 0 and goes to a minimum N 1

−2)

∼ h2 (N − 1) = − N2(N = − N2 < 0, having zero at (N +1)

n0 =

2N −1 3 .

N 1

h2 (0) N →∞ h2 (N −1)

Its special peculiarity is r = lim

=

−2 that allows one to synthesize a FIR, by saving r = −2 for arbitrary N . It is surprising that the FIR synthesized in such a way is equal to that derived in [21] for the 1-step linear prediction on a horizon [1, N ]. B. Noise-Power Gains

A. Low-Order FIRs for Unbiased Filters

h0 (i) =

3(3N 2 − 3N + 2) − 18(2N − 1)i + 30i2 , N (N + 1)(N + 2)

(32)

(33)

2 The sawtooth noise produced by SynPaQ III becomes practically delta-correlated if N = 1800 with τ = 1 s or N = 180 with τ = 10 s that corresponds to 0.5 hour of averaging.

The noise-power gain corresponding to the l-degree FIR 2 2 is specified, by (31), to be gl (N ) = σK+1−l /σv(K+1−l) = T Wl Wl . Its values associated with (32)–(35) are given below, respectively: 1 , N 2(2N − 1) g1 (N ) = , N (N + 1) 3(3N 2 − 3N + 2) g2 (N ) = , N (N + 1)(N + 2) 8(2N 3 − 3N 2 + 7N − 3) g3 (N ) = . N (N + 1)(N + 2)(N + 3) g0 (N ) =

(36) (37) (38) (39)

866

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 53, no. 5, may 2006

Fig. 4. Structure of the (K +1)-state unbiased FIR filtering algorithm for the K-degree TIE polynomial model observable with a single GPS timing receiver.

Fig. 3. Noise gains of the unbiased FIR filters: l = 0, by simple averaging (36); l = 1, by (37); k = 2, by (38); and k = 3, by (39). Dashed lines are the upper bounds calculated by (40).

Fig. 3 illustrates (36)–(39) manifesting that unbiasedness is achieved at increase of noise. Indeed, the curves for √ l > 0 trace above the lower bound 1/ N associated with simple averaging (l = 0) that produces minimum noise (among all filters). It also follows that, by large√N , the noise gain is performed by gl (N ) ∼ = (l + 1)/ N and thus traces below the upper bound:  √ (l + 1)/ N , N ≥ (l + 1)2 gl (N ) ≤ . (40) 1, N < (l + 1)2

V. An Unbiased FIR Filtering Algorithm for a Single, Multichannel GPS-Timing Receiver We now consider an important practical case in which the measurement ξ1 (n) of a TIE x1 (n) is obtained with a single, multichannel GPS-timing receiver. Here, observations of the higher-order states may be formed by increments of the lower-order estimates. An unbiased FIR filtering algorithm then is written as: ˆ λ(n) = q(n)Γ, ξk (n) =

1 [ˆ xk−1 (n) − x ˆk−1 (n − 1)], k > 1, τ

(41) (42)

where q(n) is given by (10), and the observation components for (10), k > 1, are formed by (42). The algorithm is illustrated in Fig. 4. The clock first state estimate x ˆ1 (n) is obtained with hK (i) at a horizon of NK points. The observation ξ2 (n) for the second state x2 (n) then is formed, using (42), by increments of x ˆ1 (n). Accordingly, x ˆ2 (n) is achieved with hK−1 (i) at a horizon of NK−1 points. Inherently, the first accurate value of x ˆ2 (n) appears at (NK + NK−1)th point starting from n = 0. The last state estimate x ˆK+1 (n) is calculated with h0 (i) at a

horizon of N0 points, using ξK+1 (n) that is formed in the same manner as ξ2 (n). The first correct value of x ˆK+1 (n) appears at (NK + NK−1 + · · · + N0 )th point. For the quadratic TIE model (K = 2, crystal clocks), the 3-state unbiased FIR batch algorithm becomes, by (27a) and (42): x ˆ1 (n) =

N 2 −1

h2 (i)ξ1 (n − i),

(43)

i=0

x ˆ2 (n) =

N1 −1 1  h1 (j)[ˆ x1 (n − j) − x ˆ1 (n − j − 1)], τ j=0

(44)

N0 −1 1  [ˆ x2 (n − r) − x ˆ2 (n − r − 1)], τ N0 r=0

(45)

x ˆ3 (n) =

where h2 (i) and h1 (i) are given by (34) for N = N2 and (33) for N = N1 , respectively. For the linear TIE model (K = 1, atomic clocks), (41) and (42) simplify to the 2state form of: x ˆ1 (n) =

N 1 −1

h1 (i)ξ1 (n − i),

(46)

i=0

x ˆ2 (n) =

N0 −1 1  [ˆ x1 (n − j) − x ˆ1 (n − j − 1)], τ N0 j=0

(47)

where h1 (i) is given by (33) with N = N1 . Each state also may be calculated using the matrix forms (27b) or (27c). Below, as an example of application, we use this algorithm to estimate the TIE model of an oven crystal clock embedded to the Stanford Frequency Counter SR620 (Stamford Research Systems, Inc., Sunnyvale, CA). The measurement is done with SynPaQ III and SR620 for τ = 1 s (GPS measurement). Simultaneously, to get a reference trend, the TIE of the same crystal clock is measured, by SR625 (Stamford Research Systems, Inc.), for the rubidium clock (Rb-measurement). Initial time and frequency shifts between two measurements then are eliminated statistically, and a transition to τ = 10 s is provided by the data thinning in time. At the early stage, the TIE model was identified to be quadratic, K = 2, and Nl are

shmaliy: proposed finite impulse response filter

867

TABLE I Average Error (error) and Allan Deviation (σ) of the Estimate (Est) for 9.7 Hours and τ = 10 s: F is FIR and K is Kalman. Errors are Given for Rb-Measurements.

Est F K K-F

x, ns error σx (10) 2.8313 3.1295 0.2977

1.3786 1.3627

y, 10−12 error σy (10) 1.4852 2.5698 1.0846

0.6399 0.7025

D, 10−16 /s error σD (10) 4.1660 5.3121 1.1461

1.0206 1.3881

determined for each estimate in the minimum MSE sense3 . We also compare the unbiased FIR estimates to those obtained with the 3-state Kalman filter (Appendix II). (a)

A. Several Hours Measurements In this experiment, a short-term measurement of the TIE has been done during several hours [Fig. 5(a)]. The algorithm then was run. The horizons were identified for τ = 10 s to be N1 = 155 or 0.43 hours, N2 = 950 or 2.64 hours, and N3 = 860 or 2.39 hours for the Rbmeasurements. Thereafter, we set the values of q’s in the Kalman filter (Appendix II) to obtain the minimum MSEs for the FIR estimates. Fig. 5 and Table I illustrate these studies, showing that the unbiased FIR estimates, x ˆ1 (n), x ˆ2 (n), and x ˆ3 (n), and the relevant Kalman estimates, x ˆ(n), yˆ(n), and zˆ(n), respectively, are consistent with, however, some differences. It follows from Table I that the FIR filter works accurately. Fig. 5(a) shows that x ˆ1 (n) and x ˆ(n) track the mean value of the GPS measurement and that their offsets from the Rb-measurement are coursed mostly by the GPS time uncertainty. In this experiment, a maximum estimate error of about 60 ns was indicated between the 8th and 9th hours when a time shift in the 1 PPS signal has occurred. It follows [Fig. 5(b)] that x ˆ2 (n) and yˆ(n) fit well the weighted by 1/τ increments of the Rb-measurement. Even so, there are two special ranges (dashed). In range I, the frequency shift of about 3 × 10−11 has occurred in the span between the 7th and 8th hours, and no appreciable error is indicated in a range of large time shifts (between the 8th and 9th hours) in Fig. 5(a). We associate it with the frequency shift in SR625. In range II, the Kalman filter demonstrates a brightly pronounced instability caused likely by the temporary model uncertainty, but the FIR estimate still is consistent. We watch for a bit shifted trends of x ˆ3 (n) and zˆ(n) in Fig. 5(c) that may be explained by some inconsistency between the q’s and Nl . It also is seen that zˆ(n) traces much upper x ˆ3 (n) after about 8.7 hours. We associate it with the Kalman filter instability, like the case of a range II in Fig. 5(a). 3 To identify K, x ˆ1 (n) is compared to the reference (rubidium) measurement in the MSE sense by changing N for K ∈ [0, 3]. The minimum MSE identifies K and determines N for the TIE estimate. Given K, the other Nl are determined in the same manner.

(b)

(c) Fig. 5. Short-time measurement and estimation of the crystal clock TIE model with the 3-state unbiased FIR algorithm and the 3-state Kalman filter: (a) TIE, (b) fractional frequency offset, and (c) linear fractional frequency drift rate.

The experiment was repeated for τ = 1 s. The results are presented in Table II to mention that, on the whole, the picture (Fig. 5) remains the same. The only principle point to notice is that the Allan deviations of all estimates are reduced by a large number of the points. The FIR and Kalman estimates behave here closer to each other, even though the former is still more accurate with its lower error and much lower Allan variance. B. Long-Term Measurements The same crystal clock was later examined during about 2.5 days using only the unbiased FIR filter. The measure-

868

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 53, no. 5, may 2006

TABLE II Average Error (error) and Allan Deviation (σ) of the Estimate (Est) for 9.7 Hours and τ = 1 s: F is FIR and K is Kalman. Errors are Given for Rb-Measurements.

Est F K K-F

x, ns error σx (1) 2.8127 2.8965 0.0838

0.1374 0.3956

y, 10−12 error σy (1) 1.5638 2.4924 0.9286

0.0641 0.2131

D, 10−16 /s error σD (1) 19.147 20.346 1.1990

0.1037 0.4315

ments inherently show oscillations caused by day’s variations in temperature [Fig. 6(a)] and, like the previous case, all FIR estimates fit well the Rb-measurement. Using x ˆ2 (n), the temperature drift [Fig. 6(b)] was estimated to be about 2 × 10−10 (14 to 24◦ C) and xˆ3 (n) calculates the aging rate by ˆ x3 (n) = 0.4 × 10−10 /day [Fig. 6(c)].

(a)

VI. Conclusions In this paper we presented an unbiased FIR filter for the TIE K-degree polynomial model of a local clock. In contrast to the standard Kalman filter, the proposed solution does not require the noises to be white and does not involve any knowledge about noises in the algorithm. Instead, the FIR filter needs a length Nl that is determined for a given clock using a reference source in the minimum MSE sense. The filter produces a noise with a variance that reduces as a reciprocal of Nl . We notice that timekeeping operates with large horizons, N 1, and thus one should not expect appreciable discrepancy between the optimum and unbiased estimates. The trade-off between the 3-state unbiased FIR algorithm and the 3-state Kalman algorithm has shown their consistency. However, as it was demonstrated experimentally, the FIR filter produces a smaller error and a lower Allan variance for the sawtooth noise. Moreover, it may be advanced further to be optimal in a sense of a minimum MSE that is currently under investigation.

Appendix A Coefficients of Matrix (25) The coefficients for (25) are calculated by: dm =

N −1  i=0

im =

1 [Bm+1 (N ) − Bm+1 ] , m+1

where Bn (x) is the Bernoulli polynomial and Bn = Bn (0) is the Bernoulli coefficient. For low orders, Bn (x) may be found in the reference books. For high orders, the following recurrent relation is valid: Bn (x) = n Bn−1 (x)dx + Bn .

(b)

(c) Fig. 6. Long-term measurement and estimation of the crystal clock TIE with the 3-state unbiased FIR algorithm: (a) TIE, (b) fractional frequency offset, and (c) linear fractional frequency drift rate.

Several low order coefficients dm are given below: d0 = N, N (N − 1) d1 = , 2 N (N − 1)(2N − 1) d2 = , 6 N 2 (N − 1)2 , d3 = 4 N (N − 1)(2N − 1)(3N 2 − 3N − 1) d4 = , 30 N 2 (N − 1)2 (2N 2 − 2N − 1) d5 = , 12 N (N − 1)(2N − 1)(3N 4 − 6N 3 + 3N + 1) d6 = . 42 For large horizon, N 1, the coefficients dm may be m+1 calculated by dm |N 1 ∼ = Nm+1 .

shmaliy: proposed finite impulse response filter

Appendix B Three-State Kalman Filtering Algorithm In the state space, the TIE model (1) is given by: ⎤⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ x(n − 1) w1 (n, τ ) x(n) 1 τ τ 2 /2 ⎣y(n)⎦ = ⎣0 1 τ ⎦ ⎣y(n − 1)⎦ + ⎣w2 (n, τ )⎦ 00 1 z(n − 1) w3 (n, τ ) z(n) x(n) = Ax(n − 1) + w(n, τ ), and (5) becomes, assuming a single receiver, ⎤ ⎡

 x(n) ξ(n) = 1 0 0 ⎣y(n)⎦ + v(n), z(n) ξ(n) = Cx(n) + v(n). The noises w(n, τ ) and v(n) are mean zero and jointly uncorrelated. The sawtooth noise v(n) has a uniform distribution p(v) = 1/2vmax and correlated increments. Its white Gaussian approximation has a variance V = σv2 = vmax 1 2 2 2vmax −vmax v dv = vmax /3. The autocorrelation matrix of the white Gaussian noise w(n) is given by [19]: ⎡ 2 4 3 2⎤ q1 + q23τ + q320τ q22τ + q38τ q36τ 3 2 ⎢ q3 τ ⎥ , Ψ = τ ⎣ q22τ + q38τ ⎦ q2 + q33τ 2 q3 τ 2 6

q3 τ 2

q3

in which the diffusion coefficients q’s, namely q1 , q2 , and q3 , specify the white FM noise (WHFM), white random walk FM noise (WRFM), and white random run FM noise (RRFM), respectively, in the τ -domain power law. The linear Kalman filtering algorithm reads as follows. Enter the q’s, Rn−1 , and x ˆn−1 , then calculate recursively: ˜ n = ARn−1 AT + Ψ, R ˜ n CT (CR ˜ n CT + V )−1 , Kn = R x ˆn = Aˆ xn−1 + Kn (ξn − CAˆ xn−1 ), ˜ Rn = (I − Kn C)Rn .

Acknowledgment The author would like to thank Dr. Raymond Filler of the U.S. Army Research, Development and Engineering Command (RDECOM) CERDEC; Dr. Charles Greenhall of the Jet Propulsion Laboratory (JPL), California Institute of Technology; Dr. Judah Levine of the National Institute of Standards and Technology (NIST); and two anonymous reviewers for valuable comments and remarks. References [1] W. Lewandowski, G. Petit, and C. Thomas, “Precision and accuracy of GPS time transfer,” IEEE Trans. Instrum. Meas., vol. 42, no. 2, pp. 474–479, Apr. 1993. [2] F. Meyer, “Common-view and melting-pot GPS time transfer with the UT+,” in Proc. 32nd Annu. Precise Time and Time Interval Meeting, 2000, pp. 147–156.

869 [3] D. W. Allan, J. E. Gray, and H. E. Machlan, The National Bureau of Standards Atomic Time Scale: Generation, Stability, Accuracy and Accessibility. NBS Monograph 140, Time and Frequency: Theory and Fundamentals, National Institute of Standards and Technology, 1974, pp. 205–231. [4] R. M. Hambly and T. A. Clark, “Critical evaluation of the Motorola M12+ GPS timing receiver vs. the master clock at the United States Naval Observatory, Washington, DC,” in Proc. 34th Annu. Precise Time and Time Interval Meeting, 2002, pp. 109–115. [5] J. Levine, “Time transfer using multi-channel GPS receivers,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 46, no. 2, pp. 392–398, Mar. 1999. [6] D. W. Allan and J. A. Barnes, “Optimal time and frequency transfer using GPS signals,” in Proc. 36th Annu. Freq. Contr. Symp., 1982, pp. 378–387. [7] P. V. Tryon and R. H. Jones, “Estimation of parameters in models for cesium beam atomic clocks,” J. Res. National Bureau of Standards, pp. 3–16, 1983. [8] J. W. Chaffee, “Relating the Allan variance to the diffusion coefficients of a linear stochastic differential equation model for precision oscillators,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 34, no. 6, pp. 655–658, Nov. 1987. [9] S. R. Stein and R. L. Filler, “Kalman filter analysis for real time applications of clocks and oscillators,” in Proc. 42th Annu. Freq. Contr. Symp., 1988, pp. 447–452. [10] J. A. Barnes, R. H. Jones, P. V. Tryon, and D. W. Allan, “Stochastic models for atomic clocks,” in Proc. 14th Annu. Precise Time and Time Interval Meeting, 1982, pp. 295–306. [11] L. Breakiron, “Timescale algorithms combining cesium clocks and hydrogen masers,” in Proc. 23rd Annu. Precise Time and Time Interval Meeting, 1991, pp. 297–305. [12] C. Greenhall, “Kalman plus weights: A time scale algorithm,” in Proc. 33rd Annu. Precise Time and Time Interval Meeting, 2001, pp. 445–454. [13] K. Senior, P. Koppang, and J. Ray, “Developing an IGS time scale,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 50, no. 6, pp. 585–593, Jun. 2003. [14] L. Galleani and P. Tavella, “On the use of the Kalman filter in timescales,” Metrologia, vol. 40, pp. 326–334, 2003. [15] C. J. Masreliez, “Approximate non-Gaussian filtering with linear state and observation relations,” IEEE Trans. Automat. Contr., vol. 20, pp. 107–110, 1975. [16] H. Wu and G. Chen, “Suboptimal Kalman filtering for linear systems with Gaussian-sum type of noise,” Math. Comput. Model., vol. 29, no. 5, pp. 101–125, 1999. [17] O. K. Kwon, W. H. Kwon, and K. S. Lee, “FIR filters and recursive forms for discrete-time state-space models,” Automatica, vol. 25, pp. 715–728, 1989. [18] W. H. Kwon, P. S. Kim, and P. Park, “A receding horizon Kalman FIR filter for discrete time-invariant systems,” IEEE Trans. Automat. Contr., vol. 44, no. 9, pp. 1787–1791, 1999. [19] W. H. Kwon, P. S. Kim, and S. H. Han, “A receding horizon unbiased FIR filter for discrete-time state space models,” Automatica, vol. 38, pp. 545–551, 2002. [20] Y. S. Shmaliy, “A simple optimally unbiased MA filter for timekeeping,” IEEE Trans. Ultrason., Ferroelect., Freq. Contr., vol. 49, no. 6, pp. 789–797, Jun. 2002. [21] P. Heinonen and Y. Neuvo, “FIR-median hybrid filters with predictive FIR structures,” IEEE Trans. Acoust. Speech Signal Processing, vol. 36, no. 6, pp. 892–899, Jun. 1988.

Yuriy S. Shmaliy (M’96–SM’00) was born January 2, 1953. He received the B.S., M.S., and Ph.D. degrees in 1974, 1976, and 1982, respectively, from Aviation Institute of Kharkiv, Kharkiv, Ukraine, all in electrical engineering. In 1992 he received the Doctor of Technical Sc. degree from the Railroad Academy of Kharkiv, Kharkiv, Ukraine. In March 1985, he joined the Kharkiv Military University, Kharkiv, Ukraine. He served as professor beginning in 1986. In 1999, he joined the Kharkiv National University of

870

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 53, no. 5, may 2006

Radio Electronics, Kharkiv, Ukraine. Since November 1999, he has been with the Guanajuato University of Mexico, Salamanca, Gto., Mexico, as a professor. Dr. Shmaliy has 178 papers and 80 patents. He was awarded a title, Honorary Radio Engineer of the USSR in 1991; was listed in Marquis Who’s Who in the World in 1998; and was listed in Outstanding People of the 20th Century, Cambridge, England in 1999. He is a member of several professional societies and organizing and program committees of International Symposia. His current interests include the statistical theory of precision resonators and oscillators, optimal estimation, and stochastic signal processing for frequency and time.