Include offset correction error in dispersion

This commit is contained in:
Miroslav Lichvar
2010-02-10 16:05:13 +01:00
parent 20d898d182
commit 8cb6fcad7e
7 changed files with 117 additions and 78 deletions

View File

@@ -43,6 +43,7 @@ extern RefclockDriver RCL_PPS_driver;
struct FilterSample {
double offset;
double dispersion;
struct timeval sample_time;
};
@@ -52,6 +53,7 @@ struct MedianFilter {
int used;
int last;
struct FilterSample *samples;
int *sort_array;
};
struct RCL_Instance_Record {
@@ -89,13 +91,13 @@ static int pps_stratum(RCL_Instance instance, struct timeval *tv);
static void poll_timeout(void *arg);
static void slew_samples(struct timeval *raw, struct timeval *cooked, double dfreq, double afreq,
double doffset, int is_step_change, void *anything);
static void log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset);
static void log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset, double dispersion);
static void filter_init(struct MedianFilter *filter, int length);
static void filter_fini(struct MedianFilter *filter);
static void filter_reset(struct MedianFilter *filter);
static void filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset);
static int filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset);
static void filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset, double dispersion);
static int filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion);
static int filter_get_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion);
static void filter_slew_samples(struct MedianFilter *filter, struct timeval *when, double dfreq, double doffset);
@@ -317,11 +319,12 @@ RCL_GetDriverOption(RCL_Instance instance, char *name)
int
RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset, NTP_Leap leap_status)
{
double correction, err;
double correction, dispersion;
struct timeval cooked_time;
LCL_GetOffsetCorrection(sample_time, &correction, &err);
LCL_GetOffsetCorrection(sample_time, &correction, &dispersion);
UTI_AddDoubleToTimeval(sample_time, correction, &cooked_time);
dispersion += LCL_GetSysPrecisionAsQuantum();
if (!valid_sample_time(instance, sample_time))
return 0;
@@ -331,10 +334,10 @@ RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset,
offset, offset - correction + instance->offset);
#endif
filter_add_sample(&instance->filter, &cooked_time, offset - correction + instance->offset);
filter_add_sample(&instance->filter, &cooked_time, offset - correction + instance->offset, dispersion);
instance->leap_status = leap_status;
log_sample(instance, &cooked_time, 0, offset, offset - correction + instance->offset);
log_sample(instance, &cooked_time, 0, offset, offset - correction + instance->offset, dispersion);
return 1;
}
@@ -342,18 +345,13 @@ RCL_AddSample(RCL_Instance instance, struct timeval *sample_time, double offset,
int
RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
{
double correction, err, offset;
double correction, dispersion, offset;
struct timeval cooked_time;
int rate;
struct timeval ref_time;
int is_synchronised, stratum;
double root_delay, root_dispersion, distance;
NTP_Leap leap;
unsigned long ref_id;
LCL_GetOffsetCorrection(pulse_time, &correction, &err);
LCL_GetOffsetCorrection(pulse_time, &correction, &dispersion);
UTI_AddDoubleToTimeval(pulse_time, correction, &cooked_time);
dispersion += LCL_GetSysPrecisionAsQuantum();
if (!valid_sample_time(instance, pulse_time))
return 0;
@@ -372,10 +370,10 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
if (instance->lock_ref != -1) {
struct timeval ref_sample_time;
double sample_diff, ref_offset, shift;
double sample_diff, ref_offset, ref_dispersion, shift;
if (!filter_get_last_sample(&refclocks[instance->lock_ref].filter,
&ref_sample_time, &ref_offset))
&ref_sample_time, &ref_offset, &ref_dispersion))
return 0;
UTI_DiffTimevalsToDouble(&sample_diff, &cooked_time, &ref_sample_time);
@@ -390,7 +388,7 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
offset += shift;
if (fabs(ref_offset - offset) >= 0.2 / rate)
if (fabs(ref_offset - offset) + ref_dispersion + dispersion >= 0.2 / rate)
return 0;
#if 0
@@ -398,6 +396,12 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
second, offset, ref_offset - offset, sample_diff);
#endif
} else {
struct timeval ref_time;
int is_synchronised, stratum;
double root_delay, root_dispersion, distance;
NTP_Leap leap;
unsigned long ref_id;
/* Ignore the pulse if we are not well synchronized */
REF_GetReferenceParams(&cooked_time, &is_synchronised, &leap, &stratum,
@@ -420,10 +424,10 @@ RCL_AddPulse(RCL_Instance instance, struct timeval *pulse_time, double second)
second, offset);
#endif
filter_add_sample(&instance->filter, &cooked_time, offset);
filter_add_sample(&instance->filter, &cooked_time, offset, dispersion);
instance->leap_status = LEAP_Normal;
log_sample(instance, &cooked_time, 1, second, offset);
log_sample(instance, &cooked_time, 1, second, offset, dispersion);
return 1;
}
@@ -545,7 +549,7 @@ slew_samples(struct timeval *raw, struct timeval *cooked, double dfreq, double a
}
static void
log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset)
log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double raw_offset, double cooked_offset, double dispersion)
{
char sync_stats[4] = {'N', '+', '-', '?'};
@@ -564,11 +568,11 @@ log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double
if (((logwrites++) % 32) == 0) {
fprintf(logfile,
"====================================================================\n"
" Date (UTC) Time Refid DP L P Raw offset Cooked offset\n"
"====================================================================\n");
"===============================================================================\n"
" Date (UTC) Time Refid DP L P Raw offset Cooked offset Disp.\n"
"===============================================================================\n");
}
fprintf(logfile, "%s.%06d %-5s %3d %1c %1d %13.6e %13.6e\n",
fprintf(logfile, "%s.%06d %-5s %3d %1c %1d %13.6e %13.6e %10.3e\n",
UTI_TimeToLogForm(sample_time->tv_sec),
(int)sample_time->tv_usec,
UTI_RefidToString(instance->ref_id),
@@ -576,7 +580,8 @@ log_sample(RCL_Instance instance, struct timeval *sample_time, int pulse, double
sync_stats[instance->leap_status],
pulse,
raw_offset,
cooked_offset);
cooked_offset,
dispersion);
fflush(logfile);
}
@@ -591,12 +596,14 @@ filter_init(struct MedianFilter *filter, int length)
filter->used = 0;
filter->last = -1;
filter->samples = MallocArray(struct FilterSample, filter->length);
filter->sort_array = MallocArray(int, filter->length);
}
static void
filter_fini(struct MedianFilter *filter)
{
Free(filter->samples);
Free(filter->sort_array);
}
static void
@@ -607,7 +614,7 @@ filter_reset(struct MedianFilter *filter)
}
static void
filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset)
filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, double offset, double dispersion)
{
filter->index++;
filter->index %= filter->length;
@@ -617,23 +624,30 @@ filter_add_sample(struct MedianFilter *filter, struct timeval *sample_time, doub
filter->samples[filter->index].sample_time = *sample_time;
filter->samples[filter->index].offset = offset;
filter->samples[filter->index].dispersion = dispersion;
}
static int
filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset)
filter_get_last_sample(struct MedianFilter *filter, struct timeval *sample_time, double *offset, double *dispersion)
{
if (filter->last < 0)
return 0;
*sample_time = filter->samples[filter->last].sample_time;
*offset = filter->samples[filter->last].offset;
*dispersion = filter->samples[filter->last].dispersion;
return 1;
}
static const struct FilterSample *tmp_sorted_array;
static int
sample_compare(const void *a, const void *b)
{
const struct FilterSample *s1 = a, *s2 = b;
const struct FilterSample *s1, *s2;
s1 = &tmp_sorted_array[*(int *)a];
s2 = &tmp_sorted_array[*(int *)b];
if (s1->offset < s2->offset)
return -1;
@@ -651,44 +665,66 @@ filter_get_sample(struct MedianFilter *filter, struct timeval *sample_time, doub
if (filter->used == 1) {
*sample_time = filter->samples[filter->index].sample_time;
*offset = filter->samples[filter->index].offset;
*dispersion = 0.0;
*dispersion = filter->samples[filter->index].dispersion;
} else {
int i, from, to;
double x, x1, y, d;
struct FilterSample *s;
int i, j, from, to;
double x, x1, y, d, e, min_dispersion;
/* sort samples by offset */
qsort(filter->samples, filter->used, sizeof (struct FilterSample), sample_compare);
/* average the half of the samples closest to the median */
if (filter->used > 2) {
from = (filter->used + 2) / 4;
to = filter->used - from;
} else {
from = 0;
to = filter->used;
/* find minimum dispersion */
for (i = 1, min_dispersion = filter->samples[0].dispersion; i < filter->used; i++) {
if (min_dispersion > filter->samples[i].dispersion)
min_dispersion = filter->samples[i].dispersion;
}
/* select samples with dispersion better than 1.5 * minimum */
for (i = j = 0; i < filter->used; i++) {
if (filter->samples[i].dispersion <= 1.5 * min_dispersion)
filter->sort_array[j++] = i;
}
assert(j > 0);
/* and sort their indexes by offset */
tmp_sorted_array = filter->samples;
qsort(filter->sort_array, j, sizeof (int), sample_compare);
/* select half of the samples closest to the median */
if (j > 2) {
from = (j + 2) / 4;
to = j - from;
} else {
from = 0;
to = j;
}
/* average offset and sample time */
for (i = from, x = y = 0.0; i < to; i++) {
s = &filter->samples[filter->sort_array[i]];
#if 0
LOG(LOGS_INFO, LOGF_Refclock, "refclock averaging offset %.9f [%s]",
filter->samples[i].offset, UTI_TimevalToString(&filter->samples[i].sample_time));
LOG(LOGS_INFO, LOGF_Refclock, "refclock averaging sample: offset %.9f dispersion %.9f [%s]",
s->offset, s->dispersion, UTI_TimevalToString(&filter->samples[i].sample_time));
#endif
UTI_DiffTimevalsToDouble(&x1, &filter->samples[i].sample_time, &filter->samples[0].sample_time);
UTI_DiffTimevalsToDouble(&x1, &s->sample_time, &filter->samples[0].sample_time);
x += x1;
y += filter->samples[i].offset;
y += s->offset;
}
x /= to - from;
y /= to - from;
for (i = from, d = 0.0; i < to; i++)
d += (filter->samples[i].offset - y) * (filter->samples[i].offset - y);
for (i = from, d = e = 0.0; i < to; i++) {
s = &filter->samples[filter->sort_array[i]];
d += (s->offset - y) * (s->offset - y);
e += s->dispersion;
}
d = sqrt(d / (to - from));
e /= to - from;
UTI_AddDoubleToTimeval(&filter->samples[0].sample_time, x, sample_time);
*offset = y;
*dispersion = d;
*dispersion = d + e;
}
return 1;