Thread: [Linuxptp-devel] [PATCH RFC] Replace moving average with moving median for delay filtering.

SourceForge Headquarters 1320 Columbia Street Suite 310 San Diego, CA 92101 +1 (858) 422-6466

The moving average used in the path delay calculation is very sensitive
to outliers. Replace it with moving median.

This, for instance, allows much faster recovery from an external clock
time step which happened between receiving sync message and sending
delay_req message. The measured delay includes a large error, but the
median is still a good estimate of the delay and the first step
correction applied by the servo is right.

In this implementation the median update has linear time complexity.

Signed-off-by: Miroslav Lichvar <mli...@re...>
---
 clock.c             | 18 ++++++++--------
 makefile            |  2 +-
 mave.c => mmedian.c | 60 +++++++++++++++++++++++++++++++++++++++--------------
 mave.h => mmedian.h | 14 ++++++-------
 port.c              | 16 +++++++-------
 5 files changed, 70 insertions(+), 40 deletions(-)
 rename mave.c => mmedian.c (52%)
 rename mave.h => mmedian.h (77%)

diff --git a/clock.c b/clock.c
index 419a6b7..8e7fc56 100644
--- a/clock.c
+++ b/clock.c
@@ -26,7 +26,7 @@
 #include "clock.h"
 #include "clockadj.h"
 #include "foreign.h"
-#include "mave.h"
+#include "mmedian.h"
 #include "missing.h"
 #include "msg.h"
 #include "phc.h"
@@ -40,7 +40,7 @@
 
 #define CLK_N_PORTS (MAX_PORTS + 1) /* plus one for the UDS interface */
 #define N_CLOCK_PFD (N_POLLFD + 1) /* one extra per port, for the fault timer */
-#define MAVE_LENGTH 10
+#define MMEDIAN_LENGTH 10
 #define POW2_41 ((double)(1ULL << 41))
 
 #define ARRAY_SIZE(x) (sizeof(x) / sizeof((x)[0]))
@@ -85,7 +85,7 @@ struct clock {
 	enum servo_state servo_state;
 	tmv_t master_offset;
 	tmv_t path_delay;
-	struct mave *avg_delay;
+	struct mmedian *med_delay;
 	struct freq_estimator fest;
 	struct time_status_np status;
 	double nrr;
@@ -119,7 +119,7 @@ void clock_destroy(struct clock *c)
 		phc_close(c->clkid);
 	}
 	servo_destroy(c->servo);
-	mave_destroy(c->avg_delay);
+	mmedian_destroy(c->med_delay);
 	stats_destroy(c->stats.offset);
 	stats_destroy(c->stats.freq);
 	stats_destroy(c->stats.delay);
@@ -630,9 +630,9 @@ struct clock *clock_create(int phc_index, struct interface *iface, int count,
 		return NULL;
 	}
 	c->servo_state = SERVO_UNLOCKED;
