- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

International Journal of Navigation and Observation

Volume 2008 (2008), Article ID 261384, 8 pages

http://dx.doi.org/10.1155/2008/261384

## GPS Composite Clock Analysis

Analytical Graphics, Inc., 220 Valley Creek Blvd, Exton, PA 19341, USA

Received 30 June 2007; Accepted 6 November 2007

Academic Editor: Demetrios Matsakis

Copyright © 2008 James R. Wright. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The GPS composite clock defines GPS time, the timescale used today in GPS operations. GPS time is illuminated by examination of its role in the complete estimation and control problem relative to UTC/TAI. The phase of each GPS clock is unobservable from GPS pseudorange measurements, and the mean phase of the GPS clock ensemble (GPS time) is unobservable. A new and useful *observability* definition is presented, together with new observability theorems, to demonstrate explicitly that GPS time is unobservable. Simulated GPS clock phase and frequency deviations, and simulated GPS pseudorange measurements, are used to understand GPS time in terms of Kalman filter estimation errors.

#### 1. Introduction

GPS time is *created* by processing GPS pseudorange measurements with the operational GPS Kalman
filter. Brown [2] refers to the object created by the Kalman filter as the GPS
composite clock, and to GPS time as the implicit ensemble mean phase of the GPS
composite clock. The fundamental goal by the USAF and the USNO is to control
GPS time to within a specified bound of UTC/TAI. (I refer to TAI/UTC understanding
that UTC has an accumulated discontinuity (a sum of leap seconds) when compared
to TAI. But unique two-way transformations between TAI and UTC have been in
successful operational use since 1972. I have no need herein to further
distinguish between TAI and UTC.) I present here a quantitative analysis of the
GPS composite clock, derived from detailed simulations and associated graphics.
GPS clock diffusion coefficient values used here were derived from Allan
deviation graphs presented by Oaks et al. [12] in 1998. I refer to them as
“realistic,” and in the sequel I claim “realistic” results from their use.
Figure 1 presents their diffusion coefficient values and my derivation of
associated Allan deviation lines.

My interest in the GPS composite clock derives from my interest in performing real-time orbit determination for GPS NAVSTAR spacecraft from ground receiver pseudorange measurements. (James R Wright is the architect of ODTK (Orbit Determination Tool Kit), a commercial software product offered by Analytical Graphics, Inc. (AGI).) The estimation of NAVSTAR orbits would be incomplete without the simultaneous estimation of GPS clock parameters. I use simulated GPS clock phase and frequency deviations, and simulated GPS pseudorange measurements, to study Kalman filter estimation errors.

This paper was first prepared for TimeNav'07 [20]. I am indebted to Charles Greenhall (JPL) for encouragement and help in this work.

#### 2. The Complete Estimation and Control Problem

The USNO operates two UTC/TAI master clocks, each of which provides access to an estimate of UTC/TAI in real time (1 pps). One of these clocks is maintained at the USNO, and the other is maintained at Schriever Air Force Base in Colorado Springs. This enables the USNO to compare UTC/TAI to the phase of each GPS orbital NAVSTAR clock via GPS pseudorange measurements, by using a UTC/TAI master clock in a USNO GPS ground receiver. Each GPS clock is a member of (internal to) the GPS ensemble of clocks, but the USNO master clock is external to the GPS ensemble of clocks. Because of this, the difference between UTC/TAI and the phase of each NAVSTAR GPS clock is observable. This difference can be (and is) estimated and quantified. The root mean square (RMS) on these differences quantifies the difference between UTC/TAI and GPS time. Inspection of the differences between UTC/TAI and the phase of each NAVSTAR GPS clock enables the USNO to identify GPS clocks that require particular frequency-rate control corrections. Use of this knowledge enables the USAF to adjust frequency rates of selected GPS clocks. Currently, the USAF uses an automated bang-bang controller on frequency-rate. (According to Bill Feess, an improvement in control can be achieved by replacing the existing “bang-bang controller” with a “proportional controller.”)

#### 3. Stochastic Clock Physics

