[time-nuts] Re: GPSDO Control loop autotuning

[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.