[time-nuts] Re: GPSDO Control loop autotuning

Thu Apr 21 17:25:08 UTC 2022

Kalman filters are optimal only if they work with statistical processes
that can be described with derivatives or integrals, i.e. white phase noise
(f^0), white FM (f^-2), random walk FM (f^-4) and so on.
Most of VCOs are heavily flicker FM dominated

Il giorno gio 21 apr 2022 alle ore 19:09 Erik Kaashoek <erik at kaashoek.com>
ha scritto:

> Hi Tobias,
> Using the excellent material found here:
>
> https://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb
> I went ahead and implemented the first order filter as described in
> chapter 8.
> The simulation was changed to use a real recorded PPS as measurement
> noise and a real recorded TIC of a TCXO as the to be estimated data and
> a Vtune steering an added frequency
> The STDEV of the PPS was 5.6e-9
> The kalman filter predicts the phase and frequency fairly well
> STDERR of phase is 4e-9
> STDERR of freq is 1.75e-9
>
> The loop was then closed with a PI controller tasked to minimize the
> phase error.
> Tuning was a bit of a problem but the current performance is
> STDEV phase error 1.25e-8
> STDEV frequency error 1.2e-9
>
> Unfortunately the resulting ADEV is still worse compared to a well tuned
> PI controller acting directly on the TIC of the TCXO.
> More work to do
> Erik.
>
> On 21-4-2022 9:44, Pluess, Tobias via time-nuts wrote:
> > Hi Erik,
> >
> > yes I also found this webpage, and the examples are good.
> > I had the exact same idea as you, except I thought that the model for the
> > VCO output is either phase only or phase and frequency. Yours is even a
> bit
> > more sophisticated, thats quite cool!
> > And yes, of course it will be possible to have a one dimensional filter,
> > even if your state space equations have more than one state. However, in
> > this case, A, B and C are a matrix and two vectors (sometimes also called
> > F, G and H). But the output that is measured ("observed") is just scalar,
> > i.e. 1D.
> > I am not sure whether your state extrapolation is correct, I don't see it
> > at the moment unfortunately. But I will think about it a bit more.
> >
> > To your question, about how to close the loop.
> > Well. What do you do when you don't have the alpha-beta-gamma filter?
> > you take your phase error readings and feed them into some sort of
> > controller, for example a PID controller, do you.
> > Now think of the alpha-beta-gamma filter (or also the Kalman filter) as
> > some sort of "prefilter". (do you remember when you asked me why I
> > prefilter my 1PPS measurements? :-D at that time, I had a 2pole IIR
> filter
> > with adjustable time constant as prefilter). So you feed your actual 1PPS
> > measurements into the alpha-beta-gamma filter, and the filter outputs the
> > 1PPS measurements, but with noise mostly removed (hopefully!). Then, use
> > this as the error signal for your controller.
> >
> > As how to choose the alpha, beta and gamma coefficients, you need to know
> > the variance of your input signal. i.e. determine the standard deviation
> of
> > the 1PPS readings and square that. I am not sure at the moment how to
> find
> > the coefficients from the variance, but I believe somewhere in the
> > kalmanfilter.net tutorial, there is an equation for it if I remember
> > correctly.
> > Also worth mentioning is the fact that some GPS modules also output an
> > estimate of the timing accuracy. This can also be used as an additional
> > input to calculate the coefficients and adapt them, as receiving
> conditions
> > change.
> >
> > does that help a bit?
> > I am unfortunately not yet so far as you, I have made some calculations
> > using matlab, but I am not sure whether they are correct. I do see some
> > filtering effect of the Kalman filter, though. (I used the Kalman filter,
> > with the error covariance matrices, not the alpha-beta-gamma filter, like
> > you did, but maybe this is a better idea!). Are you interested in sharing
> > your excel sheet?
> >
> > thanks,
> > best
> > Tobias
> > HB9FSX
> >
> >
> >
> >
> > On Wed, Apr 20, 2022 at 1:18 AM Erik Kaashoek <erik at kaashoek.com> wrote:
> >
> >> Tobias,
> >> Your Kalman post triggered me to study this "dummy" level introduction
> >> to Kalman filters: https://www.kalmanfilter.net/
> >> and now I'm trying to write the Kalman filter.
> >> An example of a Kalman filter described in the link above uses distance,
> >> speed and acceleration as predictors for a one dimensional Kalman filter
> >> So I wondered if it would be possible to do the same for the phase, eg
> >> use phase, frequency and drift in the prediction model and still have a
> >> one dimensional model?
> >> The Kalman filter thus uses the normalized phase as input (the measured
> >> time difference between the 1 PPS and the 1 second pulse derived from
> >> the VCO output) to predict the next phase, frequency and drift
> >> p = phase
> >> f = frequency (d Phase/d T)
> >> d = drift (d Frequency / d T)
> >> m = measured phase
> >> with a 1 second interval between observations the state extrapolation
> >> (or prediction) is
> >> p[k+1]  = p[k] + f[k] + 0.5*d[k]
> >> f[k+1] = f[k] + d[k]
> >> d[k+1] = d[k]
> >> and the state update is
> >> p[k] = p[k-1] + alpha*( m - p[k-1])
> >> f[k] = f[k-1] + beta*(m-p[k-1])
> >> d[k] = d[k-1] + gamma*(m-p[k-1])/0.5)
> >>
> >> In a small simulation (in excel) I used a VCO that initially is on the
> >> correct frequency, after some seconds Vtune has a jump and after some
> >> more seconds the Vtune starts to drift.
> >> The Vtune of the VCO can have noise(system error) and the m can have
> >> noise (measurement error).
> >> The p,f and m nicely track the input values even when measurement noise
> >> is added but after this I'm stuck as I do not yet understand how to
> >> calculate the alpha, beta and gamma values from the observed measurement
> >> error as the above website does not yet explain how to do this.
> >> And I do not understand how to close the loop to either keep the phase
> >> error or frequency error as low as possible.
> >> So much to learn.
> >> Erik.
> >>
> >> On 16-4-2022 22:44, Pluess, Tobias via time-nuts wrote:
> >>> Hallo all,
> >>>
> >>> In the meantime I had to refresh my knowledge about state-space
> >>> representation and Kalman filters, since it was quite a while ago
> since I
> >>> had this topic.
> >>>
> >>> So I looked at the equations of the Kalman filter. To my understanding,
> >> we
> >>> can use it like an observer, and instead of using the phase error and
> >>> feeding it to the PI controller, we use the output of the Kalman filter
> >> as
> >>> input to the PI controller. And the Kalman filter gets its input from
> the
> >>> phase error, but we also tell it how much variance this phase error
> has.
> >>> Luckily, the GPS module outputs an estimate of the timing accuracy, so
> I
> >>> believe one could use this (after squaring) as the estimate of the
> timing
> >>> variance, correct?
> >>>
> >>> I believe depending on how we model the VCO, we can get away with a
> >> scalar
> >>> Kalman filter and circumvent the matrix and vector operations.
> >>> I tried to simulate it in Matlab, and it kind of worked, but I noticed
> >> some
> >>> strange effects.
> >>>
> >>> a) I made a very simple VCO model, that simulates the phase error. It
> is
> >>> x[k+1] = x[k] + KVCO * u with u being the DAC code. If KVCO is chosen
> >>> correctly, this perfectly models the phase measurements. I assumed the
> >>> process noise is zero, and the 1PPS jitter contributes only to the
> >>> measurement noise.
> >>>
> >>> b) from the above model, we have a very simple state space model, if
> you
> >>> want to call it like this. We have A = 1, B=KVCO, C=1, D=0.
> >>>
> >>> c) in the "prediction phase" for the Kalman filter, the error
> covariance
> >>> (in this case, the error variance, actually) is P_new=APA' + Q, which
> >>> reduces in this case to P_new=P+Q with Q being the process noise
> >> variance,
> >>> which I believe is negligible in this case.
> >>>
> >>> d) in the "update phase" of the Kalman filter, we find the Kalman gain
> as
> >>> K=P*C*inv(C'*P*C + R), and this reduces, as everything is scalar and
> C=1,
> >>> to K=P/(P+R), with R being the measurement noise, which, I believe, is
> >>> equal to the timing accuracy estimate of the GPS module. Correct? we
> then
> >>> update the model xhat = A*xhat + b*u + K*(y-c*xhat), which simplifies
> to
> >>> xhat=xhat + Kvco*u + K*(y-xhat). Nothing special so far.
> >>>
> >>> e) now comes my weird observation. I don't know whether this is
> correct.
> >>> The error covariance is now updated according to P=(I-K*C)*P, this
> breaks
> >>> down to P=P-K*P. I now observe that P behaves very odd, first we set
> >> P=P+Q,
> >>> and then we set P=P-K*P. It does not really converge in my Matlab
> >>> Simulation, and I see that the noise is filtered somewhat, but not very
> >>> good. It could also be related to my variances not being correctly
> set, I
> >>> am not sure. Or I made some mistakes with the equations.
> >>>
> >>> Any hints?
> >>>
> >>> As soon as I see it sort of working in Matlab, I want to test it on my
> >>> GPSDO. Especially the fact that I have an estimate of the timing error
> >> (the
> >>> GPS module announces this value via a special telegram!) I find very
> >>> amazing and hope I can make use of this.
> >>>
> >>> best
> >>> Tobias
> >>> HB9FSX
> >>>
> >>>
> >>>
> >>> On Tue, Apr 12, 2022 at 11:23 PM glen english LIST <
> >> glenlist at cortexrf.com.au>
> >>> wrote:
> >>>
> >>>> For isolating noise, (for the purpose of off line analysis)  , using
> ICA
> >>>> (Independent Component Analysis) , a kind of blind source separation ,
> >>>> can assist parting out the various noises to assist understanding the
> >>>> system better . There are some Python primers around for it.
> >>>>
> >>>> fantastic discussion going on here. love it.
> >>>>
> >>>> glen
> >>>>
> >>>>
> >>>> On 12/04/2022 6:42 pm, Markus Kleinhenz via time-nuts wrote:
> >>>>> Hi Tobias,
> >>>>>
> >>>>> Am 11.04.2022 um 13:33 schrieb Pluess, Tobias via tim
> >>>> _______________________________________________
> >>>> time-nuts mailing list -- time-nuts at lists.febo.com -- To unsubscribe
> >> send
> >>>> an email to time-nuts-leave at lists.febo.com
> >>>> To unsubscribe, go to and follow the instructions there.
> >>> _______________________________________________
> >>> time-nuts mailing list -- time-nuts at lists.febo.com -- To unsubscribe
> >> send an email to time-nuts-leave at lists.febo.com
> >>> To unsubscribe, go to and follow the instructions there.
> >> _______________________________________________
> >> time-nuts mailing list -- time-nuts at lists.febo.com -- To unsubscribe
> send
> >> an email to time-nuts-leave at lists.febo.com
> >> To unsubscribe, go to and follow the instructions there.
> > _______________________________________________
> > time-nuts mailing list -- time-nuts at lists.febo.com -- To unsubscribe
> send an email to time-nuts-leave at lists.febo.com
> > To unsubscribe, go to and follow the instructions there.
> _______________________________________________
> time-nuts mailing list -- time-nuts at lists.febo.com -- To unsubscribe send
> an email to time-nuts-leave at lists.febo.com
> To unsubscribe, go to and follow the instructions there.
>