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 >> <glenl...@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@lists.febo.com -- To unsubscribe >>> send >>> an email to time-nuts-le...@lists.febo.com >>> To unsubscribe, go to and follow the instructions there. >> _______________________________________________ >> time-nuts mailing list -- time-nuts@lists.febo.com -- To unsubscribe >> send an email to time-nuts-le...@lists.febo.com >> To unsubscribe, go to and follow the instructions there. > _______________________________________________ > time-nuts mailing list -- time-nuts@lists.febo.com -- To unsubscribe > send an email to time-nuts-le...@lists.febo.com > To unsubscribe, go to and follow the instructions there. _______________________________________________ time-nuts mailing list -- time-nuts@lists.febo.com -- To unsubscribe send an email to time-nuts-le...@lists.febo.com To unsubscribe, go to and follow the instructions there.