# [ntp:questions] Re: Allan variance plots

David L. Mills mills at udel.edu
Wed Sep 17 18:37:51 UTC 2003

Piotr,

Here is a Matlab program to read loopstats files and construct Allan
deviation plots. You can modify it for another language or file format.
The key is in the first while loop.

Dave

%
% This program computes graphical statistics from
% loopstat files produced by NTP.
%
clear all;
fid = fopen('loopstats', 'r');
[table, count] = fscanf(fid, '%f %f %f %f %f %f %f');
day = table(1:7:count);
sec = table(2:7:count) + (day - day(1)) * 86400;
offset = table(3:7:count);
freq = table(4:7:count);
jitter = table(5:7:count);
wander = table(6:7:count);
poll = table(7:7:count);

y1 = offset; i = 1; d = 16;
while length(y1) >= 10
u = diff(y1) / d;
v = diff(u);
z1(i) = sqrt(mean(v .* v) / 2);
m1(i) = d;
y1 = y1(1:2:length(y1));
i = i + 1; d = d * 2;
end

loglog(m1, z1 * 1e6, '-k')
axis([1 1e5 1e-4 100]);
xlabel('Time Interval (s)');
ylabel('Allan Deviation (PPM)');
print -dtiff -r600 sim
%
% plot offset
%
clf reset; h = newplot;
set(h, 'FontSize', 16);
set(h, 'Position', [.12 .15 .85 .8]);
plot(sec / 3600, offset * 1e3, '-k');
%axis([0 24 -4 4]);
xlabel('Time (hr)');
ylabel('Offset (ms)');
print -deps offset
print -dtiff -r600 offset
%
xx = abs(offset) * 1e3;
fprintf('mean %.3f median %.3f max %.3f\n', mean(xx), median(xx),
max(xx))
%
% plot jitter
%
clf reset; h = newplot;
set(h, 'FontSize', 16);
set(h, 'Position', [.12 .15 .85 .8]);
plot(sec / 3600, jitter * 1e3, '-k');
%axis([0 24 -4 4]);
xlabel('Time (hr)');
ylabel('Jitter (ms)');
print -deps jitter
print -dtiff -r600 jitter
%
% plot wander
%
clf reset; h = newplot;
set(h, 'FontSize', 16);
set(h, 'Position', [.12 .15 .85 .8]);
plot(sec / 3600, wander, '-k');
%axis([0 24 -4 4]);
xlabel('Time (hr)');
ylabel('Wander (PPM)');
print -deps wander
print -dtiff -r600 wander
%
% plot frequency
%
clf reset; h = newplot;
set(h, 'FontSize', 16);
set(h, 'Position', [.12 .15 .85 .8]);
plot(sec / 3600, freq, '-k');
m = mean(freq);
%axis([0 24 m - .2 m + .2]);
xlabel('Time (hr)');
ylabel('Frequency (PPM)');
print -deps freq
print -dtiff -r600 freq
%
z = sort(abs(offset));
w = 1 - z / max(z);
%
clf reset; h = newplot;
set(h, 'FontSize', 16, 'Position', [.12 .15 .85 .8]);
semilogx(z * 1000, w, '-k');
xlabel('Offset (ms)');
ylabel('P[Offset < y]');
print -deps prob
print -dtiff -r600 prob

Piotr Trojanek wrote:
>
> Dear NTP experts,
>
> I have been reading NTP papers for a long time and now I wanted
> to make Allan variance plot of my machine's clock. I have collected
> time differences measured beetwen free running system clock and
> PPS signal from GPS receiver using modified ppsctl.c from PPSkit.
>
> My data is sampled once per second, so I have offsets like:
>
> .007540355
> .007513540
> .007486726
> .007460749
> .007433934
> .007412148
> .007380305
>
> Now I'm confused how to process this data to plots like those in
> NTP papers. I have now about samples from few days of non-stop operation.
>
> I feed this data directly to dataplot package (NIST software) and use
> 'AV PLOT' command in log-log scale -- result is straight line with
> slope +1, so now I think it is wrong method. I can't figure out the
> right way.
>
> Thanks in advance for any help!
>
> --
> Piotr Trojanek