[Linuxptp-devel] [PATCH] Fix drift calculation in PI servo with large values.

[Miroslav Lichvar]

When the drift value is adjusted by the newly measured frequency offset,
multiply the frequencies instead of adding the measured offset to the
old value to get accurate result even when updating a large drift.

Signed-off-by: Miroslav Lichvar <***@redhat.com>
---
pi.c | 6 ++++--
1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/pi.c b/pi.c
index bd78e40..ccad98e 100644
--- a/pi.c
+++ b/pi.c
@@ -107,8 +107,10 @@ static double pi_sample(struct servo *servo,
break;
}

- s->drift += (s->offset[1] - s->offset[0]) * 1e9 /
- (s->local[1] - s->local[0]);
+ /* Adjust drift by the measured frequency offset. */
+ s->drift = 1e9 - (1e9 - s->drift) *
+ (1.0 - (double)(s->offset[1] - s->offset[0]) /
+ (s->local[1] - s->local[0]));
if (s->drift < -s->maxppb)
s->drift = -s->maxppb;
else if (s->drift > s->maxppb)

[Richard Cochran]

Post by Miroslav Lichvar
When the drift value is adjusted by the newly measured frequency offset,
multiply the frequencies instead of adding the measured offset to the
old value to get accurate result even when updating a large drift.

In other words, when running with a programmed rate correction, it is
incorrect to add the rates, since the local time stamps are relative
to the adjusted time base. The local time values must be first scaled
to the original rate, right?

Thanks,
Richard

[Miroslav Lichvar]

Post by Richard Cochran

Post by Miroslav Lichvar
When the drift value is adjusted by the newly measured frequency offset,
multiply the frequencies instead of adding the measured offset to the
old value to get accurate result even when updating a large drift.

In other words, when running with a programmed rate correction, it is
incorrect to add the rates, since the local time stamps are relative
to the adjusted time base. The local time values must be first scaled
to the original rate, right?

Right.

For example, when a clock which doesn't need any frequency correction
is set to run at 200% (i.e. +100% frequency offset), the measured
drift will be 50% and the frequency of the clock needs to be set to
200% * (100% - 50%) instead of 200% - 50% to get back to the original
100% (0% offset).

The PI update in the locked state could use multiplication too to get
a response that doesn't depend on the absolute value of the frequency,
but it's probably not worth the extra computational cost.

[Richard Cochran]

Post by Miroslav Lichvar
+ /* Adjust drift by the measured frequency offset. */
+ s->drift = 1e9 - (1e9 - s->drift) *
+ (1.0 - (double)(s->offset[1] - s->offset[0]) /
+ (s->local[1] - s->local[0]));

I tried to understand this computation, but I could not see how you
derived it. I came up with the following (equivalent?) form.

s->drift += (1e9 - s->drift) * (s->offset[1] - s->offset[0]) /
(s->local[1] - s->local[0]);

(Notice the += assignment.)

I came up with this by starting with the simple drift calculation.

drift' = 1e9 * (s->offset[1] - s->offset[0]) /
(s->local[1] - s->local[0]);

If we convert the interval (s->local[1] - s->local[0]) into the
original time base, then the result will be the partial rate offset to
be added to the original offset. So we divide by the adjusted rate.

((s->local[1] - s->local[0]) / X)

where

X = (1 - s->drift / 1e9)

Since we divide by X in a denominator, we can move into the numerator.

drift' = 1e9 * (1 - s->drift / 1e9) * (s->offset[1] - s->offset[0]) /
(s->local[1] - s->local[0]);

Simplifying the first two terms and adding this to original drift
yields the first equation, above.

Thoughts?

Thanks,
Richard

[Miroslav Lichvar]

Post by Richard Cochran

Post by Miroslav Lichvar
+ /* Adjust drift by the measured frequency offset. */
+ s->drift = 1e9 - (1e9 - s->drift) *
+ (1.0 - (double)(s->offset[1] - s->offset[0]) /
+ (s->local[1] - s->local[0]));

I tried to understand this computation, but I could not see how you
derived it. I came up with the following (equivalent?) form.
s->drift += (1e9 - s->drift) * (s->offset[1] - s->offset[0]) /
(s->local[1] - s->local[0]);

Cool. That does look nicer than mine. I got it this way:

Original frequency of the clock:

old_freq = 1 - s->drift / 1e9

Ratio between the remote and local frequency:

freq_ratio = 1 - (s->offset[1] - s->offset[0]) /
(s->local[1] - s->local[0])
New frequency:

new_freq = old_freq * freq_ratio

New s->drift:

new_drift = 1e9 - 1e9 * new_freq

I put that together and simplified it a bit. But it seems I should
have gone further:

r = (s->offset[1] - s->offset[0]) / (s->local[1] - s->local[0])

new_drift = 1e9 - 1e9 * (1 - s->drift / 1e9) * (1 - r)
new_drift = 1e9 - (1e9 - s->drift) * (1 - r)
new_drift = 1e9 - (1e9 - 1e9 * r - s->drift + r * s->drift)
new_drift = 1e9 * r + s->drift - r * s->drift
new_drift = s->drift + r * (1e9 - s->drift)

And here is your form. :)

[Richard Cochran]

Post by Miroslav Lichvar
And here is your form. :)

Okay, thanks for that confirmation. I had a hunch that the two
expressions were equivalent. If you don't mind, I will commit your
patch, but with the simpler expression.

Thanks,
Richard

[Richard Cochran]

Post by Richard Cochran
Okay, thanks for that confirmation. I had a hunch that the two
expressions were equivalent. If you don't mind, I will commit your
patch, but with the simpler expression.

Like this?

---
commit 907635e62cb2774736a174fde603344f7958cd1b
Author: Miroslav Lichvar <***@redhat.com>
Date: Fri Dec 20 16:51:35 2013 +0100

Fix drift calculation in PI servo with large values.

When the drift value is adjusted by the newly measured frequency offset,
multiply the frequencies instead of adding the measured offset to the
old value to get accurate result even when updating a large drift.

Signed-off-by: Miroslav Lichvar <***@redhat.com>

diff --git a/pi.c b/pi.c
index bd78e40..689b232 100644
--- a/pi.c
+++ b/pi.c
@@ -107,8 +107,9 @@ static double pi_sample(struct servo *servo,
break;
}

- s->drift += (s->offset[1] - s->offset[0]) * 1e9 /
- (s->local[1] - s->local[0]);
+ /* Adjust drift by the measured frequency offset. */
+ s->drift += (1e9 - s->drift) * (s->offset[1] - s->offset[0]) /
+ (s->local[1] - s->local[0]);
if (s->drift < -s->maxppb)
s->drift = -s->maxppb;
else if (s->drift > s->maxppb)

[Miroslav Lichvar]

Post by Richard Cochran

Post by Richard Cochran
Okay, thanks for that confirmation. I had a hunch that the two
expressions were equivalent. If you don't mind, I will commit your
patch, but with the simpler expression.

Like this?

Ok, but I think you should get credit for finding that simpler form.
I'd suggest to put it in a separate commit or at least add a sign off
with comment to this commit.

Thanks.

Post by Richard Cochran
---
commit 907635e62cb2774736a174fde603344f7958cd1b
Date: Fri Dec 20 16:51:35 2013 +0100
Fix drift calculation in PI servo with large values.
When the drift value is adjusted by the newly measured frequency offset,
multiply the frequencies instead of adding the measured offset to the
old value to get accurate result even when updating a large drift.