The most significant stochastic clock physics are understood in terms of Wiener processes and their integrals. Clock physics are characterized by particular values of clock-dependent diffusion coefficients, and are conveniently studied with aid of a relevant clock model that relates diffusion coefficient values to their underlying Wiener processes. For my presentation here I have selected “The clock model and its relationship with the Allan and related variances" presented as an IEEE paper by Zucca and Tavella [19] in 2005. Except for FM flicker noise, this model captures the most significant physics for all GPS clocks. I simulate and validate GPS pseudorange measurements using simulated phase deviations and simulated frequency deviations, according to Zucca and Tavella.

#### 4. Kalman Filters

I present my approach for the *optimal* sequential estimation of clock deviation states and their error covariance
functions. Sequential state estimates are generated recursively from two
multidimensional stochastic update functions, the *time update* (TU) and
the *measurement update* (MU). The TU moves the state estimate and
covariance forward with time, accumulating integrals of random clock deviation
process noise in the covariance. The MU is performed at a fixed measurement
time where the state estimate and covariance are corrected with new observation
information.

The sequential estimation of GPS clock deviations requires the development of a linear TU and nonlinear MU. The nonlinear MU must be linearized locally to enable application of the linear Kalman MU. Kalman's MU [9] derives from Sherman's theorem [11, 15, 16], Sherman's theorem derives from Anderson's theorem [1], and Anderson's theorem derives from the Brunn-Minkowki inequality theorem [5, 17]. The theoretical foundation for my linearized MU derives from these theorems.

##### 4.1. Initial Conditions

Initialization of all sequential estimators requires the use of an initial state estimate column matrix and an intial state estimate error covariance matrix for time .

##### 4.2. Kalman Filter: Linear TU and Linear MU

Derivation and calculation for the discrete-time Kalman filter, linear in both TU and MU, is best presented by Meditch [11, Chapter 5].

##### 4.3. Linear TU and Nonlinear MU

The simultaneous sequential estimation of GPS clock phase and frequency deviation parameters can be studied with the development of a linear TU and nonlinear MU for the clock state estimate subset. This is useful to study clock parameter estimation, as demonstrated in Section 6.

Let denote an column matrix of state estimate components, where the left subscript denotes state epoch and the right subscript denotes time-tag for the last observation processed, where . Let denote an associated square symmetric state estimate error covariance matrix (positive eigenvalues).

###### 4.3.1. Linear TU

For , the propagation of the true unknown matrix state is given bywhere is called the process noise matrix. Propagation of the known matrix state estimate is given bybecause the conditional mean of is zero. Propagation of the known matrix state estimate error covariance matrix is given bywhere the matrix is called the process noise covariance matrix (see [19] for concrete clock examples of and ).

###### 4.3.2. Nonlinear MU

Calculate the matrix filter gain : The filter measurement update state estimate matrix , due to the observation , is calculated withwhere is the scalar variance on the observation residual , and is a nonlinear function of . Define the error in :Define the state estimate error covariance matrix with Bucy and Joseph [3, page 141] recommend that should be calculated withwhereEquations (8) and (9) reduce to the form given by Kalman:Calculation of by (8) and (9) is numerically stable, whereas the Kalman form calculation is not.

##### 4.4. Nonlinear TU and Nonlinear MU

Refer to Section 4.3.2 for the nonlinear MU.

###### 4.4.1. Nonlinear TU

The nonlinear TU always spans a nonempty time interval and requires the use of a numerical state estimate integrator . Given an initial time , a final time , and a force model , then propagates the state estimate from to using forces to get . That is,This can be shortened to writewhere the use of forces is tacitly implied. Thus, is a column matrix with elements:

##### 4.5. Kalman Filter Advantage

Severe computational problems are incurred in any attempt to estimate unobservable states using iterated batch least squares methods or iterated maximum likelihood methods for navigation, because state-sized inversions of singular matrices are required. Here the Kalman filter is distinguished in that estimates of unobservable states can be created and used without matrix inversion problems because the Kalman filter MU is free of state-sized matrix inversions.

By design, one typically estimates observable states. But the Kalman filter enables one to create unobservable states. The USAF chose to create unobservable GPS clock parameter states for construction of GPS time.

#### 5. Observability

