From a67c21f093202f142438689d3f7cfbdf4ea82eea Mon Sep 17 00:00:00 2001 From: Lennart Poettering Date: Sun, 28 Oct 2007 19:13:50 +0000 Subject: merge 'lennart' branch back into trunk. git-svn-id: file:///home/lennart/svn/public/pulseaudio/trunk@1971 fefdeb5f-60dc-0310-8127-8f9354f1896f --- src/pulsecore/time-smoother.c | 378 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 378 insertions(+) create mode 100644 src/pulsecore/time-smoother.c (limited to 'src/pulsecore/time-smoother.c') 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 +#endif + +#include + +#include +#include + +#include + +#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; +} -- cgit