[time-nuts] Re: GPSDO Control loop autotuning

[Carsten Andrich]

Hi Erik,

when I familiarized myself with the topic back in 2018, I found this paper:
https://ieeexplore.ieee.org/document/7324685 (IEEE paywall)
https://www.researchgate.net/publication/308730123_Synchronization_robustness_using_Kalman-based_clock_servos

They found that Kalman filters and PI controllers perform equally when 
used digital clock servos if tuned appropriately. While the paper is 
about network clock synchronization, the problem is very similar to OCXO 
disciplining. Instead of GNSS PPS jitter you have packet delay 
variation, which is more complex than simple jitter. I haven't checked 
any scientific literature after 2018.

Best regards,
Carsten

On 21.04.22 11:00, Erik Kaashoek wrote:
> 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.

-- 
M.Sc. Carsten Andrich

Technische Universität Ilmenau
Fachgebiet Elektronische Messtechnik und Signalverarbeitung (EMS)
Helmholtzplatz 2
98693 Ilmenau
T +49 3677 69-4269