Hardware Clock Drift Measurement with Python + NTP

README.md

Hardware Clock Drift Measurement with Python + NTP

I've been tinkering with microcontroller code and have learned about "real time clocks" and how computers attempt to keep accurate time.

More recently, I've become loosely familiar with NTP at the UDP level and realized that I could measure the clock drift over time just by repeatedly fetching hardware time and NTP time and comparing their differences over time.

I found that my desktop machine's RTC drifted by about 6 seconds in 33 hours. This kind of error is acceptable for some applications and unacceptable for others.

I did this on a Windows 10 machine, but I first

• Turned off the automatic system clock update. Presumably, this is periodically fetching NTP to stay accurate.
• Set the machie to never turn off

The scripts require numpy, matplotlib, and ntplib. To measure your own setup, simply run

python clock_drift_record.py

And you can visualize the data recorded to a .npy or .bin file by running

python clock_drift_plot.py clock_drift_0123456789ABCDEF.bin

Here's an example plot generated by these scripts

clock_drift_example.png

clock_drift_plot.py

import os, sys, time, uuid, matplotlib
from numpy import array, zeros, frombuffer, save, load, reshape, diff
from os.path import isfile, isdir, join
from matplotlib.pyplot import plot, xlabel, ylabel, title, grid, show, subplots
matplotlib.rc('font', family='sans-serif', weight='bold', size=22)

if len(sys.argv) < 2:
    print(f"Intended usage:\n  > python clock_drift_plot.py drift.npy") 
    sys.exit()
if not isfile(sys.argv[1]):
    print(f"First argument is not a file. First argument is\n  > {sys.argv[1]}")
    sys.exit()

fpath = sys.argv[1]
if fpath.endswith("npy"):
    a = load(fpath)
elif fpath.endswith("bin"):
    f = open(fpath, "rb")
    a = frombuffer(f.read(), dtype='float64')
    a = reshape(a, (int(len(a)/2), 2))

da = a - a[0, :] # a - first time stamp of a
dda = diff(da, axis=1)
t = a[:, 1] # Let's plot against NTP time on x-axis
t = t-t[0]

year70 = 2208988800
plot(array(t)/3600, dda, "o", linewidth=5, markersize=10)
xlabel("test time [hour]")
ylabel("error [sec]    (ntp - hw)")
title(f"RTC clock drift from\n{time.ctime(a[0,1]-year70)} to {time.ctime(a[-1, 1]-year70)}")
grid()
show()

clock_drift_record.py

import ntplib, uuid
import time
from numpy import ndarray, zeros, frombuffer, float64, concatenate, save

# make a prefix for saving files
prefix = f"clock_drift_{uuid.uuid1()}"
f=open(f"{prefix}.bin", "wb")

# make an array with two columns of data
a=zeros((0,2), dtype=float64)

n=0
c = ntplib.NTPClient()
start_time = time.time()
while True:
    try:
        # capture both times sequentially
        hwsec = float64(time.time())
        ntpr = c.request('time-a-g.nist.gov')
        n += 1

        # interpret ntp time to float64 seconds
        intsec = float64( frombuffer( ntpr.to_data()[40:44], dtype='uint32' ).byteswap()[0] )
        fracsec = float64(frombuffer( ntpr.to_data()[44:48], dtype='uint32' ).byteswap()[0])/(2**32)
        ntpsec = intsec + fracsec

        # save to file and flush so we don't lose results if something goes wrong
        f.write(hwsec.tobytes())
        f.write(ntpsec.tobytes())
        f.flush()

        arow = zeros((1,2), dtype='float64')
        arow[0,0] = hwsec
        arow[0,1] = ntpsec
        a = concatenate((a, arow))

        print(f"\nStarted @ {time.ctime(start_time)}")
        print(f" > {int(hwsec-start_time)} seconds later, took a measurement @ {time.ctime(hwsec)}")

        if a.shape[0] >= 2:
            hwdiff = a[-1, 0] - a[0,0]
            ntpdiff = a[-1, 1] - a[0,1]
            ntpdiff_hwdiff = ntpdiff - hwdiff
            print(f"TOTAL: hwdiff = {hwdiff}, ntpdiff = {ntpdiff}")
            print(f"TOTAL: ntpdiff - hwdiff = {ntpdiff_hwdiff}")
        time.sleep(60*10)
        print("")
    except KeyboardInterrupt:
        break
    except ntplib.NTPException:
        continue
# end while True

# save the array, a, to the file, f
save(f"{prefix}", a)
f.close()