I have defined *observability* in terms of a Kalman
filter formulation, and I have proved simple theorems related thereto. My
definition of observability is different than Kalman's definition and, unlike
Kalman's definition, is directly applicable to covariance matrices derived from
a Kalman filter.

##### 5.1. Definition

If the state estimate error variance of a particular
state estimate component is reduced by processing an observation, then that
state estimate component is *observable* to that observation. Otherwise,
that state estimate component is *not observable* (*unobservable*) to
that observation.

Theorem 1. *If every component of the row matrix of
measurement-state partial derivatives is zero at time , then every component of the state estimate is unobservable at time .*

*Proof. * implies that according to (10). Thus none of the variances of are reduced due
to processing the observation . Then by definition, is unobservable in every component.

Theorem 2. *Given values for scalars , , at time , and given that , then the scalar state estimate is observable at time .*

*Proof. *The obvious inequality implies thatMultiply through by : Add 1: Multiply through by :
Now use (4) and (10) to write
Insert (18) into the inequality (17) to get the result:
Thus the variance is reduced due
to processing the observation . Then the scalar state is observable
by definition.

##### 5.2. Theoretical Foundation

These theorems are referred to expressions given by Kalman for filter gain and covariance , see (4) and (10). Kalman's expressions are derived from the rigorous theorem chain provided by Sherman, Anderson, and Brunn-Minkowski—the theoretical foundation is deep.

##### 5.3. Determine Observability Directly

