]>
code.delx.au - pulseaudio/blob - src/pulsecore/time-smoother.c
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 50
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
;
66 pa_usec_t time_offset
;
68 pa_usec_t px
, py
; /* Point p, where we want to reach stability */
69 double dp
; /* Gradient we want at point p */
71 pa_usec_t ex
, ey
; /* Point e, which we estimated before and need to smooth to */
72 double de
; /* Gradient we estimated for point e */
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 */
81 /* Cached parameters for our interpolation polynomial y=ax^3+b^2+cx */
89 pa_smoother
* pa_smoother_new(pa_usec_t adjust_time
, pa_usec_t history_time
, pa_bool_t monotonic
) {
92 pa_assert(adjust_time
> 0);
93 pa_assert(history_time
> 0);
95 s
= pa_xnew(pa_smoother
, 1);
96 s
->adjust_time
= adjust_time
;
97 s
->history_time
= history_time
;
99 s
->monotonic
= monotonic
;
112 s
->abc_valid
= FALSE
;
119 void pa_smoother_free(pa_smoother
* s
) {
125 static void drop_old(pa_smoother
*s
, pa_usec_t x
) {
128 /* First drop items from history which are too old, but make sure
129 * to always keep two entries in the history */
131 for (j
= s
->n_history
; j
> 2; j
--) {
133 if (s
->history_x
[s
->history_idx
] + s
->history_time
>= x
) {
134 /* This item is still valid, and thus all following ones
135 * are too, so let's quit this loop */
139 /* Item is too old, let's drop it */
141 while (s
->history_idx
>= HISTORY_MAX
)
142 s
->history_idx
-= HISTORY_MAX
;
148 static void add_to_history(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
154 /* Calculate position for new entry */
155 j
= s
->history_idx
+ s
->n_history
;
156 while (j
>= HISTORY_MAX
)
166 /* And make sure we don't store more entries than fit in */
167 if (s
->n_history
>= HISTORY_MAX
) {
168 s
->history_idx
+= s
->n_history
- HISTORY_MAX
;
169 s
->n_history
= HISTORY_MAX
;
173 static double avg_gradient(pa_smoother
*s
, pa_usec_t x
) {
174 unsigned i
, j
, c
= 0;
175 int64_t ax
= 0, ay
= 0, k
, t
;
180 /* First, calculate average of all measurements */
182 for (j
= s
->n_history
; j
> 0; j
--) {
184 ax
+= s
->history_x
[i
];
185 ay
+= s
->history_y
[i
];
189 while (i
>= HISTORY_MAX
)
193 /* Too few measurements, assume gradient of 1 */
200 /* Now, do linear regression */
204 for (j
= s
->n_history
; j
> 0; j
--) {
207 dx
= (int64_t) s
->history_x
[i
] - ax
;
208 dy
= (int64_t) s
->history_y
[i
] - ay
;
214 while (i
>= HISTORY_MAX
)
220 return s
->monotonic
&& r
< 0 ? 0 : r
;
223 static void estimate(pa_smoother
*s
, pa_usec_t x
, pa_usec_t
*y
, double *deriv
) {
230 /* The requested point is right of the point where we wanted
231 * to be on track again, thus just linearly estimate */
233 t
= (int64_t) s
->py
+ (int64_t) (s
->dp
* (x
- s
->px
));
246 pa_usec_t ex
, ey
, px
, py
;
250 /* Ok, we're not yet on track, thus let's interpolate, and
251 * make sure that the first derivative is smooth */
253 /* We have two points: (ex|ey) and (px|py) with two gradients
254 * at these points de and dp. We do a polynomial interpolation
255 * of degree 3 with these 6 values */
257 ex
= s
->ex
; ey
= s
->ey
;
258 px
= s
->px
; py
= s
->py
;
259 de
= s
->de
; dp
= s
->dp
;
263 /* To increase the dynamic range and symplify calculation, we
264 * move these values to the origin */
265 kx
= (int64_t) px
- (int64_t) ex
;
266 ky
= (int64_t) py
- (int64_t) ey
;
268 /* Calculate a, b, c for y=ax^3+b^2+cx */
270 s
->b
= (((double) (3*ky
)/kx
- dp
- 2*de
)) / kx
;
271 s
->a
= (dp
/kx
- 2*s
->b
- de
/kx
) / (3*kx
);
280 *y
= (pa_usec_t
) ((double) x
* (s
->c
+ (double) x
* (s
->b
+ (double) x
* s
->a
)));
282 /* Move back from origin */
287 *deriv
= s
->c
+ ((double) x
* (s
->b
*2 + (double) x
* s
->a
*3));
290 /* Guarantee monotonicity */
298 if (deriv
&& *deriv
< 0)
303 void pa_smoother_put(pa_smoother
*s
, pa_usec_t x
, pa_usec_t y
) {
308 pa_assert(x
>= s
->time_offset
);
316 pa_assert(x
>= s
->ex
);
318 /* First, we calculate the position we'd estimate for x, so that
319 * we can adjust our position smoothly from this one */
320 estimate(s
, x
, &ney
, &nde
);
321 s
->ex
= x
; s
->ey
= ney
; s
->de
= nde
;
323 /* Then, we add the new measurement to our history */
324 add_to_history(s
, x
, y
);
326 /* And determine the average gradient of the history */
327 s
->dp
= avg_gradient(s
, x
);
329 /* And calculate when we want to be on track again */
330 s
->px
= x
+ s
->adjust_time
;
331 s
->py
= y
+ s
->dp
*s
->adjust_time
;
333 s
->abc_valid
= FALSE
;
336 pa_usec_t
pa_smoother_get(pa_smoother
*s
, pa_usec_t x
) {
340 pa_assert(x
>= s
->time_offset
);
348 pa_assert(x
>= s
->ex
);
350 estimate(s
, x
, &y
, NULL
);
354 void pa_smoother_set_time_offset(pa_smoother
*s
, pa_usec_t offset
) {
357 s
->time_offset
= offset
;
360 void pa_smoother_pause(pa_smoother
*s
, pa_usec_t x
) {
370 void pa_smoother_resume(pa_smoother
*s
, pa_usec_t x
) {
377 s
->time_offset
+= x
- s
->pause_time
;