[time-nuts] Re: GPSDO Control loop autotuning

[Markus Kleinhenz]

Erik,

https://filterpy.readthedocs.io/en/latest/kalman/KalmanFilter.html

has the documentation for a python implementation of kalman filters and
some good introduction.

Markus.

Am 19.04.2022 um 21:03 schrieb Erik Kaashoek:
> 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.