]>
code.delx.au - pulseaudio/blob - src/pulsecore/time-smoother.c
2 This file is part of PulseAudio.
4 Copyright 2007 Lennart Poettering
6 PulseAudio is free software; you can redistribute it and/or modify
7 it under the terms of the GNU Lesser General Public License as
8 published by the Free Software Foundation; either version 2.1 of the
9 License, or (at your option) any later version.
11 PulseAudio is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
16 You should have received a copy of the GNU Lesser General Public
17 License along with PulseAudio; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
28 #include <pulse/sample.h>
29 #include <pulse/xmalloc.h>
31 #include <pulsecore/macro.h>
33 #include "time-smoother.h"
35 #define HISTORY_MAX 64
38 * Implementation of a time smoothing algorithm to synchronize remote
39 * clocks to a local one. Evens out noise, adjusts to clock skew and
40 * allows cheap estimations of the remote time while clock updates may
41 * be seldom and recieved in non-equidistant intervals.
43 * Basically, we estimate the gradient of received clock samples in a
44 * certain history window (of size 'history_time') with linear
45 * regression. With that info we estimate the remote time in
46 * 'adjust_time' ahead and smoothen our current estimation function
47 * towards that point with a 3rd order polynomial interpolation with
48 * fitting derivatives. (more or less a b-spline)
50 * The larger 'history_time' is chosen the better we will surpress
51 * noise -- but we'll adjust to clock skew slower..
53 * The larger 'adjust_time' is chosen the smoother our estimation
54 * function will be -- but we'll adjust to clock skew slower, too.
56 * If 'monotonic' is TRUE the resulting estimation function is
57 * guaranteed to be monotonic.
61 pa_usec_t adjust_time
, history_time
;
63 pa_usec_t time_offset
;
65 pa_usec_t px
, py
; /* Point p, where we want to reach stability */
66 double dp
; /* Gradient we want at point p */
68 pa_usec_t ex
, ey
; /* Point e, which we estimated before and need to smooth to */
69 double de
; /* Gradient we estimated for point e */
70 pa_usec_t ry
; /* The original y value for ex */
72 /* History of last measurements */
73 pa_usec_t history_x
[HISTORY_MAX
], history_y
[HISTORY_MAX
];
74 unsigned history_idx
, n_history
;
76 /* To even out for monotonicity */
77 pa_usec_t last_y
, last_x
;
79 /* Cached parameters for our interpolation polynomial y=ax^3+b^2+cx */
81 pa_bool_t abc_valid
:1;
83 pa_bool_t monotonic
:1;
85 pa_bool_t smoothing
:1; /* If FALSE we skip the polonyomial interpolation step */
92 pa_smoother
* pa_smoother_new(
93 pa_usec_t adjust_time
,
94 pa_usec_t history_time
,
98 pa_usec_t time_offset
,
103 pa_assert(adjust_time
> 0);
104 pa_assert(history_time
> 0);
105 pa_assert(min_history
>= 2);
106 pa_assert(min_history
<= HISTORY_MAX
);
108 s
= pa_xnew(pa_smoother
, 1);
109 s
->adjust_time
= adjust_time
;
110 s
->history_time
= history_time
;
111 s
->min_history
= min_history
;
112 s
->monotonic
= monotonic
;
113 s
->smoothing
= smoothing
;
115 pa_smoother_reset(s
, time_offset
, paused
);
120 void pa_smoother_free(pa_smoother
* s
) {
128 x = (x) % HISTORY_MAX; \
131 #define REDUCE_INC(x) \
133 x = ((x)+1) % HISTORY_MAX; \
137 static void drop_old(pa_smoother
*s
, pa_usec_t x
) {
139 /* Drop items from history which are too old, but make sure to
140 * always keep min_history in the history */
142 while (s
->n_history
> s
->min_history
) {
144 if (s
->history_x
[s
->history_idx
] + s
->history_time
>= x
)
145 /* This item is still valid, and thus all following ones
146 * are too, so let's quit this loop */
149 /* Item is too old, let's drop it */
150 REDUCE_INC(s
->history_idx
);
156 static void add_to_history(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
160 /* First try to update an existing history entry */
162 for (j
= s
->n_history
; j
> 0; j
--) {
164 if (s
->history_x
[i
] == x
) {
172 /* Drop old entries */
175 /* Calculate position for new entry */
176 j
= s
->history_idx
+ s
->n_history
;
186 /* And make sure we don't store more entries than fit in */
187 if (s
->n_history
> HISTORY_MAX
) {
188 s
->history_idx
+= s
->n_history
- HISTORY_MAX
;
189 REDUCE(s
->history_idx
);
190 s
->n_history
= HISTORY_MAX
;
194 static double avg_gradient(pa_smoother
*s
, pa_usec_t x
) {
195 unsigned i
, j
, c
= 0;
196 int64_t ax
= 0, ay
= 0, k
, t
;
199 /* Too few measurements, assume gradient of 1 */
200 if (s
->n_history
< s
->min_history
)
203 /* First, calculate average of all measurements */
205 for (j
= s
->n_history
; j
> 0; j
--) {
207 ax
+= (int64_t) s
->history_x
[i
];
208 ay
+= (int64_t) s
->history_y
[i
];
214 pa_assert(c
>= s
->min_history
);
218 /* Now, do linear regression */
222 for (j
= s
->n_history
; j
> 0; j
--) {
225 dx
= (int64_t) s
->history_x
[i
] - ax
;
226 dy
= (int64_t) s
->history_y
[i
] - ay
;
234 r
= (double) k
/ (double) t
;
236 return (s
->monotonic
&& r
< 0) ? 0 : r
;
239 static void calc_abc(pa_smoother
*s
) {
240 pa_usec_t ex
, ey
, px
, py
;
249 /* We have two points: (ex|ey) and (px|py) with two gradients at
250 * these points de and dp. We do a polynomial
251 * interpolation of degree 3 with these 6 values */
253 ex
= s
->ex
; ey
= s
->ey
;
254 px
= s
->px
; py
= s
->py
;
255 de
= s
->de
; dp
= s
->dp
;
259 /* To increase the dynamic range and symplify calculation, we
260 * move these values to the origin */
261 kx
= (int64_t) px
- (int64_t) ex
;
262 ky
= (int64_t) py
- (int64_t) ey
;
264 /* Calculate a, b, c for y=ax^3+bx^2+cx */
266 s
->b
= (((double) (3*ky
)/ (double) kx
- dp
- (double) (2*de
))) / (double) kx
;
267 s
->a
= (dp
/(double) kx
- 2*s
->b
- de
/(double) kx
) / (double) (3*kx
);
272 static void estimate(pa_smoother
*s
, pa_usec_t x
, pa_usec_t
*y
, double *deriv
) {
277 /* Linear interpolation right from px */
280 /* The requested point is right of the point where we wanted
281 * to be on track again, thus just linearly estimate */
283 t
= (int64_t) s
->py
+ (int64_t) llrint(s
->dp
* (double) (x
- s
->px
));
293 } else if (x
<= s
->ex
) {
294 /* Linear interpolation left from ex */
297 t
= (int64_t) s
->ey
- (int64_t) llrint(s
->de
* (double) (s
->ex
- x
));
308 /* Spline interpolation between ex and px */
311 /* Ok, we're not yet on track, thus let's interpolate, and
312 * make sure that the first derivative is smooth */
319 tx
-= (double) s
->ex
;
322 ty
= (tx
* (s
->c
+ tx
* (s
->b
+ tx
* s
->a
)));
324 /* Move back from origin */
325 ty
+= (double) s
->ey
;
327 *y
= ty
>= 0 ? (pa_usec_t
) llrint(ty
) : 0;
331 *deriv
= s
->c
+ (tx
* (s
->b
*2 + tx
* s
->a
*3));
334 /* Guarantee monotonicity */
337 if (deriv
&& *deriv
< 0)
342 void pa_smoother_put(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
353 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
358 /* First, we calculate the position we'd estimate for x, so that
359 * we can adjust our position smoothly from this one */
360 estimate(s
, x
, &ney
, &nde
);
361 s
->ex
= x
; s
->ey
= ney
; s
->de
= nde
;
365 /* Then, we add the new measurement to our history */
366 add_to_history(s
, x
, y
);
368 /* And determine the average gradient of the history */
369 s
->dp
= avg_gradient(s
, x
);
371 /* And calculate when we want to be on track again */
373 s
->px
= s
->ex
+ s
->adjust_time
;
374 s
->py
= s
->ry
+ (pa_usec_t
) llrint(s
->dp
* (double) s
->adjust_time
);
380 s
->abc_valid
= FALSE
;
383 pa_log_debug("%p, put(%llu | %llu) = %llu", s
, (unsigned long long) (x
+ s
->time_offset
), (unsigned long long) x
, (unsigned long long) y
);
387 pa_usec_t
pa_smoother_get(pa_smoother
*s
, pa_usec_t x
) {
396 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
402 estimate(s
, x
, &y
, NULL
);
406 /* Make sure the querier doesn't jump forth and back. */
416 pa_log_debug("%p, get(%llu | %llu) = %llu", s
, (unsigned long long) (x
+ s
->time_offset
), (unsigned long long) x
, (unsigned long long) y
);
422 void pa_smoother_set_time_offset(pa_smoother
*s
, pa_usec_t offset
) {
425 s
->time_offset
= offset
;
428 pa_log_debug("offset(%llu)", (unsigned long long) offset
);
432 void pa_smoother_pause(pa_smoother
*s
, pa_usec_t x
) {
439 pa_log_debug("pause(%llu)", (unsigned long long) x
);
446 void pa_smoother_resume(pa_smoother
*s
, pa_usec_t x
, pa_bool_t fix_now
) {
452 if (x
< s
->pause_time
)
456 pa_log_debug("resume(%llu)", (unsigned long long) x
);
460 s
->time_offset
+= x
- s
->pause_time
;
463 pa_smoother_fix_now(s
);
466 void pa_smoother_fix_now(pa_smoother
*s
) {
473 pa_usec_t
pa_smoother_translate(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y_delay
) {
483 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
485 estimate(s
, x
, &ney
, &nde
);
487 /* Play safe and take the larger gradient, so that we wakeup
488 * earlier when this is used for sleeping */
493 pa_log_debug("translate(%llu) = %llu (%0.2f)", (unsigned long long) y_delay
, (unsigned long long) ((double) y_delay
/ nde
), nde
);
496 return (pa_usec_t
) llrint((double) y_delay
/ nde
);
499 void pa_smoother_reset(pa_smoother
*s
, pa_usec_t time_offset
, pa_bool_t paused
) {
505 s
->ex
= s
->ey
= s
->ry
= 0;
511 s
->last_y
= s
->last_x
= 0;
513 s
->abc_valid
= FALSE
;
516 s
->time_offset
= s
->pause_time
= time_offset
;
519 pa_log_debug("reset()");