-	c->avg_delay = mave_create(MAVE_LENGTH);
-	if (!c->avg_delay) {
-		pr_err("Failed to create moving average");
+	c->med_delay = mmedian_create(MMEDIAN_LENGTH);
+	if (!c->med_delay) {
+		pr_err("Failed to create moving median");
 		return NULL;
 	}
 	c->stats_interval = dds->stats_interval;
@@ -983,7 +983,7 @@ void clock_path_delay(struct clock *c, struct timespec req, struct timestamp rx,
 		pr_warning("c3 %10lld", c3);
 	}
 
-	c->path_delay = mave_accumulate(c->avg_delay, pd);
+	c->path_delay = mmedian_accumulate(c->med_delay, pd);
 
 	c->cur.meanPathDelay = tmv_to_TimeInterval(c->path_delay);
 
@@ -1154,7 +1154,7 @@ static void handle_state_decision_event(struct clock *c)
 
 	if (!cid_eq(&best_id, &c->best_id)) {
 		clock_freq_est_reset(c);
-		mave_reset(c->avg_delay);
+		mmedian_reset(c->med_delay);
 		c->path_delay = 0;
 		fresh_best = 1;
 	}
diff --git a/makefile b/makefile
index d79a1df..71119a1 100644
--- a/makefile
+++ b/makefile
@@ -32,7 +32,7 @@ VER     = -DVER=$(version)
 CFLAGS	= -Wall $(VER) $(INC) $(DEBUG) $(FEAT_CFLAGS) $(EXTRA_CFLAGS)
 LDLIBS	= -lm -lrt $(EXTRA_LDFLAGS)
 PRG	= ptp4l pmc phc2sys hwstamp_ctl
-OBJ     = bmc.o clock.o clockadj.o config.o fault.o fsm.o ptp4l.o mave.o \
+OBJ     = bmc.o clock.o clockadj.o config.o fault.o fsm.o ptp4l.o mmedian.o \
  msg.o phc.o pi.o port.o print.o raw.o servo.o sk.o stats.o tlv.o tmtab.o \
  transport.o udp.o udp6.o uds.o util.o version.o
 
diff --git a/mave.c b/mmedian.c
similarity index 52%
rename from mave.c
rename to mmedian.c
index 7742f05..f55dad0 100644
--- a/mave.c
+++ b/mmedian.c
@@ -1,5 +1,5 @@
 /**
- * @file mave.c
+ * @file mmedian.c
  * @note Copyright (C) 2011 Richard Cochran <ric...@gm...>
  *
  * This program is free software; you can redistribute it and/or modify
@@ -19,25 +19,33 @@
 #include <stdlib.h>
 #include <string.h>
 
-#include "mave.h"
+#include "mmedian.h"
 
-struct mave {
+struct mmedian {
 	int cnt;
 	int len;
 	int index;
-	tmv_t sum;
+	/* Indices sorted by value. */
+	int *order;
+	/* Values stored in circular buffer. */
 	tmv_t *val;
 };
 