Given an optimal sequential estimator,
given a particular collection of applicable observations (real or simulated),
and given realistic state estimate error covariance matrices and at each time , apply the definition of observability directly (note
that this is impossible using Kalman's definition of observability) to
distinguish between *observable* and *unobservable* state elements.
An optimal sequential estimator is designed to eliminate significant aliasing
between estimated state elements, and thus enables
this distinction.

#### 6. Unobservable GPS Clock States

*GPS time* is created by the operational USAF Kalman filter by processing GPS pseudorange
observations. *GPS time* is the mean phase of an ensemble of many GPS
clocks, and yet the clock phase of every operational GPS clock is unobservable
from GPS pseudorange observations, as demonstrated below. GPS NAVSTAR orbit
parameters are observable from GPS pseudorange observations. The USAF Kalman
filter simultaneously estimates orbit parameters and clock parameters from GPS
pseudorange observations, so the state estimate is partitioned in this manner
into a subset of unobservable clock parameters and a subset of observable orbit
parameters. This partition is performed by application of Sherman's theorem in
the MU.

##### 6.1. GPS Pseudorange Representation

Let denote time of radio wave transmission for the th NAVSTAR clock, and let denote time of radio wave receipt for the th ground station clock. (Refer all times to a coordinate time, e.g., to GPS time. Appropriate transformations between proper time and coordinate time must be performed in the operational algorithms, but state estimate observability is independent of relativity, so observability can be defined and discussed independent of relativity.) Let and denote Kalman filter estimation errors in clock phase for and . Define time of transmission difference and time of receipt difference :Thus,Equation (21) present and as additive combinations of deterministic times and and Kalman filter estimation errors in clock phase and . Define the one-way GPS pseudorange measurement : Insert (21) into (22):where is speed of light in vacuum. Define Then,where is deterministic and is random.

##### 6.2. Partition of Kalman Filter Estimation Errors

Let denote the
phase component of Kalman filter estimation error that is *common* to
every GPS ensemble clock, when it exists. Define phase differences and with for ground station and NAVSTAR . Then Kalman filter estimation errors , , for ground station clocks and , , for NAVSTAR clocks have the additive partition:(This partition
was introduced by Brown [2].)

##### 6.3. The Common Random Phase Component Is Unobservable

Insert (28) into (25): Insert (29) into (26):Thus, the random phase component that is common to the Kalman filter estimation error for every ensemble clock has vanished in the range representation . Variations in cannot cause variations in :because the partial derivative is zero: An application of Theorem 1 to (32) demonstrates that is unobservable from .

But the architect who designs the complete estimator
must design an *optimal* NAVSTAR orbit estimator to prevent aliasing from
NAVSTAR orbit estimation errors into . It helps to know that there is no coupling between
the orbit and in the complete
state transition function. I have provided a new method herein to identify this
aliasing, and I have provided suggestions on where to look for inadequate
modeling that would be the source of this aliasing. See Section 9.

##### 6.4. Independent Random Phase Components Are Observable

The independent phase deviations and are observable to GPS pseudorange observations because their partial derivatives are nonzero:Estimation of and by the Kalman filter will reduce their error variances.

##### 6.5. Partition of KF1 Estimation Errors

Subtract estimated clock deviations from simulated
(true) clock deviations to define and quantify Kalman filter (KF1) estimation
errors. Adopt Brown's additive partition of KF1 estimation errors into two
components. I refer to the first component as the unobservable error common to
each clock (UECC), and to the second component as the observable error
independent for each clock (OEIC). (Observability is meaningful here only when
processing *simulated* GPS pseudorange data.) See (28). is the UECC, is the OEIC for
ground station clocks and is the OEIC for
NAVSTAR clocks. On initialization of KF1, the variances on the UECC and OEIC
are identical. On processing the first GPS pseudorange measurements with KF1
the variances on both fall quickly. But with continued measurement processing
the variances on the UECC increase without bound while the variances on the
OEIC appoach zero asymptotically.

For simulated GPS pseudorange data I create an optimal sequential estimate of the UECC by application of a second Kalman filter KF2 to pseudomeasurements defined by the phase components of KF1 estimation errors.

Since there is no physical process noise on the UECC, an estimate of the UECC can also be achieved using a batch least squares estimation algorithm on the phase components of KF1 estimation errors–-demonstrated previously by Greenhall [7]. (I apply sufficient process noise covariance for KF2 to mask the effects of double-precision computer word truncation. Without this, KF2 does diverge.)

##### 6.6. Unobservable Error Common to Each Clock

There are at least four techniques to estimate the UECC when simulating GPS pseudorange data. First, one could take the sample mean of KF1 estimation errors across the clock ensemble at each time and form a sample variance about the mean; this would yield a sequential sampling procedure, but where each mean and variance is sequentially unconnected. Second, one can employ Ken Brown's implicit ensemble mean (IEM) and covariance; this is a batch procedure requiring an inversion of the KF1 covariance matrix followed by a second matrix inversion of the modified covariance matrix inverse; this is not a sequential procedure. Third, one can adopt the new procedure by Greenhall [7] wherein KF1 phase estimation errors are treated as pseudomeasurements, and are processed by a batch least squares estimator to obtain optimal batch estimates and covariance matrices for the UECC. Fourth, one can treat the KF1 phase estimation errors as pseudomeasurements, invoke a second Kalman filter (KF2), and process these phase pseudomeasurements with KF2 to obtain optimal sequential estimates and variances for the UECC. I have been successful with this approach. Figure 3 presents an ensemble of “realistic" KF1 phase estimation errors, overlaid with “realistic" KF2 sequential estimates of UECC in phase. (By “realistic" I refer to realistic clock diffusion coefficient values.)

##### 6.7. Observable Error Independent for Each Clock

At each applicable time subtract the estimate of the UECC from the KF1 phase deviation estimate, for each particular GPS clock, to estimate the OEIC in phase for that clock. During measurement processing, the OEIC is contained within an envelope of a few parts of a nanosecond (see Figure 4).

Figure 4 presents a graph of two cases of the OEIC for ground station clock S1. For the blue line of Figure 4, intervals of link visibility and KF1 range measurement processing are clearly distinguished from propagation intervals with no measurements. During measurement processing, the observable component of KF1 estimation error is contained within an envelope of a few parts of a nanosecond.

Calculation of the sequential covariance for the OEIC requires a matrix value for the cross-covariance between the KF1 phase deviation estimation error and the UECC estimation error at each time. I have not yet been able to calculate this cross-covariance.

#### 7. Allan Variance and PPN Relations

##### 7.1. Allan Coefficients versus Diffusion Coefficients

Denote as clock averaging time, as Allan variance, as Allan's FMWN coefficient, as Allan's FMRW coefficient, as the FMWN diffusion coefficient, and as the FMRW diffusion coefficient. Then,where

##### 7.2. Proportionate Process Noise (PPN)

Let denote a variable to identify each GPS clock in an ensemble of clocks. For each clock define the ratio between diffusion coefficients and :Then PPN is defined when, for each GPS clock and each associated ratio , we have

##### 7.3. Case 12

The calculation of , , according to the diffusion coefficient values
presented in Figure 1 shows that PPN is *not* satisfied for Case 12:

#### 8. Kalman Filters KF1 and KF2

I have simulated GPS pseudorange measurements for two GPS ground station clocks S1 and S2, and for two GPS NAVSTAR clocks N1 and N2. Here I set simulated measurement time granularity to 30 s for the set of all visible link intervals. Visible and non-visible intervals are clearly evident in the blue line of Figure 4. I set the scalar root-variance for both measurement simulations and Kalman filter KF1 to cm. Typically m for GPS pseudorange, but when carrier phase measurements are processed simultaneously with pseudorange, the root-variance is reduced by two orders of magnitude. So the use of cm enables me to quantify lower performance bounds for the simultaneous processing of both measurement types.

##### 8.1. Create GPS Clock Ensemble

Typically, one processes measurements with a Kalman
filter to derive sequential estimates of a multidimensional *observable* state. Instead, here I imitate the GPS operational procedure and process
simulated GPS pseudorange measurements with KF1 to *create* a sequence of
unobservable multidimensional clock state estimates. Clock state components are
unobservable from GPS pseudorange measurements. See Figure 2 for an example of an
ensemble of estimated unobservable clock phase deviation state components
created by KF1.

###### 8.1.1. Sherman's Theorem

GPS time, the unobservable GPS clock ensemble mean
phase, is created by the use of Sherman's theorem [11, 18] in the USAF Kalman
filter measurement update algorithm on GPS range measurements. Satisfaction of
Sherman's Theorem guarantees that the mean-squared state estimate error on each *observable* state estimate component is minimized. But the mean-squared
state estimate error on each *unobservable* state estimate component is
not reduced. Thus the unobservable clock phase deviation state estimate
component common to every GPS clock is isolated by application of Sherman's
theorem. An ensemble of unobservable state estimate components is thus created
by Sherman's theorem—see Figure 3 for an example.

##### 8.2. Initial Condition Errors

A significant result emerges due to the modeling of Kalman filter (KF1) initial condition errors in phase and frequency. Initial estimated clock phase deviations are significantly displaced by the KF1 initial condition errors in phase. As time evolves estimated clock phase deviation magnitudes diverge continuously and increasingly when referred to true (simulated) phase deviations, and this is due to filter initial condition errors in frequency. See Figure 2 for an example.

#### 9. Identify Nonclock Modeling Errors

My interest in the GPS NAVSTAR (SV) orbit
determination problem, combined with that of the clock parameter estimation
problem, has enabled the identification of a useful diagnostic tool: given
realistic values for diffusion coefficients for each of the real GPS clocks,
then quantitative upper bounds can be calculated on OEIC magnitudes. These
calculations require the use of a rigorous *simulator*. Existence of
significant cross-correlations between GPS clock phase errors and other
nonclock GPS estimation modeling errors enables significant aliasing into GPS
clock phase estimates during operation of KF1 on *real* data. But given
rigorous quantitative upper bounds on OEIC magnitudes, then significant
violation of these bounds when processing *real* GPS pseudorange and
carrier phase data identifies nonclock modeling errors related to the GPS
estimation model. Modeling error candidates here include NAVSTAR orbit force
modeling errors, ground antenna modeling errors (multipath), and tropospheric
modeling errors. NAVSTAR orbit force modeling errors include those of solar
photon pressure, albedo, thermal dump, and propellant outgassing. The accuracy
of this diagnostic tool depends on the use of realistic clock diffusion
coefficient values and a rigorous clock model simulation capability.

#### 10. Observable Clocks

In an earlier version of my paper, I reported on KF1 validation results where clock S1 was specified as a TAI/UTC clock, external to the GPS clock ensemble consisting of S2, N1, and N2. This brought observablity (see Sections 5 and 6 herein) to S2, N1, and N2 clock states from GPS pseudorange measurements, drove clocks S2, N1, and N2 immediately to the TAI/UTC timescale, and enabled a clean validation of my filter implementation. Also it raised the question: why not the same thing for the real GPS clock ensemble? Discussions with Ed Powers (USNO) and Bill Feess (Aerospace Corporation) reveal that this approach was tried and discarded after the difficulty in recovery from an uplink hardware failure was blamed on the use of a single TAI/UTC Master Clock. This issue was resolved with Kenneth Brown’s introduction of the implicit ensemble mean. The mean phase (GPS time) of the GPS clock ensemble will remain unobservable to GPS pseudorange measurements in the USAF Kalman filter for the foreseeable future.

#### References

- T. W. Anderson, “The integral of a symmetric unimodal function over a symmetric convex set and some probability inequalities,”
*Proceedings of the American Mathematical Society*, vol. 6, no. 2, pp. 170–176, 1955. View at Publisher · View at Google Scholar - K. R. Brown, “The theory of the GPS composite clock,” in
*Proceedings of the 4th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS '91)*, pp. 223–241, Albuquerque, NM, USA, September 1991. View at Google Scholar - R. S. Bucy and P. D. Joseph,
*Filtering for Stochastic Processes with Applications to Guidance*, Interscience, New York, NY, USA, 1968. - W. Feess, “The Aerospace Corporation,” 2006, private Communications. View at Google Scholar
- R. J. Gardner, “The Brunn-Minkowski inequality,”
*Bulletin of the American Mathematical Society*, vol. 39, no. 3, pp. 355–405, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - C. A. Greenhall, 2006-2007, private Communications.
- C. A. Greenhall, “A Kalman filter clock ensemble algorithm that admits measurement noise,”
*Metrologia*, vol. 43, no. 4, pp. S311–S321, 2006. View at Publisher · View at Google Scholar - S. T. Hutsell, “Relating the hadamard variance to MCS Kalman filter clock estimation,” in
*Proceedings of the 27th Annual Precise Time and Time Interval (PTTI) Applications and Planning Meeting*, p. 293, San Diego, Calif, USA, December 1995. View at Google Scholar - R. E. Kalman, “New methods in wiener filtering theory,” in
*Proceedings of the 1st Symposium on Engineering Applications of Random Function Theory and Probability*, J. L. Bogdanoff and F. Kozin, Eds., John Wiley & Sons, New York, NY, USA, 1963. View at Google Scholar - D. Matsakis, November 2007, private Communication.
- J. S. Meditch,
*Stochastic Optimal Linear Estimation and Control*, McGraw-Hill, New York, NY, USA, 1969. - O. J. Oaks, T. B. McCaskill, M. M. Largay, W. G. Reid, and J. A. Buisson, “Performance of GPS on-orbit NAVSTAR frequency standards and monitor station time references,” in
*Proceedings of the 30th Annual Precise Time and Time Interval (PTTI) Meeting*, pp. 135–143, Reston, Va, USA, December 1998. View at Google Scholar - E. Powers, 2006, private Communications.
- W. Riley, December 2006, private Communications, PTTI meeting.
- S. Sherman, “A theorem on convex sets with applications,”
*The Annals of Mathematical Statistics*, vol. 26, no. 4, pp. 763–767, 1955. View at Google Scholar - S. Sherman, “Non-mean-square error criteria,”
*IEEE Transactions on Information Theory*, vol. 4, no. 3, pp. 125–126, 1958. View at Google Scholar - E. M. Stein and R. Shakarchi,
*Real Analysis*, Princeton University Press, Princeton, NJ, USA, 2005. - J. R. Wright, “Sherman's theorem,” in
*The Malcolm D. Shuster Astronautics Symposium (AAS '05)*, Grand Island, NY, USA, June 2005. View at Google Scholar - C. Zucca and P. Tavella, “The clock model and its relationship with the allan and related variances,”
*IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control*, vol. 52, no. 2, pp. 289–295, 2005. View at Publisher · View at Google Scholar - J. R. Wright, “GPS composite clock analysis,” in
*IEEE International Frequency Control Symposium, European Frequency and Time Forum*, pp. 523–528, Geneva, Switzerland, June 2007. View at Google Scholar