summaryrefslogtreecommitdiffstats
path: root/src/pulsecore/time-smoother.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/pulsecore/time-smoother.c')
-rw-r--r--src/pulsecore/time-smoother.c378
1 files changed, 378 insertions, 0 deletions
diff --git a/src/pulsecore/time-smoother.c b/src/pulsecore/time-smoother.c
new file mode 100644
index 00000000..6bda3df0
--- /dev/null
+++ b/src/pulsecore/time-smoother.c
@@ -0,0 +1,378 @@
+/* $Id$ */
+
+/***
+ This file is part of PulseAudio.
+
+ Copyright 2007 Lennart Poettering
+
+ PulseAudio is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as
+ published by the Free Software Foundation; either version 2.1 of the
+ License, or (at your option) any later version.
+
+ PulseAudio is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with PulseAudio; if not, write to the Free Software
+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+ USA.
+***/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+
+#include <pulse/sample.h>
+#include <pulse/xmalloc.h>
+
+#include <pulsecore/macro.h>
+
+#include "time-smoother.h"
+
+#define HISTORY_MAX 50
+
+/*
+ * Implementation of a time smoothing algorithm to synchronize remote
+ * clocks to a local one. Evens out noise, adjusts to clock skew and
+ * allows cheap estimations of the remote time while clock updates may
+ * be seldom and recieved in non-equidistant intervals.
+ *
+ * Basically, we estimate the gradient of received clock samples in a
+ * certain history window (of size 'history_time') with linear
+ * regression. With that info we estimate the remote time in
+ * 'adjust_time' ahead and smoothen our current estimation function
+ * towards that point with a 3rd order polynomial interpolation with
+ * fitting derivatives. (more or less a b-spline)
+ *
+ * The larger 'history_time' is chosen the better we will surpress
+ * noise -- but we'll adjust to clock skew slower..
+ *
+ * The larger 'adjust_time' is chosen the smoother our estimation
+ * function will be -- but we'll adjust to clock skew slower, too.
+ *
+ * If 'monotonic' is TRUE the resulting estimation function is
+ * guaranteed to be monotonic.
+ */
+
+struct pa_smoother {
+ pa_usec_t adjust_time, history_time;
+ pa_bool_t monotonic;
+
+ pa_usec_t time_offset;
+
+ pa_usec_t px, py; /* Point p, where we want to reach stability */
+ double dp; /* Gradient we want at point p */
+
+ pa_usec_t ex, ey; /* Point e, which we estimated before and need to smooth to */
+ double de; /* Gradient we estimated for point e */
+
+ /* History of last measurements */
+ pa_usec_t history_x[HISTORY_MAX], history_y[HISTORY_MAX];
+ unsigned history_idx, n_history;
+
+ /* To even out for monotonicity */
+ pa_usec_t last_y;
+
+ /* Cached parameters for our interpolation polynomial y=ax^3+b^2+cx */
+ double a, b, c;
+ pa_bool_t abc_valid;
+
+ pa_bool_t paused;
+ pa_usec_t pause_time;
+};
+
+pa_smoother* pa_smoother_new(pa_usec_t adjust_time, pa_usec_t history_time, pa_bool_t monotonic) {
+ pa_smoother *s;
+
+ pa_assert(adjust_time > 0);
+ pa_assert(history_time > 0);
+
+ s = pa_xnew(pa_smoother, 1);
+ s->adjust_time = adjust_time;
+ s->history_time = history_time;
+ s->time_offset = 0;
+ s->monotonic = monotonic;
+
+ s->px = s->py = 0;
+ s->dp = 1;
+
+ s->ex = s->ey = 0;
+ s->de = 1;
+
+ s->history_idx = 0;
+ s->n_history = 0;
+
+ s->last_y = 0;
+
+ s->abc_valid = FALSE;
+
+ s->paused = FALSE;
+
+ return s;
+}
+
+void pa_smoother_free(pa_smoother* s) {
+ pa_assert(s);
+
+ pa_xfree(s);
+}
+
+static void drop_old(pa_smoother *s, pa_usec_t x) {
+ unsigned j;
+
+ /* First drop items from history which are too old, but make sure
+ * to always keep two entries in the history */
+
+ for (j = s->n_history; j > 2; j--) {
+
+ if (s->history_x[s->history_idx] + s->history_time >= x) {
+ /* This item is still valid, and thus all following ones
+ * are too, so let's quit this loop */
+ break;
+ }
+
+ /* Item is too old, let's drop it */
+ s->history_idx ++;
+ while (s->history_idx >= HISTORY_MAX)
+ s->history_idx -= HISTORY_MAX;
+
+ s->n_history --;
+ }
+}
+
+static void add_to_history(pa_smoother *s, pa_usec_t x, pa_usec_t y) {
+ unsigned j;
+ pa_assert(s);
+
+ drop_old(s, x);
+
+ /* Calculate position for new entry */
+ j = s->history_idx + s->n_history;
+ while (j >= HISTORY_MAX)
+ j -= HISTORY_MAX;
+
+ /* Fill in entry */
+ s->history_x[j] = x;
+ s->history_y[j] = y;
+
+ /* Adjust counter */
+ s->n_history ++;
+
+ /* And make sure we don't store more entries than fit in */
+ if (s->n_history >= HISTORY_MAX) {
+ s->history_idx += s->n_history - HISTORY_MAX;
+ s->n_history = HISTORY_MAX;
+ }
+}
+
+static double avg_gradient(pa_smoother *s, pa_usec_t x) {
+ unsigned i, j, c = 0;
+ int64_t ax = 0, ay = 0, k, t;
+ double r;
+
+ drop_old(s, x);
+
+ /* First, calculate average of all measurements */
+ i = s->history_idx;
+ for (j = s->n_history; j > 0; j--) {
+
+ ax += s->history_x[i];
+ ay += s->history_y[i];
+ c++;
+
+ i++;
+ while (i >= HISTORY_MAX)
+ i -= HISTORY_MAX;
+ }
+
+ /* Too few measurements, assume gradient of 1 */
+ if (c < 2)
+ return 1;
+
+ ax /= c;
+ ay /= c;
+
+ /* Now, do linear regression */
+ k = t = 0;
+
+ i = s->history_idx;
+ for (j = s->n_history; j > 0; j--) {
+ int64_t dx, dy;
+
+ dx = (int64_t) s->history_x[i] - ax;
+ dy = (int64_t) s->history_y[i] - ay;
+
+ k += dx*dy;
+ t += dx*dx;
+
+ i++;
+ while (i >= HISTORY_MAX)
+ i -= HISTORY_MAX;
+ }
+
+ r = (double) k / t;
+
+ return s->monotonic && r < 0 ? 0 : r;
+}
+
+static void estimate(pa_smoother *s, pa_usec_t x, pa_usec_t *y, double *deriv) {
+ pa_assert(s);
+ pa_assert(y);
+
+ if (x >= s->px) {
+ int64_t t;
+
+ /* The requested point is right of the point where we wanted
+ * to be on track again, thus just linearly estimate */
+
+ t = (int64_t) s->py + (int64_t) (s->dp * (x - s->px));
+
+ if (t < 0)
+ t = 0;
+
+ *y = (pa_usec_t) t;
+
+ if (deriv)
+ *deriv = s->dp;
+
+ } else {
+
+ if (!s->abc_valid) {
+ pa_usec_t ex, ey, px, py;
+ int64_t kx, ky;
+ double de, dp;
+
+ /* Ok, we're not yet on track, thus let's interpolate, and
+ * make sure that the first derivative is smooth */
+
+ /* We have two points: (ex|ey) and (px|py) with two gradients
+ * at these points de and dp. We do a polynomial interpolation
+ * of degree 3 with these 6 values */
+
+ ex = s->ex; ey = s->ey;
+ px = s->px; py = s->py;
+ de = s->de; dp = s->dp;
+
+ pa_assert(ex < px);
+
+ /* To increase the dynamic range and symplify calculation, we
+ * move these values to the origin */
+ kx = (int64_t) px - (int64_t) ex;
+ ky = (int64_t) py - (int64_t) ey;
+
+ /* Calculate a, b, c for y=ax^3+b^2+cx */
+ s->c = de;
+ s->b = (((double) (3*ky)/kx - dp - 2*de)) / kx;
+ s->a = (dp/kx - 2*s->b - de/kx) / (3*kx);
+
+ s->abc_valid = TRUE;
+ }
+
+ /* Move to origin */
+ x -= s->ex;
+
+ /* Horner scheme */
+ *y = (pa_usec_t) ((double) x * (s->c + (double) x * (s->b + (double) x * s->a)));
+
+ /* Move back from origin */
+ *y += s->ey;
+
+ /* Horner scheme */
+ if (deriv)
+ *deriv = s->c + ((double) x * (s->b*2 + (double) x * s->a*3));
+ }
+
+ /* Guarantee monotonicity */
+ if (s->monotonic) {
+
+ if (*y < s->last_y)
+ *y = s->last_y;
+ else
+ s->last_y = *y;
+
+ if (deriv && *deriv < 0)
+ *deriv = 0;
+ }
+}
+
+void pa_smoother_put(pa_smoother *s, pa_usec_t x, pa_usec_t y) {
+ pa_usec_t ney;
+ double nde;
+
+ pa_assert(s);
+ pa_assert(x >= s->time_offset);
+
+ /* Fix up x value */
+ if (s->paused)
+ x = s->pause_time;
+ else
+ x -= s->time_offset;
+
+ pa_assert(x >= s->ex);
+
+ /* First, we calculate the position we'd estimate for x, so that
+ * we can adjust our position smoothly from this one */
+ estimate(s, x, &ney, &nde);
+ s->ex = x; s->ey = ney; s->de = nde;
+
+ /* Then, we add the new measurement to our history */
+ add_to_history(s, x, y);
+
+ /* And determine the average gradient of the history */
+ s->dp = avg_gradient(s, x);
+
+ /* And calculate when we want to be on track again */
+ s->px = x + s->adjust_time;
+ s->py = y + s->dp *s->adjust_time;
+
+ s->abc_valid = FALSE;
+}
+
+pa_usec_t pa_smoother_get(pa_smoother *s, pa_usec_t x) {
+ pa_usec_t y;
+
+ pa_assert(s);
+ pa_assert(x >= s->time_offset);
+
+ /* Fix up x value */
+ if (s->paused)
+ x = s->pause_time;
+ else
+ x -= s->time_offset;
+
+ pa_assert(x >= s->ex);
+
+ estimate(s, x, &y, NULL);
+ return y;
+}
+
+void pa_smoother_set_time_offset(pa_smoother *s, pa_usec_t offset) {
+ pa_assert(s);
+
+ s->time_offset = offset;
+}
+
+void pa_smoother_pause(pa_smoother *s, pa_usec_t x) {
+ pa_assert(s);
+
+ if (s->paused)
+ return;
+
+ s->paused = TRUE;
+ s->pause_time = x;
+}
+
+void pa_smoother_resume(pa_smoother *s, pa_usec_t x) {
+ pa_assert(s);
+
+ if (!s->paused)
+ return;
+
+ s->paused = FALSE;
+ s->time_offset += x - s->pause_time;
+}