-struct mave *mave_create(int length)
+struct mmedian *mmedian_create(int length)
 {
-	struct mave *m;
+	struct mmedian *m;
 	m = calloc(1, sizeof(*m));
-	if (!m) {
+	if (!m || length < 1) {
+		return NULL;
+	}
+	m->order = calloc(1, length * sizeof(*m->order));
+	if (!m->order) {
+		free(m);
 		return NULL;
 	}
 	m->val = calloc(1, length * sizeof(*m->val));
 	if (!m->val) {
+		free(m->order);
 		free(m);
 		return NULL;
 	}
@@ -45,28 +53,50 @@ struct mave *mave_create(int length)
 	return m;
 }
 
-void mave_destroy(struct mave *m)
+void mmedian_destroy(struct mmedian *m)
 {
+	free(m->order);
 	free(m->val);
 	free(m);
 }
 
-tmv_t mave_accumulate(struct mave *m, tmv_t val)
+tmv_t mmedian_accumulate(struct mmedian *m, tmv_t val)
 {
-	m->sum = tmv_sub(m->sum, m->val[m->index]);
+	int i;
+
 	m->val[m->index] = val;
-	m->index = (1 + m->index) % m->len;
-	m->sum = tmv_add(m->sum, val);
 	if (m->cnt < m->len) {
 		m->cnt++;
+	} else {
+		/* Remove index of the replaced value from order. */
+		for (i = 0; i < m->cnt; i++)
+			if (m->order[i] == m->index)
+				break;
+		for (; i + 1 < m->cnt; i++)
+			m->order[i] = m->order[i + 1];
 	}
-	return tmv_div(m->sum, m->cnt);
+
+	/* Insert index of the new value to order. */
+	for (i = m->cnt - 1; i > 0; i--) {
+		if (m->val[m->order[i - 1]] <= m->val[m->index])
+			break;
+		m->order[i] = m->order[i - 1];
+	}
+	m->order[i] = m->index;
+
+	m->index = (1 + m->index) % m->len;
+
+	if (m->cnt % 2)
+		return m->val[m->order[m->cnt / 2]];
+	else
+		return tmv_div(tmv_add(m->val[m->order[m->cnt / 2 - 1]],
+				       m->val[m->order[m->cnt / 2]]), 2);
 }
 
-void mave_reset(struct mave *m)
+void mmedian_reset(struct mmedian *m)
 {
 	m->cnt = 0;
 	m->index = 0;
-	m->sum = 0;
+	memset(m->order, 0, m->len * sizeof(*m->order));
 	memset(m->val, 0, m->len * sizeof(*m->val));
 }
diff --git a/mave.h b/mmedian.h
similarity index 77%
rename from mave.h
rename to mmedian.h
index 84241d4..f0bd8ac 100644
--- a/mave.h
+++ b/mmedian.h
@@ -1,6 +1,6 @@
 /**
- * @file mave.h
- * @brief Implements a moving average.
+ * @file mmedian.h
+ * @brief Implements a moving median.
  * @note Copyright (C) 2011 Richard Cochran <ric...@gm...>
  *
  * This program is free software; you can redistribute it and/or modify
@@ -22,14 +22,14 @@
 
 #include "tmv.h"
 
-struct mave;
+struct mmedian;
 
-struct mave *mave_create(int length);
+struct mmedian *mmedian_create(int length);
 
-void mave_destroy(struct mave *m);
+void mmedian_destroy(struct mmedian *m);
 
-tmv_t mave_accumulate(struct mave *m, tmv_t val);
+tmv_t mmedian_accumulate(struct mmedian *m, tmv_t val);
 
-void mave_reset(struct mave *m);
+void mmedian_reset(struct mmedian *m);
 
 #endif
diff --git a/port.c b/port.c
index 3fc262a..e6ef5b1 100644
--- a/port.c
+++ b/port.c
@@ -25,7 +25,7 @@
 
 #include "bmc.h"
 #include "clock.h"
-#include "mave.h"
+#include "mmedian.h"
 #include "missing.h"
 #include "msg.h"
 #include "port.h"
@@ -37,7 +37,7 @@
 #include "util.h"
 
 #define ALLOWED_LOST_RESPONSES 3
-#define PORT_MAVE_LENGTH 10
+#define PORT_MMEDIAN_LENGTH 10
 
 enum syfu_state {
 	SF_EMPTY,
@@ -82,7 +82,7 @@ struct port {
 	} seqnum;
 	struct tmtab tmtab;
 	tmv_t peer_delay;
-	struct mave *avg_delay;
+	struct mmedian *med_delay;
 	int log_sync_interval;
 	struct nrate_estimator nrate;
 	unsigned int pdr_missing;
@@ -1712,7 +1712,7 @@ calc:
 	pd = tmv_sub(pd, c2);
 	pd = tmv_div(pd, 2);
 
-	p->peer_delay = mave_accumulate(p->avg_delay, pd);
+	p->peer_delay = mmedian_accumulate(p->med_delay, pd);
 
 	p->peerMeanPathDelay = tmv_to_TimeInterval(p->peer_delay);
 
@@ -1841,7 +1841,7 @@ void port_close(struct port *p)
 		port_disable(p);
 	}
 	transport_destroy(p->trp);
-	mave_destroy(p->avg_delay);
+	mmedian_destroy(p->med_delay);
 	free(p);
 }
 
@@ -2299,9 +2299,9 @@ struct port *port_open(int phc_index,
 	p->delayMechanism = interface->dm;
 	p->versionNumber = PTP_VERSION;
 
-	p->avg_delay = mave_create(PORT_MAVE_LENGTH);
-	if (!p->avg_delay) {
-		pr_err("Failed to create moving average");
+	p->med_delay = mmedian_create(PORT_MMEDIAN_LENGTH);
+	if (!p->med_delay) {
+		pr_err("Failed to create moving median");
 		transport_destroy(p->trp);
 		free(p);
 		return NULL;
-- 
1.8.3.1





Thread: [Linuxptp-devel] [PATCH RFC] Replace moving average with moving median for delay filtering.

PTP IEEE 1588 stack for Linux

linuxptp-devel