1 /* _______ ____ __ ___ ___
2 * \ _ \ \ / \ / \ \ / / ' ' '
3 * | | \ \ | | || | \/ | . .
4 * | | | | | | || ||\ /| |
5 * | | | | | | || || \/ | | ' ' '
6 * | | | | | | || || | | . .
7 * | |_/ / \ \__// || | |
8 * /_______/ynamic \____/niversal /__\ /____\usic /| . . ibliotheque
9 * / \
10 * / . \
11 * resample.c - Resampling helper. / / \ \
12 * | < / \_
13 * By Bob and entheh. | \/ /\ /
14 * \_ / > /
15 * In order to find a good trade-off between | \ / /
16 * speed and accuracy in this code, some tests | ' /
17 * were carried out regarding the behaviour of \__/
18 * long long ints with gcc. The following code
19 * was tested:
20 *
21 * int a, b, c;
22 * c = ((long long)a * b) >> 16;
23 *
24 * DJGPP GCC Version 3.0.3 generated the following assembly language code for
25 * the multiplication and scaling, leaving the 32-bit result in EAX.
26 *
27 * movl -8(%ebp), %eax ; read one int into EAX
28 * imull -4(%ebp) ; multiply by the other; result goes in EDX:EAX
29 * shrdl $16, %edx, %eax ; shift EAX right 16, shifting bits in from EDX
30 *
31 * Note that a 32*32->64 multiplication is performed, allowing for high
32 * accuracy. On the Pentium 2 and above, shrdl takes two cycles (generally),
33 * so it is a minor concern when four multiplications are being performed
34 * (the cubic resampler). On the Pentium MMX and earlier, it takes four or
35 * more cycles, so this method is unsuitable for use in the low-quality
36 * resamplers.
37 *
38 * Since "long long" is a gcc-specific extension, we use LONG_LONG instead,
39 * defined in dumb.h. We may investigate later what code MSVC generates, but
40 * if it seems too slow then we suggest you use a good compiler.
41 *
42 * FIXME: these comments are somewhat out of date now.
43 */
44
45 #include <math.h>
46 #include "dumb.h"
47
48
49
50 /* Compile with -DHEAVYDEBUG if you want to make sure the pick-up function is
51 * called when it should be. There will be a considerable performance hit,
52 * since at least one condition has to be tested for every sample generated.
53 */
54 #ifdef HEAVYDEBUG
55 #define HEAVYASSERT(cond) ASSERT(cond)
56 #else
57 #define HEAVYASSERT(cond)
58 #endif
59
60
61
62 //#define MULSC(a, b) ((int)((LONG_LONG)(a) * (b) >> 16))
63 //#define MULSC(a, b) ((a) * ((b) >> 2) >> 14)
64 #define MULSC(a, b) ((int)((LONG_LONG)((a) << 4) * ((b) << 12) >> 32))
65
66
67
68 /* A global variable for controlling resampling quality wherever a local
69 * specification doesn't override it. The following values are valid:
70 *
71 * 0 - DUMB_RQ_ALIASING - fastest
72 * 1 - DUMB_RQ_LINEAR
73 * 2 - DUMB_RQ_CUBIC - nicest
74 *
75 * Values outside the range 0-2 will behave the same as the nearest
76 * value within the range.
77 */
78 int dumb_resampling_quality = 2;
79
80
81
dumb_reset_resampler(DUMB_RESAMPLER * resampler,sample_t * src,long pos,long start,long end)82 void dumb_reset_resampler(DUMB_RESAMPLER *resampler, sample_t *src, long pos, long start, long end)
83 {
84 resampler->src = src;
85 resampler->pos = pos;
86 resampler->subpos = 0;
87 resampler->start = start;
88 resampler->end = end;
89 resampler->dir = 1;
90 resampler->pickup = NULL;
91 resampler->pickup_data = NULL;
92 resampler->min_quality = 0;
93 resampler->max_quality = DUMB_RQ_N_LEVELS - 1;
94 resampler->x[2] = resampler->x[1] = resampler->x[0] = 0;
95 resampler->overshot = -1;
96 }
97
98
99
dumb_start_resampler(sample_t * src,long pos,long start,long end)100 DUMB_RESAMPLER *dumb_start_resampler(sample_t *src, long pos, long start, long end)
101 {
102 DUMB_RESAMPLER *resampler = malloc(sizeof(*resampler));
103 if (!resampler) return NULL;
104 dumb_reset_resampler(resampler, src, pos, start, end);
105 return resampler;
106 }
107
108
109
110 /* For convenience, returns nonzero on stop. */
process_pickup(DUMB_RESAMPLER * resampler)111 static int process_pickup(DUMB_RESAMPLER *resampler)
112 {
113 if (resampler->overshot < 0) {
114 resampler->overshot = 0;
115 dumb_resample(resampler, NULL, 2, 0, 1.0f);
116 resampler->x[0] = resampler->x[1];
117 }
118
119 for (;;) {
120 if (resampler->dir < 0) {
121 if (resampler->overshot >= 3 && resampler->pos+3 >= resampler->start) resampler->x[0] = resampler->src[resampler->pos+3];
122 if (resampler->overshot >= 2 && resampler->pos+2 >= resampler->start) resampler->x[1] = resampler->src[resampler->pos+2];
123 if (resampler->overshot >= 1 && resampler->pos+1 >= resampler->start) resampler->x[2] = resampler->src[resampler->pos+1];
124 resampler->overshot = resampler->start - resampler->pos - 1;
125 } else {
126 if (resampler->overshot >= 3 && resampler->pos-3 < resampler->end) resampler->x[0] = resampler->src[resampler->pos-3];
127 if (resampler->overshot >= 2 && resampler->pos-2 < resampler->end) resampler->x[1] = resampler->src[resampler->pos-2];
128 if (resampler->overshot >= 1 && resampler->pos-1 < resampler->end) resampler->x[2] = resampler->src[resampler->pos-1];
129 resampler->overshot = resampler->pos - resampler->end;
130 }
131
132 if (resampler->overshot < 0) {
133 resampler->overshot = 0;
134 return 0;
135 }
136
137 if (!resampler->pickup) {
138 resampler->dir = 0;
139 return 1;
140 }
141 (*resampler->pickup)(resampler, resampler->pickup_data);
142 if (resampler->dir == 0) return 1;
143 ASSERT(resampler->dir == -1 || resampler->dir == 1);
144 }
145 }
146
147
148
149 /* Executes the content 'iterator' times.
150 * Clobbers the 'iterator' variable.
151 * The loop is unrolled by four.
152 */
153 #define LOOP4(iterator, CONTENT) \
154 { \
155 if ((iterator) & 2) { \
156 CONTENT; \
157 CONTENT; \
158 } \
159 if ((iterator) & 1) { \
160 CONTENT; \
161 } \
162 (iterator) >>= 2; \
163 while (iterator) { \
164 CONTENT; \
165 CONTENT; \
166 CONTENT; \
167 CONTENT; \
168 (iterator)--; \
169 } \
170 }
171
172
173
dumb_resample(DUMB_RESAMPLER * resampler,sample_t * dst,long dst_size,float volume,float delta)174 long dumb_resample(DUMB_RESAMPLER *resampler, sample_t *dst, long dst_size, float volume, float delta)
175 {
176 int dt;
177 int vol;
178 long done;
179 long todo;
180 int quality;
181
182 if (!resampler || resampler->dir == 0) return 0;
183 ASSERT(resampler->dir == -1 || resampler->dir == 1);
184
185 done = 0;
186 dt = (int)(delta * 65536.0 + 0.5);
187 vol = (int)floor(volume * 65536.0 + 0.5);
188
189 if (vol == 0) dst = NULL;
190
191 quality = dumb_resampling_quality;
192 if (quality > resampler->max_quality) quality = resampler->max_quality;
193 else if (quality < resampler->min_quality) quality = resampler->min_quality;
194
195 while (done < dst_size) {
196 if (process_pickup(resampler)) return done;
197
198 if ((resampler->dir ^ dt) < 0)
199 dt = -dt;
200
201 if (resampler->dir < 0)
202 todo = (long)((((LONG_LONG)(resampler->pos - resampler->start) << 16) + resampler->subpos - dt) / -dt);
203 else
204 todo = (long)((((LONG_LONG)(resampler->end - resampler->pos) << 16) - resampler->subpos - 1 + dt) / dt);
205
206 if (todo < 0)
207 todo = 0;
208 else if (todo > dst_size - done)
209 todo = dst_size - done;
210
211 done += todo;
212
213 {
214 sample_t *src = resampler->src;
215 long pos = resampler->pos;
216 int subpos = resampler->subpos;
217 long diff = pos;
218 long overshot;
219 if (resampler->dir < 0) {
220 if (!dst) {
221 /* Silence or simulation */
222 LONG_LONG new_subpos = subpos + dt * todo;
223 pos += (long)(new_subpos >> 16);
224 subpos = (long)new_subpos & 65535;
225 } else if (quality <= DUMB_RQ_ALIASING) {
226 /* Aliasing, backwards */
227 sample_t xbuf[2];
228 sample_t *x = &xbuf[0];
229 sample_t *xstart;
230 xbuf[0] = resampler->x[1];
231 xbuf[1] = resampler->x[2];
232 while (todo && x < &xbuf[2]) {
233 HEAVYASSERT(pos >= resampler->start);
234 *dst++ += MULSC(x[0], vol);
235 subpos += dt;
236 pos += subpos >> 16;
237 x -= subpos >> 16;
238 subpos &= 65535;
239 todo--;
240 }
241 x = xstart = &src[pos];
242 LOOP4(todo,
243 *dst++ += MULSC(x[2], vol);
244 subpos += dt;
245 x += subpos >> 16;
246 subpos &= 65535;
247 );
248 pos += x - xstart;
249 } else if (quality <= DUMB_RQ_LINEAR) {
250 /* Linear interpolation, backwards */
251 sample_t xbuf[3];
252 sample_t *x = &xbuf[1];
253 xbuf[0] = resampler->x[1];
254 xbuf[1] = resampler->x[2];
255 xbuf[2] = src[pos];
256 while (todo && x < &xbuf[3]) {
257 HEAVYASSERT(pos >= resampler->start);
258 *dst++ += MULSC(x[0] + MULSC(x[-1] - x[0], subpos), vol);
259 subpos += dt;
260 pos += subpos >> 16;
261 x -= subpos >> 16;
262 subpos &= 65535;
263 todo--;
264 }
265 x = &src[pos];
266 LOOP4(todo,
267 HEAVYASSERT(pos >= resampler->start);
268 *dst++ += MULSC(x[1] + MULSC(x[2] - x[1], subpos), vol);
269 subpos += dt;
270 pos += subpos >> 16;
271 x += subpos >> 16;
272 subpos &= 65535;
273 );
274 } else {
275 /* Cubic interpolation, backwards */
276 sample_t xbuf[6];
277 sample_t *x = &xbuf[3];
278 sample_t *lastx = NULL;
279 int a = 0, b = 0, c = 0;
280 xbuf[0] = resampler->x[0];
281 xbuf[1] = resampler->x[1];
282 xbuf[2] = resampler->x[2];
283 xbuf[3] = src[pos];
284 if (pos-1 >= resampler->start) xbuf[4] = src[pos-1];
285 if (pos-2 >= resampler->start) xbuf[5] = src[pos-2];
286 while (todo && x < &xbuf[6]) {
287 HEAVYASSERT(pos >= resampler->start);
288 if (lastx != x) {
289 lastx = x;
290 a = (((x[-1] - x[-2]) << 1) + (x[-1] - x[-2]) + (x[-3] - x[0])) >> 1;
291 b = (x[-2] << 1) + x[0] - ((5 * x[-1] + x[-3]) >> 1);
292 c = (x[-2] - x[0]) >> 1;
293 }
294 *dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-1], vol);
295 subpos += dt;
296 pos += subpos >> 16;
297 x -= subpos >> 16;
298 subpos &= 65535;
299 todo--;
300 }
301 x = &src[pos];
302 lastx = NULL;
303 LOOP4(todo,
304 HEAVYASSERT(pos >= resampler->start);
305 if (lastx != x) {
306 lastx = x;
307 a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
308 b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
309 c = (x[2] - x[0]) >> 1;
310 }
311 *dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[1], vol);
312 subpos += dt;
313 pos += subpos >> 16;
314 x += subpos >> 16;
315 subpos &= 65535;
316 );
317 }
318 diff = diff - pos;
319 overshot = resampler->start - pos - 1;
320 if (diff >= 3) {
321 resampler->x[0] = overshot >= 3 ? 0 : src[pos+3];
322 resampler->x[1] = overshot >= 2 ? 0 : src[pos+2];
323 resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
324 } else if (diff >= 2) {
325 resampler->x[0] = resampler->x[2];
326 resampler->x[1] = overshot >= 2 ? 0 : src[pos+2];
327 resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
328 } else if (diff >= 1) {
329 resampler->x[0] = resampler->x[1];
330 resampler->x[1] = resampler->x[2];
331 resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
332 }
333 } else {
334 if (!dst) {
335 /* Silence or simulation */
336 LONG_LONG new_subpos = subpos + dt * todo;
337 pos += (long)(new_subpos >> 16);
338 subpos = (long)new_subpos & 65535;
339 } else if (dumb_resampling_quality <= DUMB_RQ_ALIASING) {
340 /* Aliasing, forwards */
341 sample_t xbuf[2];
342 sample_t *x = &xbuf[0];
343 sample_t *xstart;
344 xbuf[0] = resampler->x[1];
345 xbuf[1] = resampler->x[2];
346 while (todo && x < &xbuf[2]) {
347 HEAVYASSERT(pos < resampler->end);
348 *dst++ += MULSC(x[0], vol);
349 subpos += dt;
350 pos += subpos >> 16;
351 x += subpos >> 16;
352 subpos &= 65535;
353 todo--;
354 }
355 x = xstart = &src[pos];
356 LOOP4(todo,
357 *dst++ += MULSC(x[-2], vol);
358 subpos += dt;
359 x += subpos >> 16;
360 subpos &= 65535;
361 );
362 pos += x - xstart;
363 } else if (dumb_resampling_quality <= DUMB_RQ_LINEAR) {
364 /* Linear interpolation, forwards */
365 sample_t xbuf[3];
366 sample_t *x = &xbuf[1];
367 xbuf[0] = resampler->x[1];
368 xbuf[1] = resampler->x[2];
369 xbuf[2] = src[pos];
370 while (todo && x < &xbuf[3]) {
371 HEAVYASSERT(pos < resampler->end);
372 *dst++ += MULSC(x[-1] + MULSC(x[0] - x[-1], subpos), vol);
373 subpos += dt;
374 pos += subpos >> 16;
375 x += subpos >> 16;
376 subpos &= 65535;
377 todo--;
378 }
379 x = &src[pos];
380 LOOP4(todo,
381 HEAVYASSERT(pos < resampler->end);
382 *dst++ += MULSC(x[-2] + MULSC(x[-1] - x[-2], subpos), vol);
383 subpos += dt;
384 pos += subpos >> 16;
385 x += subpos >> 16;
386 subpos &= 65535;
387 );
388 } else {
389 /* Cubic interpolation, forwards */
390 sample_t xbuf[6];
391 sample_t *x = &xbuf[3];
392 sample_t *lastx = NULL;
393 int a = 0, b = 0, c = 0;
394 xbuf[0] = resampler->x[0];
395 xbuf[1] = resampler->x[1];
396 xbuf[2] = resampler->x[2];
397 xbuf[3] = src[pos];
398 if (pos+1 < resampler->end) xbuf[4] = src[pos+1];
399 if (pos+2 < resampler->end) xbuf[5] = src[pos+2];
400 while (todo && x < &xbuf[6]) {
401 HEAVYASSERT(pos < resampler->end);
402 if (lastx != x) {
403 lastx = x;
404 a = (((x[-2] - x[-1]) << 1) + (x[-2] - x[-1]) + (x[0] - x[-3])) >> 1;
405 b = (x[-1] << 1) + x[-3] - ((5 * x[-2] + x[0]) >> 1);
406 c = (x[-1] - x[-3]) >> 1;
407 }
408 *dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-2], vol);
409 subpos += dt;
410 pos += subpos >> 16;
411 x += subpos >> 16;
412 subpos &= 65535;
413 todo--;
414 }
415 x = &src[pos];
416 lastx = NULL;
417 LOOP4(todo,
418 HEAVYASSERT(pos < resampler->end);
419 if (lastx != x) {
420 lastx = x;
421 a = (((x[-2] - x[-1]) << 1) + (x[-2] - x[-1]) + (x[0] - x[-3])) >> 1;
422 b = (x[-1] << 1) + x[-3] - ((5 * x[-2] + x[0]) >> 1);
423 c = (x[-1] - x[-3]) >> 1;
424 }
425 *dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-2], vol);
426 subpos += dt;
427 pos += subpos >> 16;
428 x += subpos >> 16;
429 subpos &= 65535;
430 );
431 }
432 diff = pos - diff;
433 overshot = pos - resampler->end;
434 if (diff >= 3) {
435 resampler->x[0] = overshot >= 3 ? 0 : src[pos-3];
436 resampler->x[1] = overshot >= 2 ? 0 : src[pos-2];
437 resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
438 } else if (diff >= 2) {
439 resampler->x[0] = resampler->x[2];
440 resampler->x[1] = overshot >= 2 ? 0 : src[pos-2];
441 resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
442 } else if (diff >= 1) {
443 resampler->x[0] = resampler->x[1];
444 resampler->x[1] = resampler->x[2];
445 resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
446 }
447 }
448 resampler->pos = pos;
449 resampler->subpos = subpos;
450 }
451 }
452
453 return done;
454 }
455
456
457
dumb_resample_get_current_sample(DUMB_RESAMPLER * resampler,float volume)458 sample_t dumb_resample_get_current_sample(DUMB_RESAMPLER *resampler, float volume)
459 {
460 int vol;
461 sample_t *src;
462 long pos;
463 int subpos;
464 int quality;
465
466 if (!resampler || resampler->dir == 0) return 0;
467 ASSERT(resampler->dir == -1 || resampler->dir == 1);
468
469 if (process_pickup(resampler)) return 0;
470
471 vol = (int)floor(volume * 65536.0 + 0.5);
472 if (vol == 0) return 0;
473
474 quality = dumb_resampling_quality;
475 if (quality > resampler->max_quality) quality = resampler->max_quality;
476 else if (quality < resampler->min_quality) quality = resampler->min_quality;
477
478 src = resampler->src;
479 pos = resampler->pos;
480 subpos = resampler->subpos;
481
482 if (resampler->dir < 0) {
483 HEAVYASSERT(pos >= resampler->start);
484 if (dumb_resampling_quality <= 0) {
485 /* Aliasing, backwards */
486 return MULSC(src[pos], vol);
487 } else if (quality <= DUMB_RQ_LINEAR) {
488 /* Linear interpolation, backwards */
489 return MULSC(resampler->x[2] + MULSC(resampler->x[1] - resampler->x[2], subpos), vol);
490 } else {
491 /* Cubic interpolation, backwards */
492 sample_t *x = resampler->x;
493 int a, b, c;
494 a = (((x[2] - x[1]) << 1) + (x[2] - x[1]) + (x[0] - src[pos])) >> 1;
495 b = (x[1] << 1) + src[pos] - ((5 * x[2] + x[0]) >> 1);
496 c = (x[1] - src[pos]) >> 1;
497 return MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[2], vol);
498 }
499 } else {
500 HEAVYASSERT(pos < resampler->end);
501 if (dumb_resampling_quality <= 0) {
502 /* Aliasing */
503 return MULSC(src[pos], vol);
504 } else if (dumb_resampling_quality <= DUMB_RQ_LINEAR) {
505 /* Linear interpolation, forwards */
506 return MULSC(resampler->x[1] + MULSC(resampler->x[2] - resampler->x[1], subpos), vol);
507 } else {
508 /* Cubic interpolation, forwards */
509 sample_t *x = resampler->x;
510 int a, b, c;
511 a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (src[pos] - x[0])) >> 1;
512 b = (x[2] << 1) + x[0] - ((5 * x[1] + src[pos]) >> 1);
513 c = (x[2] - x[0]) >> 1;
514 return MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[1], vol);
515 }
516 }
517 }
518
519
520
dumb_end_resampler(DUMB_RESAMPLER * resampler)521 void dumb_end_resampler(DUMB_RESAMPLER *resampler)
522 {
523 if (resampler)
524 free(resampler);
525 }
526
527
528
529 #if 0
530 /* The following macro is used to overcome the fact that most C
531 * compilers (including gcc and MSVC) can't correctly multiply signed
532 * integers outside the range -32768 to 32767. i86 assembler versions
533 * don't need to use this method, since the processor does in fact
534 * have instructions to multiply large numbers correctly - which
535 * means using assembly language could make a significant difference
536 * to the speed.
537 *
538 * The basic method is as follows. We halve the subposition (how far
539 * we are between samples), so it never exceeds 32767. We also halve
540 * the delta, which is the amount to be added to the subposition each
541 * time. Then we unroll the loop twofold, so that we can add the lost
542 * one every other time if necessary (since the halving may have
543 * resulted in rounding down).
544 *
545 * This method doesn't incur any cumulative inaccuracies. There is a
546 * very slight loss of quality, which I challenge anyone to notice -
547 * but the position will advance at *exactly* the same rate as it
548 * would if we didn't use this method. This also means the pitch is
549 * exactly the same, which may even make a difference to trained
550 * musicians when resampling down a lot :)
551 *
552 * Each time this macro is invoked, DO_RESAMPLE(inc) must be defined
553 * to calculate the samples by the appropriate equation (linear,
554 * cubic, etc.). See the individual cases for examples of how this is
555 * done.
556 */
557 #define MAKE_RESAMPLER() \
558 { \
559 if (dt & 1) { \
560 long todo2; \
561 \
562 dt >>= 1; \
563 \
564 if (src_subpos & 1) { \
565 src_subpos >>= 1; \
566 DO_RESAMPLE(1); \
567 todo--; \
568 } else \
569 src_subpos >>= 1; \
570 \
571 todo2 = todo >> 1; \
572 \
573 while (todo2) { \
574 DO_RESAMPLE(0); \
575 DO_RESAMPLE(1); \
576 todo2--; \
577 } \
578 \
579 if (todo & 1) { \
580 DO_RESAMPLE(0); \
581 src_subpos = (src_subpos << 1) | 1; \
582 } else \
583 src_subpos <<= 1; \
584 \
585 todo = 0; \
586 dt = (dt << 1) | 1; \
587 } else { \
588 long subposbit = src_subpos & 1; \
589 dt >>= 1; \
590 src_subpos >>= 1; \
591 \
592 if (todo & 1) { \
593 DO_RESAMPLE(0); \
594 } \
595 \
596 todo >>= 1; \
597 \
598 while (todo) { \
599 DO_RESAMPLE(0); \
600 DO_RESAMPLE(0); \
601 todo--; \
602 } \
603 \
604 src_subpos = (src_subpos << 1) | subposbit; \
605 dt <<= 1; \
606 } \
607 }
608
609
610
611 sample_t dumb_resample_get_current_sample(
612 sample_t *src, long *_src_pos, int *_src_subpos,
613 long src_start, long src_end,
614 float volume, int *_dir,
615 DUMB_RESAMPLE_PICKUP pickup, void *pickup_data
616 )
617 {
618 long src_pos = *_src_pos;
619 int src_subpos = *_src_subpos;
620 int dir = _dir ? *_dir : 1;
621
622 sample_t value = 0;
623
624 if (dir == 0)
625 return 0;
626
627 ASSERT(dir == 1 || dir == -1);
628
629 if (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end)) {
630
631 /* If there's no pick-up function, we stop. */
632 if (!pickup) {
633 dir = 0;
634 goto end;
635 }
636
637 /* Process the pick-up. It may need invoking more than once. */
638 do {
639 dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);
640
641 if (dir == 0)
642 goto end;
643
644 ASSERT(dir == 1 || dir == -1);
645 } while (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end));
646 }
647
648 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end));
649
650 if (dumb_resampling_quality == 0) {
651 /* Aliasing (coarse) */
652 int volume_fact = (int)(volume * 16384.0);
653 value = (src[src_pos] * volume_fact) >> 14;
654 } else if (dumb_resampling_quality <= 2) {
655 /* Linear interpolation */
656 int volume_fact = (int)(volume * 16384.0);
657 int subpos = src_subpos >> 1;
658 value = ((src[src_pos] + ((((src[src_pos + 1] - src[src_pos]) >> 1) * subpos) >> 14)) * volume_fact) >> 14;
659 } else if (dumb_resampling_quality == 3) {
660 /* Quadratic interpolation */
661 int volume_fact = (int)(volume * 16384.0);
662 int a, b;
663 sample_t *x;
664 int subpos = src_subpos >> 1;
665 x = &src[src_pos];
666 a = ((x[0] + x[2]) >> 1) - x[1];
667 b = ((x[2] - x[0]) >> 1) - (a << 1);
668 value = (((((((a * subpos) >> 15) + b) * subpos) >> 15) + x[0]) * volume_fact) >> 14;
669 } else {
670 /* Cubic interpolation */
671 int volume_fact = (int)(volume * 16384.0);
672 int a, b, c;
673 sample_t *x;
674 int subpos = src_subpos >> 1;
675 x = &src[src_pos];
676 a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
677 b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
678 c = (x[2] - x[0]) >> 1;
679 value = (((int)(((LONG_LONG)((int)(((LONG_LONG)((int)(((LONG_LONG)a * subpos) >> 15) + b) * subpos) >> 15) + c) * subpos) >> 15) + x[1]) * volume_fact) >> 14;
680 }
681
682 end:
683
684 *_src_pos = src_pos;
685 *_src_subpos = src_subpos;
686 if (_dir) *_dir = dir;
687
688 return value;
689 }
690
691
692
693 long dumb_resample(
694 sample_t *src, long *_src_pos, int *_src_subpos,
695 long src_start, long src_end,
696 sample_t *dst, long dst_size,
697 float volume, float delta, int *_dir,
698 DUMB_RESAMPLE_PICKUP pickup, void *pickup_data
699 )
700 {
701 int dt = (int)(delta * 65536.0 + 0.5);
702 long s = 0; /* Current position in the destination buffer */
703
704 long src_pos = *_src_pos;
705 int src_subpos = *_src_subpos;
706 int dir = _dir ? *_dir : 1;
707
708 int linear_average;
709
710 if (dir == 0)
711 return 0;
712
713 ASSERT(dir == 1 || dir == -1);
714
715 linear_average = dst && dumb_resampling_quality >= 2 && dt > 65536;
716
717 if (dir < 0) dt = -dt;
718
719 if (linear_average)
720 volume /= delta;
721
722 while (s < dst_size) {
723
724 long todo;
725
726 /* Process pick-ups first, just in case. */
727
728 if (linear_average) {
729
730 /* For linear average, the pick-up point could split a sum into
731 * two parts. We handle this by putting the pick-up code inside
732 * the summing loop. Note that this code is only executed when we
733 * know that a pick-up is necessary somewhere during this sum
734 * (although it is always executed once for the first sample).
735 * We use a separate loop further down when we know we won't have
736 * to do a pick-up, so the condition does not need testing inside
737 * the loop.
738 */
739
740 float sum;
741 long i;
742 int advance;
743 int x[3];
744
745 advance = src_subpos + dt;
746
747 /* Make these negative. Then they stay within the necessary
748 * range for integer multiplication, -32768 to 32767 ;)
749 */
750 x[0] = ~(src_subpos >> 1); /* = -1 - (src_subpos >> 1) */
751 x[2] = x[0] ^ 0x7FFF; /* = -32768 + (src_subpos >> 1) */
752
753 sum = (float)(-((src[src_pos] * (x+1)[dir]) >> 15));
754
755 i = src_pos + (advance >> 16);
756 src_pos += dir;
757 src_subpos = (dir >> 1) & 65535; /* changes 1,-1 to 0,65535 */
758
759 advance &= 65535;
760
761 /* i is the index of the first sample NOT to sum fully,
762 * regardless of the direction of resampling.
763 */
764
765 while (dir < 0 ? (i < src_start) : (i >= src_end)) {
766 if (dir < 0) {
767 while (src_pos >= src_start)
768 sum += src[src_pos--];
769 } else {
770 while (src_pos < src_end)
771 sum += src[src_pos++];
772 }
773
774 i -= src_pos;
775 /* i is now the number of samples left to sum fully, except
776 * it's negative if we're going backwards.
777 */
778
779 if (!pickup) {
780 dir = 0;
781 goto endsum;
782 }
783
784 dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);
785
786 if (dir == 0)
787 goto endsum;
788
789 ASSERT(dir == 1 || dir == -1);
790
791 if ((dir ^ dt) < 0) {
792 dt = -dt;
793 advance ^= 65535;
794 i = -i;
795 }
796
797 i += src_pos;
798 /* There, i is back to normal. */
799 }
800
801 for (; src_pos != i; src_pos += dir)
802 sum += src[src_pos];
803
804 src_subpos = advance;
805
806 x[2] = src_subpos >> 1;
807 x[0] = x[2] ^ 0x7FFF; /* = 32767 - (src_subpos >> 1) */
808
809 sum += (src[src_pos] * (x+1)[dir]) >> 15;
810
811 endsum:
812
813 sum *= volume;
814 dst[s] += (int)sum;
815
816 s++;
817
818 if (dir == 0)
819 break;
820
821 } else if (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end)) {
822
823 /* If there's no pick-up function, we stop. */
824 if (!pickup) {
825 dir = 0;
826 break;
827 }
828
829 /* Process the pick-up. It may need invoking more than once. */
830 do {
831 dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);
832
833 if (dir == 0)
834 goto end;
835
836 ASSERT(dir == 1 || dir == -1);
837 } while (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end));
838
839 /* Update sign of dt to match that of dir. */
840 if ((dir ^ dt) < 0)
841 dt = -dt;
842 }
843
844 /* Work out how many contiguous samples we can now render. */
845 if (dir < 0)
846 todo = (long)((((LONG_LONG)(src_pos - src_start) << 16) + src_subpos) / -dt);
847 else
848 todo = (long)((((LONG_LONG)(src_end - src_pos) << 16) - src_subpos - 1) / dt);
849
850 /* The above equations work out how many complete dt-sized
851 * intervals there are between the current position and the loop
852 * point (provided there is a little fractional extra). The linear
853 * average function needs complete intervals - but the other
854 * resamplers only read a sample from the beginning of each interval,
855 * so they can process one extra sample in their main loops (so we
856 * increment todo in a moment).
857 *
858 * The linear average function makes up the extra sample using the
859 * specialised pick-up code above.
860 *
861 * Note that our above pick-up process should have absolutely ensured
862 * that the result of this function will be nonnegative.
863 */
864
865 ASSERT(todo >= 0);
866
867 if (!linear_average)
868 todo++;
869
870 /* Of course we don't want to overrun the output buffer! */
871 if (todo > dst_size - s)
872 todo = dst_size - s;
873
874 if (!dst) {
875
876 LONG_LONG t = src_subpos + (LONG_LONG)dt * todo;
877 src_pos += (long)(t >> 16);
878 src_subpos = (int)t & 0xFFFFl;
879
880 s += todo;
881
882 } else if (linear_average) {
883
884 float sum;
885 long i;
886 int advance;
887 int x[3];
888
889 while (todo) {
890
891 advance = src_subpos + dt;
892
893 /* Make these negative. Then they stay within the necessary
894 * range for integer multiplication, -32768 to 32767 ;)
895 */
896 x[0] = ~(src_subpos >> 1); /* = -1 - (src_subpos >> 1) */
897 x[2] = x[0] ^ 0x7FFF; /* = -32768 + (src_subpos >> 1) */
898
899 sum = (float)(-((src[src_pos] * (x+1)[dir]) >> 15));
900
901 i = src_pos + (advance >> 16);
902 src_pos += dir;
903 src_subpos = (dir >> 1) & 65535; /* changes 1,-1 to 0,65535 */
904
905 advance &= 65535;
906
907 /* i is the index of the first sample NOT to sum fully,
908 * regardless of the direction of resampling.
909 */
910
911 HEAVYASSERT(dir < 0 ? (i >= src_start) : (i < src_end));
912
913 for (; src_pos != i; src_pos += dir)
914 sum += src[src_pos];
915
916 src_subpos = advance;
917
918 x[2] = src_subpos >> 1;
919 x[0] = x[2] ^ 0x7FFF; /* = 32767 - (src_subpos >> 1) */
920
921 sum += (src[src_pos] * (x+1)[dir]) >> 15;
922
923 sum *= volume;
924 dst[s] += (int)sum;
925
926 s++;
927 todo--;
928 }
929
930 } else if (dumb_resampling_quality == 0 || (dumb_resampling_quality == 1 && delta >= 1.0)) {
931
932 /* Aliasing (coarse) */
933 int volume_fact = (int)(volume * 16384.0);
934
935 do {
936 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end));
937 dst[s] += ((src[src_pos] * volume_fact) >> 14);
938 src_subpos += dt;
939 src_pos += src_subpos >> 16;
940 src_subpos &= 0xFFFFl;
941 s++;
942 } while (--todo);
943
944 } else if (dumb_resampling_quality <= 2) {
945
946 /* Linear interpolation */
947 int volume_fact = (int)(volume * 16384.0);
948
949 #define DO_RESAMPLE(inc) \
950 { \
951 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
952 \
953 dst[s] += (((src[src_pos] + ((((src[src_pos + 1] - src[src_pos]) >> 1) * src_subpos) >> 14)) * volume_fact) >> 14); \
954 \
955 src_subpos += dt + inc; \
956 src_pos += src_subpos >> 15; \
957 src_subpos &= 0x7FFFl; \
958 s++; \
959 }
960
961 MAKE_RESAMPLER();
962
963 #undef DO_RESAMPLE
964
965 } else if (dumb_resampling_quality == 3) {
966
967 /* Quadratic interpolation */
968
969 int volume_fact = (int)(volume * 16384.0);
970 int a = 0, b = 0;
971 sample_t *x = NULL;
972 int last_src_pos = -1;
973
974 /* AIM: no integer multiplicands must transcend the range -32768 to 32767.
975 * This limitation is imposed by most compilers, including gcc and MSVC.
976 *
977 * a = 0.5 * (s0 + s2) - s1
978 * b = -1.5 * s0 + 2 * s1 - 0.5 * s2
979 * c = s0
980 *
981 * s = (a * t + b) * t + c
982 *
983 * In fixed-point:
984 *
985 * a = ((s0 + s2) >> 1) - s1
986 * b = ((-3 * s0 - s2) >> 1) + (s1 << 1)
987 *
988 * s = (((((a * t) >> 16) + b) * t) >> 16) + s0
989 *
990 * With t halved (since t can reach 65535):
991 *
992 * s = (((((a * t) >> 15) + b) * t) >> 15) + s0
993 *
994 * a currently reaches 65536
995 * b currently reaches 131072
996 *
997 * So we must use aon2
998 *
999 * s = (((((aon2 * t) >> 14) + b) * t) >> 15) + s0
1000 *
1001 * ((aon2 * t) >> 14) + b is 5 times too big
1002 * so we must divide by 8
1003 *
1004 * s = (((((aon2 * t) >> 17) + bon8) * t) >> 12) + s0
1005 *
1006 * aon2 = ((s0 + s2) >> 2) - (s1 >> 1)
1007 * bon8 = ((-3 * s0 - s2) >> 4) + (s1 >> 2)
1008 * or:
1009 * bon8 = ((s2 - s0) >> 4) - (aon2 >> 1)
1010 */
1011
1012 /* Unh4x0r3d version:
1013 #define DO_RESAMPLE(inc) \
1014 { \
1015 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
1016 \
1017 if (src_pos != last_src_pos) { \
1018 last_src_pos = src_pos; \
1019 x = &src[src_pos]; \
1020 a = ((x[0] + x[2]) >> 2) - (x[1] >> 1); \
1021 b = ((x[2] - x[0]) >> 4) - (a >> 1); \
1022 } \
1023 \
1024 dst[s] += ((((((((a * src_subpos) >> 17) + b) * src_subpos) >> 12) + x[0]) * volume_fact) >> 14); \
1025 \
1026 src_subpos += dt + inc; \
1027 src_pos += src_subpos >> 15; \
1028 src_subpos &= 0x7FFFl; \
1029 s++; \
1030 }
1031 */
1032
1033 /* H4x0r3d version: */
1034 #define DO_RESAMPLE(inc) \
1035 { \
1036 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
1037 \
1038 if (src_pos != last_src_pos) { \
1039 last_src_pos = src_pos; \
1040 x = &src[src_pos]; \
1041 a = ((x[0] + x[2]) >> 1) - x[1]; \
1042 b = ((x[2] - x[0]) >> 1) - (a << 1); \
1043 } \
1044 \
1045 dst[s] += ((((((((a * src_subpos) >> 15) + b) * src_subpos) >> 15) + x[0]) * volume_fact) >> 14); \
1046 \
1047 src_subpos += dt + inc; \
1048 src_pos += src_subpos >> 15; \
1049 src_subpos &= 0x7FFFl; \
1050 s++; \
1051 }
1052
1053 MAKE_RESAMPLER();
1054
1055 #undef DO_RESAMPLE
1056
1057 } else {
1058
1059 /* Cubic interpolation */
1060
1061 int volume_fact = (int)(volume * 16384.0);
1062 int a = 0, b = 0, c = 0;
1063 sample_t *x = NULL;
1064 int last_src_pos = -1;
1065
1066 /* AIM: never multiply integers outside the range -32768 to 32767.
1067 *
1068 * a = 1.5f * (x[1] - x[2]) + (x[3] - x[0]) * 0.5f;
1069 * b = 2.0f * x[2] + x[0] - 2.5f * x[1] - x[3] * 0.5f;
1070 * c = (x[2] - x[0]) * 0.5f;
1071 *
1072 * s = ((a * t + b) * t + c) * t + x[1];
1073 *
1074 * Fixed-point version:
1075 *
1076 * a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
1077 * b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
1078 * c = (x[2] - x[0]) >> 1;
1079 *
1080 * s = ((((((((a * t) >> 15) + b) * t) >> 15) + c) * t) >> 15) + x[1];
1081 * (with t already halved, maximum 32767)
1082 *
1083 * a is in (((1+1)*2)+(1+1)+(1+1))/2 = 8 times the required range
1084 * b is in (1*2)+1+((5*1+1)/2) = 6 times
1085 * c is in the required range
1086 *
1087 * We must use aon8
1088 *
1089 * s = ((((((((aon8 * t) >> 12) + b) * t) >> 15) + c) * t) >> 15) + x[1];
1090 *
1091 * But ((aon8 * t) >> 12) is in 2^(15+15-12) = 2^18 = 8 times
1092 * b is in 6 times
1093 * so we divide both ((aon8 * t) >> 12) and b by 16
1094 *
1095 * s = ((((((((aon8 * t) >> 16) + bon16) * t) >> 11) + c) * t) >> 15) + x[1];
1096 *
1097 * ((... + bon16) * t) >> 11 is 16 times too big
1098 * c is in the correct range
1099 * we must divide both by 32
1100 *
1101 * s = ((((((((aon8 * t) >> 16) + bon16) * t) >> 16) + con32) * t) >> 10) + x[1];
1102 *
1103 * aon8 = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 4;
1104 * bon16 = ((x[2] << 2) + (x[0] << 1) - (5 * x[1] + x[3])) >> 5;
1105 * con32 = (x[2] - x[0]) >> 6;
1106 *
1107 * A lot of accuracy is lost here. It is quite likely that some
1108 * of the above would cancel anyway, so the scaling down wouldn't
1109 * have to be so severe. However, I'm not in the mood to work it
1110 * out now :P
1111 *
1112 * It may also be worth investigating whether doing this stuff
1113 * in floats would be faster.
1114 */
1115
1116 /* Unh4x0r3d version:
1117 #define DO_RESAMPLE(inc) \
1118 { \
1119 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
1120 \
1121 if (src_pos != last_src_pos) { \
1122 last_src_pos = src_pos; \
1123 x = &src[src_pos]; \
1124 a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 4; \
1125 b = ((x[2] << 2) + (x[0] << 1) - (5 * x[1] + x[3])) >> 5; \
1126 c = (x[2] - x[0]) >> 6; \
1127 } \
1128 \
1129 dst[s] += ((((((((((a * src_subpos) >> 16) + b) * src_subpos) >> 16) + c) * src_subpos) >> 10) + x[1]) * volume_fact) >> 14; \
1130 \
1131 src_subpos += dt + inc; \
1132 src_pos += src_subpos >> 15; \
1133 src_subpos &= 0x7FFFl; \
1134 s++; \
1135 }
1136 */
1137
1138 /* H4x0r3d version: */
1139 #define DO_RESAMPLE(inc) \
1140 { \
1141 HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
1142 \
1143 if (src_pos != last_src_pos) { \
1144 last_src_pos = src_pos; \
1145 x = &src[src_pos]; \
1146 a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1; \
1147 b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1); \
1148 c = (x[2] - x[0]) >> 1; \
1149 } \
1150 \
1151 dst[s] += (((int)(((LONG_LONG)((int)(((LONG_LONG)((int)(((LONG_LONG)a * src_subpos) >> 15) + b) * src_subpos) >> 15) + c) * src_subpos) >> 15) + x[1]) * volume_fact) >> 14; \
1152 \
1153 src_subpos += dt + inc; \
1154 src_pos += src_subpos >> 15; \
1155 src_subpos &= 0x7FFFl; \
1156 s++; \
1157 }
1158
1159 MAKE_RESAMPLER();
1160
1161 #undef DO_RESAMPLE
1162
1163 }
1164
1165 }
1166
1167 end:
1168
1169 ASSERT(s <= dst_size);
1170
1171 *_src_pos = src_pos;
1172 *_src_subpos = src_subpos;
1173 if (_dir) *_dir = dir;
1174
1175 return s;
1176 }
1177 #endif
1178