]>
code.delx.au - pulseaudio/blob - src/pulsecore/time-smoother.c
9b4be29fa811be5573822ee3e2038472312602f3
4 This file is part of PulseAudio.
6 Copyright 2007 Lennart Poettering
8 PulseAudio is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as
10 published by the Free Software Foundation; either version 2.1 of the
11 License, or (at your option) any later version.
13 PulseAudio is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
18 You should have received a copy of the GNU Lesser General Public
19 License along with PulseAudio; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
30 #include <pulse/sample.h>
31 #include <pulse/xmalloc.h>
33 #include <pulsecore/macro.h>
35 #include "time-smoother.h"
37 #define HISTORY_MAX 64
40 * Implementation of a time smoothing algorithm to synchronize remote
41 * clocks to a local one. Evens out noise, adjusts to clock skew and
42 * allows cheap estimations of the remote time while clock updates may
43 * be seldom and recieved in non-equidistant intervals.
45 * Basically, we estimate the gradient of received clock samples in a
46 * certain history window (of size 'history_time') with linear
47 * regression. With that info we estimate the remote time in
48 * 'adjust_time' ahead and smoothen our current estimation function
49 * towards that point with a 3rd order polynomial interpolation with
50 * fitting derivatives. (more or less a b-spline)
52 * The larger 'history_time' is chosen the better we will surpress
53 * noise -- but we'll adjust to clock skew slower..
55 * The larger 'adjust_time' is chosen the smoother our estimation
56 * function will be -- but we'll adjust to clock skew slower, too.
58 * If 'monotonic' is TRUE the resulting estimation function is
59 * guaranteed to be monotonic.
63 pa_usec_t adjust_time
, history_time
;
65 pa_usec_t time_offset
;
67 pa_usec_t px
, py
; /* Point p, where we want to reach stability */
68 double dp
; /* Gradient we want at point p */
70 pa_usec_t ex
, ey
; /* Point e, which we estimated before and need to smooth to */
71 double de
; /* Gradient we estimated for point e */
72 pa_usec_t ry
; /* The original y value for ex */
74 /* History of last measurements */
75 pa_usec_t history_x
[HISTORY_MAX
], history_y
[HISTORY_MAX
];
76 unsigned history_idx
, n_history
;
78 /* To even out for monotonicity */
79 pa_usec_t last_y
, last_x
;
81 /* Cached parameters for our interpolation polynomial y=ax^3+b^2+cx */
85 pa_bool_t monotonic
:1;
93 pa_smoother
* pa_smoother_new(pa_usec_t adjust_time
, pa_usec_t history_time
, pa_bool_t monotonic
, unsigned min_history
) {
96 pa_assert(adjust_time
> 0);
97 pa_assert(history_time
> 0);
98 pa_assert(min_history
>= 2);
99 pa_assert(min_history
<= HISTORY_MAX
);
101 s
= pa_xnew(pa_smoother
, 1);
102 s
->adjust_time
= adjust_time
;
103 s
->history_time
= history_time
;
105 s
->monotonic
= monotonic
;
110 s
->ex
= s
->ey
= s
->ry
= 0;
116 s
->last_y
= s
->last_x
= 0;
118 s
->abc_valid
= FALSE
;
122 s
->min_history
= min_history
;
127 void pa_smoother_free(pa_smoother
* s
) {
135 x = (x) % HISTORY_MAX; \
138 #define REDUCE_INC(x) \
140 x = ((x)+1) % HISTORY_MAX; \
144 static void drop_old(pa_smoother
*s
, pa_usec_t x
) {
146 /* Drop items from history which are too old, but make sure to
147 * always keep min_history in the history */
149 while (s
->n_history
> s
->min_history
) {
151 if (s
->history_x
[s
->history_idx
] + s
->history_time
>= x
)
152 /* This item is still valid, and thus all following ones
153 * are too, so let's quit this loop */
156 /* Item is too old, let's drop it */
157 REDUCE_INC(s
->history_idx
);
163 static void add_to_history(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
167 /* First try to update an existing history entry */
169 for (j
= s
->n_history
; j
> 0; j
--) {
171 if (s
->history_x
[i
] == x
) {
179 /* Drop old entries */
182 /* Calculate position for new entry */
183 j
= s
->history_idx
+ s
->n_history
;
193 /* And make sure we don't store more entries than fit in */
194 if (s
->n_history
> HISTORY_MAX
) {
195 s
->history_idx
+= s
->n_history
- HISTORY_MAX
;
196 REDUCE(s
->history_idx
);
197 s
->n_history
= HISTORY_MAX
;
201 static double avg_gradient(pa_smoother
*s
, pa_usec_t x
) {
202 unsigned i
, j
, c
= 0;
203 int64_t ax
= 0, ay
= 0, k
, t
;
206 /* Too few measurements, assume gradient of 1 */
207 if (s
->n_history
< s
->min_history
)
210 /* First, calculate average of all measurements */
212 for (j
= s
->n_history
; j
> 0; j
--) {
214 ax
+= s
->history_x
[i
];
215 ay
+= s
->history_y
[i
];
221 pa_assert(c
>= s
->min_history
);
225 /* Now, do linear regression */
229 for (j
= s
->n_history
; j
> 0; j
--) {
232 dx
= (int64_t) s
->history_x
[i
] - ax
;
233 dy
= (int64_t) s
->history_y
[i
] - ay
;
243 return (s
->monotonic
&& r
< 0) ? 0 : r
;
246 static void calc_abc(pa_smoother
*s
) {
247 pa_usec_t ex
, ey
, px
, py
;
256 /* We have two points: (ex|ey) and (px|py) with two gradients at
257 * these points de and dp. We do a polynomial
258 * interpolation of degree 3 with these 6 values */
260 ex
= s
->ex
; ey
= s
->ey
;
261 px
= s
->px
; py
= s
->py
;
262 de
= s
->de
; dp
= s
->dp
;
266 /* To increase the dynamic range and symplify calculation, we
267 * move these values to the origin */
268 kx
= (int64_t) px
- (int64_t) ex
;
269 ky
= (int64_t) py
- (int64_t) ey
;
271 /* Calculate a, b, c for y=ax^3+bx^2+cx */
273 s
->b
= (((double) (3*ky
)/kx
- dp
- 2*de
)) / kx
;
274 s
->a
= (dp
/kx
- 2*s
->b
- de
/kx
) / (3*kx
);
279 static void estimate(pa_smoother
*s
, pa_usec_t x
, pa_usec_t
*y
, double *deriv
) {
286 /* The requested point is right of the point where we wanted
287 * to be on track again, thus just linearly estimate */
289 t
= (int64_t) s
->py
+ (int64_t) (s
->dp
* (x
- s
->px
));
301 /* Ok, we're not yet on track, thus let's interpolate, and
302 * make sure that the first derivative is smooth */
310 *y
= (pa_usec_t
) ((double) x
* (s
->c
+ (double) x
* (s
->b
+ (double) x
* s
->a
)));
312 /* Move back from origin */
317 *deriv
= s
->c
+ ((double) x
* (s
->b
*2 + (double) x
* s
->a
*3));
320 /* Guarantee monotonicity */
323 if (deriv
&& *deriv
< 0)
328 void pa_smoother_put(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
339 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
344 /* First, we calculate the position we'd estimate for x, so that
345 * we can adjust our position smoothly from this one */
346 estimate(s
, x
, &ney
, &nde
);
347 s
->ex
= x
; s
->ey
= ney
; s
->de
= nde
;
352 /* Then, we add the new measurement to our history */
353 add_to_history(s
, x
, y
);
355 /* And determine the average gradient of the history */
356 s
->dp
= avg_gradient(s
, x
);
358 /* And calculate when we want to be on track again */
359 s
->px
= s
->ex
+ s
->adjust_time
;
360 s
->py
= s
->ry
+ s
->dp
*s
->adjust_time
;
362 s
->abc_valid
= FALSE
;
364 /* pa_log_debug("put(%llu | %llu) = %llu", (unsigned long long) (x + s->time_offset), (unsigned long long) x, (unsigned long long) y); */
367 pa_usec_t
pa_smoother_get(pa_smoother
*s
, pa_usec_t x
) {
376 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
378 estimate(s
, x
, &y
, NULL
);
382 /* Make sure the querier doesn't jump forth and back. */
383 pa_assert(x
>= s
->last_x
);
392 /* pa_log_debug("get(%llu | %llu) = %llu", (unsigned long long) (x + s->time_offset), (unsigned long long) x, (unsigned long long) y); */
397 void pa_smoother_set_time_offset(pa_smoother
*s
, pa_usec_t offset
) {
400 s
->time_offset
= offset
;
402 /* pa_log_debug("offset(%llu)", (unsigned long long) offset); */
405 void pa_smoother_pause(pa_smoother
*s
, pa_usec_t x
) {
411 /* pa_log_debug("pause(%llu)", (unsigned long long) x); */
417 void pa_smoother_resume(pa_smoother
*s
, pa_usec_t x
) {
423 pa_assert(x
>= s
->pause_time
);
425 /* pa_log_debug("resume(%llu)", (unsigned long long) x); */
428 s
->time_offset
+= x
- s
->pause_time
;
431 pa_usec_t
pa_smoother_translate(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y_delay
) {
441 x
= PA_LIKELY(x
>= s
->time_offset
) ? x
- s
->time_offset
: 0;
443 estimate(s
, x
, &ney
, &nde
);
445 /* Play safe and take the larger gradient, so that we wakeup
446 * earlier when this is used for sleeping */
450 /* pa_log_debug("translate(%llu) = %llu (%0.2f)", (unsigned long long) y_delay, (unsigned long long) ((double) y_delay / nde), nde); */
452 return (pa_usec_t
) ((double) y_delay
/ nde
);