path: root/apps/codecs/dumb/src/helpers/resample.c

/*  _______         ____    __         ___    ___
 * \    _  \       \    /  \  /       \   \  /   /       '   '  '
 *  |  | \  \       |  |    ||         |   \/   |         .      .
 *  |  |  |  |      |  |    ||         ||\  /|  |
 *  |  |  |  |      |  |    ||         || \/ |  |         '  '  '
 *  |  |  |  |      |  |    ||         ||    |  |         .      .
 *  |  |_/  /        \  \__//          ||    |  |
 * /_______/ynamic    \____/niversal  /__\  /____\usic   /|  .  . ibliotheque
 *                                                      /  \
 *                                                     / .  \
 * resample.c - Resampling helper.                    / / \  \
 *                                                   | <  /   \_
 * By Bob and entheh.                                |  \/ /\   /
 *                                                    \_  /  > /
 * In order to find a good trade-off between            | \ / /
 * speed and accuracy in this code, some tests          |  ' /
 * were carried out regarding the behaviour of           \__/
 * long long ints with gcc. The following code
 * was tested:
 *
 * int a, b, c;
 * c = ((long long)a * b) >> 16;
 *
 * DJGPP GCC Version 3.0.3 generated the following assembly language code for
 * the multiplication and scaling, leaving the 32-bit result in EAX.
 *
 * movl  -8(%ebp), %eax    ; read one int into EAX
 * imull -4(%ebp)          ; multiply by the other; result goes in EDX:EAX
 * shrdl $16, %edx, %eax   ; shift EAX right 16, shifting bits in from EDX
 *
 * Note that a 32*32->64 multiplication is performed, allowing for high
 * accuracy. On the Pentium 2 and above, shrdl takes two cycles (generally),
 * so it is a minor concern when four multiplications are being performed
 * (the cubic resampler). On the Pentium MMX and earlier, it takes four or
 * more cycles, so this method is unsuitable for use in the low-quality
 * resamplers.
 *
 * Since "long long" is a gcc-specific extension, we use LONG_LONG instead,
 * defined in dumb.h. We may investigate later what code MSVC generates, but
 * if it seems too slow then we suggest you use a good compiler.
 *
 * FIXME: these comments are somewhat out of date now.
 */

#include <math.h>
#include "dumb.h"



/* Compile with -DHEAVYDEBUG if you want to make sure the pick-up function is
 * called when it should be. There will be a considerable performance hit,
 * since at least one condition has to be tested for every sample generated.
 */
#ifdef HEAVYDEBUG
#define HEAVYASSERT(cond) ASSERT(cond)
#else
#define HEAVYASSERT(cond)
#endif



//#define MULSC(a, b) ((int)((LONG_LONG)(a) * (b) >> 16))
//#define MULSC(a, b) ((a) * ((b) >> 2) >> 14)
#define MULSC(a, b) ((int)((LONG_LONG)((a) << 4) * ((b) << 12) >> 32))



/* A global variable for controlling resampling quality wherever a local
 * specification doesn't override it. The following values are valid:
 *
 *  0 - DUMB_RQ_ALIASING - fastest
 *  1 - DUMB_RQ_LINEAR
 *  2 - DUMB_RQ_CUBIC    - nicest
 *
 * Values outside the range 0-2 will behave the same as the nearest
 * value within the range.
 */
int dumb_resampling_quality = 2;



void dumb_reset_resampler(DUMB_RESAMPLER *resampler, sample_t *src, long pos, long start, long end)
{
	resampler->src = src;
	resampler->pos = pos;
	resampler->subpos = 0;
	resampler->start = start;
	resampler->end = end;
	resampler->dir = 1;
	resampler->pickup = NULL;
	resampler->pickup_data = NULL;
	resampler->min_quality = 0;
	resampler->max_quality = DUMB_RQ_N_LEVELS - 1;
	resampler->x[2] = resampler->x[1] = resampler->x[0] = 0;
	resampler->overshot = -1;
}



DUMB_RESAMPLER *dumb_start_resampler(sample_t *src, long pos, long start, long end)
{
	DUMB_RESAMPLER *resampler = malloc(sizeof(*resampler));
	if (!resampler) return NULL;
	dumb_reset_resampler(resampler, src, pos, start, end);
	return resampler;
}



/* For convenience, returns nonzero on stop. */
static int process_pickup(DUMB_RESAMPLER *resampler)
{
	if (resampler->overshot < 0) {
		resampler->overshot = 0;
		dumb_resample(resampler, NULL, 2, 0, 1.0f);
		resampler->x[0] = resampler->x[1];
	}

	for (;;) {
		if (resampler->dir < 0) {
			if (resampler->overshot >= 3 && resampler->pos+3 >= resampler->start) resampler->x[0] = resampler->src[resampler->pos+3];
			if (resampler->overshot >= 2 && resampler->pos+2 >= resampler->start) resampler->x[1] = resampler->src[resampler->pos+2];
			if (resampler->overshot >= 1 && resampler->pos+1 >= resampler->start) resampler->x[2] = resampler->src[resampler->pos+1];
			resampler->overshot = resampler->start - resampler->pos - 1;
		} else {
			if (resampler->overshot >= 3 && resampler->pos-3 < resampler->end) resampler->x[0] = resampler->src[resampler->pos-3];
			if (resampler->overshot >= 2 && resampler->pos-2 < resampler->end) resampler->x[1] = resampler->src[resampler->pos-2];
			if (resampler->overshot >= 1 && resampler->pos-1 < resampler->end) resampler->x[2] = resampler->src[resampler->pos-1];
			resampler->overshot = resampler->pos - resampler->end;
		}

		if (resampler->overshot < 0) {
			resampler->overshot = 0;
			return 0;
		}

		if (!resampler->pickup) {
			resampler->dir = 0;
			return 1;
		}
		(*resampler->pickup)(resampler, resampler->pickup_data);
		if (resampler->dir == 0) return 1;
		ASSERT(resampler->dir == -1 || resampler->dir == 1);
	}
}



/* Executes the content 'iterator' times.
 * Clobbers the 'iterator' variable.
 * The loop is unrolled by four.
 */
#define LOOP4(iterator, CONTENT) \
{ \
	if ((iterator) & 2) { \
		CONTENT; \
		CONTENT; \
	} \
	if ((iterator) & 1) { \
		CONTENT; \
	} \
	(iterator) >>= 2; \
	while (iterator) { \
		CONTENT; \
		CONTENT; \
		CONTENT; \
		CONTENT; \
		(iterator)--; \
	} \
}



long dumb_resample(DUMB_RESAMPLER *resampler, sample_t *dst, long dst_size, float volume, float delta)
{
	int dt;
	int vol;
	long done;
	long todo;
	int quality;

	if (!resampler || resampler->dir == 0) return 0;
	ASSERT(resampler->dir == -1 || resampler->dir == 1);

	done = 0;
	dt = (int)(delta * 65536.0 + 0.5);
	vol = (int)floor(volume * 65536.0 + 0.5);

	if (vol == 0) dst = NULL;

	quality = dumb_resampling_quality;
	if (quality > resampler->max_quality) quality = resampler->max_quality;
	else if (quality < resampler->min_quality) quality = resampler->min_quality;

	while (done < dst_size) {
		if (process_pickup(resampler)) return done;

		if ((resampler->dir ^ dt) < 0)
			dt = -dt;

		if (resampler->dir < 0)
			todo = (long)((((LONG_LONG)(resampler->pos - resampler->start) << 16) + resampler->subpos - dt) / -dt);
		else
			todo = (long)((((LONG_LONG)(resampler->end - resampler->pos) << 16) - resampler->subpos - 1 + dt) / dt);

		if (todo < 0)
			todo = 0;
		else if (todo > dst_size - done)
			todo = dst_size - done;

		done += todo;

		{
			sample_t *src = resampler->src;
			long pos = resampler->pos;
			int subpos = resampler->subpos;
			long diff = pos;
			long overshot;
			if (resampler->dir < 0) {
				if (!dst) {
					/* Silence or simulation */
					LONG_LONG new_subpos = subpos + dt * todo;
					pos += (long)(new_subpos >> 16);
					subpos = (long)new_subpos & 65535;
				} else if (quality <= DUMB_RQ_ALIASING) {
					/* Aliasing, backwards */
					sample_t xbuf[2];
					sample_t *x = &xbuf[0];
					sample_t *xstart;
					xbuf[0] = resampler->x[1];
					xbuf[1] = resampler->x[2];
					while (todo && x < &xbuf[2]) {
						HEAVYASSERT(pos >= resampler->start);
						*dst++ += MULSC(x[0], vol);
						subpos += dt;
						pos += subpos >> 16;
						x -= subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = xstart = &src[pos];
					LOOP4(todo,
						*dst++ += MULSC(x[2], vol);
						subpos += dt;
						x += subpos >> 16;
						subpos &= 65535;
					);
					pos += x - xstart;
				} else if (quality <= DUMB_RQ_LINEAR) {
					/* Linear interpolation, backwards */
					sample_t xbuf[3];
					sample_t *x = &xbuf[1];
					xbuf[0] = resampler->x[1];
					xbuf[1] = resampler->x[2];
					xbuf[2] = src[pos];
					while (todo && x < &xbuf[3]) {
						HEAVYASSERT(pos >= resampler->start);
						*dst++ += MULSC(x[0] + MULSC(x[-1] - x[0], subpos), vol);
						subpos += dt;
						pos += subpos >> 16;
						x -= subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = &src[pos];
					LOOP4(todo,
						HEAVYASSERT(pos >= resampler->start);
						*dst++ += MULSC(x[1] + MULSC(x[2] - x[1], subpos), vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
					);
				} else {
					/* Cubic interpolation, backwards */
					sample_t xbuf[6];
					sample_t *x = &xbuf[3];
					sample_t *lastx = NULL;
					int a = 0, b = 0, c = 0;
					xbuf[0] = resampler->x[0];
					xbuf[1] = resampler->x[1];
					xbuf[2] = resampler->x[2];
					xbuf[3] = src[pos];
					if (pos-1 >= resampler->start) xbuf[4] = src[pos-1];
					if (pos-2 >= resampler->start) xbuf[5] = src[pos-2];
					while (todo && x < &xbuf[6]) {
						HEAVYASSERT(pos >= resampler->start);
						if (lastx != x) {
							lastx = x;
							a = (((x[-1] - x[-2]) << 1) + (x[-1] - x[-2]) + (x[-3] - x[0])) >> 1;
							b = (x[-2] << 1) + x[0] - ((5 * x[-1] + x[-3]) >> 1);
							c = (x[-2] - x[0]) >> 1;
						}
						*dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-1], vol);
						subpos += dt;
						pos += subpos >> 16;
						x -= subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = &src[pos];
					lastx = NULL;
					LOOP4(todo,
						HEAVYASSERT(pos >= resampler->start);
						if (lastx != x) {
							lastx = x;
							a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
							b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
							c = (x[2] - x[0]) >> 1;
						}
						*dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[1], vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
					);
				}
				diff = diff - pos;
				overshot = resampler->start - pos - 1;
				if (diff >= 3) {
					resampler->x[0] = overshot >= 3 ? 0 : src[pos+3];
					resampler->x[1] = overshot >= 2 ? 0 : src[pos+2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
				} else if (diff >= 2) {
					resampler->x[0] = resampler->x[2];
					resampler->x[1] = overshot >= 2 ? 0 : src[pos+2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
				} else if (diff >= 1) {
					resampler->x[0] = resampler->x[1];
					resampler->x[1] = resampler->x[2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos+1];
				}
			} else {
				if (!dst) {
					/* Silence or simulation */
					LONG_LONG new_subpos = subpos + dt * todo;
					pos += (long)(new_subpos >> 16);
					subpos = (long)new_subpos & 65535;
				} else if (dumb_resampling_quality <= DUMB_RQ_ALIASING) {
					/* Aliasing, forwards */
					sample_t xbuf[2];
					sample_t *x = &xbuf[0];
					sample_t *xstart;
					xbuf[0] = resampler->x[1];
					xbuf[1] = resampler->x[2];
					while (todo && x < &xbuf[2]) {
						HEAVYASSERT(pos < resampler->end);
						*dst++ += MULSC(x[0], vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = xstart = &src[pos];
					LOOP4(todo,
						*dst++ += MULSC(x[-2], vol);
						subpos += dt;
						x += subpos >> 16;
						subpos &= 65535;
					);
					pos += x - xstart;
				} else if (dumb_resampling_quality <= DUMB_RQ_LINEAR) {
					/* Linear interpolation, forwards */
					sample_t xbuf[3];
					sample_t *x = &xbuf[1];
					xbuf[0] = resampler->x[1];
					xbuf[1] = resampler->x[2];
					xbuf[2] = src[pos];
					while (todo && x < &xbuf[3]) {
						HEAVYASSERT(pos < resampler->end);
						*dst++ += MULSC(x[-1] + MULSC(x[0] - x[-1], subpos), vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = &src[pos];
					LOOP4(todo,
						HEAVYASSERT(pos < resampler->end);
						*dst++ += MULSC(x[-2] + MULSC(x[-1] - x[-2], subpos), vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
					);
				} else {
					/* Cubic interpolation, forwards */
					sample_t xbuf[6];
					sample_t *x = &xbuf[3];
					sample_t *lastx = NULL;
					int a = 0, b = 0, c = 0;
					xbuf[0] = resampler->x[0];
					xbuf[1] = resampler->x[1];
					xbuf[2] = resampler->x[2];
					xbuf[3] = src[pos];
					if (pos+1 < resampler->end) xbuf[4] = src[pos+1];
					if (pos+2 < resampler->end) xbuf[5] = src[pos+2];
					while (todo && x < &xbuf[6]) {
						HEAVYASSERT(pos < resampler->end);
						if (lastx != x) {
							lastx = x;
							a = (((x[-2] - x[-1]) << 1) + (x[-2] - x[-1]) + (x[0] - x[-3])) >> 1;
							b = (x[-1] << 1) + x[-3] - ((5 * x[-2] + x[0]) >> 1);
							c = (x[-1] - x[-3]) >> 1;
						}
						*dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-2], vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
						todo--;
					}
					x = &src[pos];
					lastx = NULL;
					LOOP4(todo,
						HEAVYASSERT(pos < resampler->end);
						if (lastx != x) {
							lastx = x;
							a = (((x[-2] - x[-1]) << 1) + (x[-2] - x[-1]) + (x[0] - x[-3])) >> 1;
							b = (x[-1] << 1) + x[-3] - ((5 * x[-2] + x[0]) >> 1);
							c = (x[-1] - x[-3]) >> 1;
						}
						*dst++ += MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[-2], vol);
						subpos += dt;
						pos += subpos >> 16;
						x += subpos >> 16;
						subpos &= 65535;
					);
				}
				diff = pos - diff;
				overshot = pos - resampler->end;
				if (diff >= 3) {
					resampler->x[0] = overshot >= 3 ? 0 : src[pos-3];
					resampler->x[1] = overshot >= 2 ? 0 : src[pos-2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
				} else if (diff >= 2) {
					resampler->x[0] = resampler->x[2];
					resampler->x[1] = overshot >= 2 ? 0 : src[pos-2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
				} else if (diff >= 1) {
					resampler->x[0] = resampler->x[1];
					resampler->x[1] = resampler->x[2];
					resampler->x[2] = overshot >= 1 ? 0 : src[pos-1];
				}
			}
			resampler->pos = pos;
			resampler->subpos = subpos;
		}
	}

	return done;
}



sample_t dumb_resample_get_current_sample(DUMB_RESAMPLER *resampler, float volume)
{
	int vol;
	sample_t *src;
	long pos;
	int subpos;
	int quality;

	if (!resampler || resampler->dir == 0) return 0;
	ASSERT(resampler->dir == -1 || resampler->dir == 1);

	if (process_pickup(resampler)) return 0;

	vol = (int)floor(volume * 65536.0 + 0.5);
	if (vol == 0) return 0;

	quality = dumb_resampling_quality;
	if (quality > resampler->max_quality) quality = resampler->max_quality;
	else if (quality < resampler->min_quality) quality = resampler->min_quality;

	src = resampler->src;
	pos = resampler->pos;
	subpos = resampler->subpos;

	if (resampler->dir < 0) {
		HEAVYASSERT(pos >= resampler->start);
		if (dumb_resampling_quality <= 0) {
			/* Aliasing, backwards */
			return MULSC(src[pos], vol);
		} else if (quality <= DUMB_RQ_LINEAR) {
			/* Linear interpolation, backwards */
			return MULSC(resampler->x[2] + MULSC(resampler->x[1] - resampler->x[2], subpos), vol);
		} else {
			/* Cubic interpolation, backwards */
			sample_t *x = resampler->x;
			int a, b, c;
			a = (((x[2] - x[1]) << 1) + (x[2] - x[1]) + (x[0] - src[pos])) >> 1;
			b = (x[1] << 1) + src[pos] - ((5 * x[2] + x[0]) >> 1);
			c = (x[1] - src[pos]) >> 1;
			return MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[2], vol);
		}
	} else {
		HEAVYASSERT(pos < resampler->end);
		if (dumb_resampling_quality <= 0) {
			/* Aliasing */
			return MULSC(src[pos], vol);
		} else if (dumb_resampling_quality <= DUMB_RQ_LINEAR) {
			/* Linear interpolation, forwards */
			return MULSC(resampler->x[1] + MULSC(resampler->x[2] - resampler->x[1], subpos), vol);
		} else {
			/* Cubic interpolation, forwards */
			sample_t *x = resampler->x;
			int a, b, c;
			a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (src[pos] - x[0])) >> 1;
			b = (x[2] << 1) + x[0] - ((5 * x[1] + src[pos]) >> 1);
			c = (x[2] - x[0]) >> 1;
			return MULSC(MULSC(MULSC(MULSC(a, subpos) + b, subpos) + c, subpos) + x[1], vol);
		}
	}
}



void dumb_end_resampler(DUMB_RESAMPLER *resampler)
{
	if (resampler)
		free(resampler);
}



#if 0
/* The following macro is used to overcome the fact that most C
 * compilers (including gcc and MSVC) can't correctly multiply signed
 * integers outside the range -32768 to 32767. i86 assembler versions
 * don't need to use this method, since the processor does in fact
 * have instructions to multiply large numbers correctly - which
 * means using assembly language could make a significant difference
 * to the speed.
 *
 * The basic method is as follows. We halve the subposition (how far
 * we are between samples), so it never exceeds 32767. We also halve
 * the delta, which is the amount to be added to the subposition each
 * time. Then we unroll the loop twofold, so that we can add the lost
 * one every other time if necessary (since the halving may have
 * resulted in rounding down).
 *
 * This method doesn't incur any cumulative inaccuracies. There is a
 * very slight loss of quality, which I challenge anyone to notice -
 * but the position will advance at *exactly* the same rate as it
 * would if we didn't use this method. This also means the pitch is
 * exactly the same, which may even make a difference to trained
 * musicians when resampling down a lot :)
 *
 * Each time this macro is invoked, DO_RESAMPLE(inc) must be defined
 * to calculate the samples by the appropriate equation (linear,
 * cubic, etc.). See the individual cases for examples of how this is
 * done.
 */
#define MAKE_RESAMPLER()							\
{													\
	if (dt & 1) {									\
		long todo2;									\
													\
		dt >>= 1;									\
													\
		if (src_subpos & 1) {						\
			src_subpos >>= 1;						\
			DO_RESAMPLE(1);							\
			todo--;									\
		} else										\
			src_subpos >>= 1;						\
													\
		todo2 = todo >> 1;							\
													\
		while (todo2) {								\
			DO_RESAMPLE(0);							\
			DO_RESAMPLE(1);							\
			todo2--;								\
		}											\
													\
		if (todo & 1) {								\
			DO_RESAMPLE(0);							\
			src_subpos = (src_subpos << 1) | 1;		\
		} else										\
			src_subpos <<= 1;						\
													\
		todo = 0;									\
		dt = (dt << 1) | 1;							\
	} else {										\
		long subposbit = src_subpos & 1;			\
		dt >>= 1;									\
		src_subpos >>= 1;							\
													\
		if (todo & 1) {								\
			DO_RESAMPLE(0);							\
		}											\
													\
		todo >>= 1;									\
													\
		while (todo) {								\
			DO_RESAMPLE(0);							\
			DO_RESAMPLE(0);							\
			todo--;									\
		}											\
													\
		src_subpos = (src_subpos << 1) | subposbit; \
		dt <<= 1;									\
	}												\
}



sample_t dumb_resample_get_current_sample(
	sample_t *src, long *_src_pos, int *_src_subpos,
	long src_start, long src_end,
	float volume, int *_dir,
	DUMB_RESAMPLE_PICKUP pickup, void *pickup_data
)
{
	long src_pos = *_src_pos;
	int src_subpos = *_src_subpos;
	int dir = _dir ? *_dir : 1;

	sample_t value = 0;

	if (dir == 0)
		return 0;

	ASSERT(dir == 1 || dir == -1);

	if (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end)) {

		/* If there's no pick-up function, we stop. */
		if (!pickup) {
			dir = 0;
			goto end;
		}

		/* Process the pick-up. It may need invoking more than once. */
		do {
			dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);

			if (dir == 0)
				goto end;

			ASSERT(dir == 1 || dir == -1);
		} while (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end));
	}

	HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end));

	if (dumb_resampling_quality == 0) {
		/* Aliasing (coarse) */
		int volume_fact = (int)(volume * 16384.0);
		value = (src[src_pos] * volume_fact) >> 14;
	} else if (dumb_resampling_quality <= 2) {
		/* Linear interpolation */
		int volume_fact = (int)(volume * 16384.0);
		int subpos = src_subpos >> 1;
		value = ((src[src_pos] + ((((src[src_pos + 1] - src[src_pos]) >> 1) * subpos) >> 14)) * volume_fact) >> 14;
	} else if (dumb_resampling_quality == 3) {
		/* Quadratic interpolation */
		int volume_fact = (int)(volume * 16384.0);
		int a, b;
		sample_t *x;
		int subpos = src_subpos >> 1;
		x = &src[src_pos];
		a = ((x[0] + x[2]) >> 1) - x[1];
		b = ((x[2] - x[0]) >> 1) - (a << 1);
		value = (((((((a * subpos) >> 15) + b) * subpos) >> 15) + x[0]) * volume_fact) >> 14;
	} else {
		/* Cubic interpolation */
		int volume_fact = (int)(volume * 16384.0);
		int a, b, c;
		sample_t *x;
		int subpos = src_subpos >> 1;
		x = &src[src_pos];
		a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
		b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
		c = (x[2] - x[0]) >> 1;
		value = (((int)(((LONG_LONG)((int)(((LONG_LONG)((int)(((LONG_LONG)a * subpos) >> 15) + b) * subpos) >> 15) + c) * subpos) >> 15) + x[1]) * volume_fact) >> 14;
	}

	end:

	*_src_pos = src_pos;
	*_src_subpos = src_subpos;
	if (_dir) *_dir = dir;

	return value;
}



long dumb_resample(
	sample_t *src, long *_src_pos, int *_src_subpos,
	long src_start, long src_end,
	sample_t *dst, long dst_size,
	float volume, float delta, int *_dir,
	DUMB_RESAMPLE_PICKUP pickup, void *pickup_data
)
{
	int dt = (int)(delta * 65536.0 + 0.5);
	long s = 0; /* Current position in the destination buffer */

	long src_pos = *_src_pos;
	int src_subpos = *_src_subpos;
	int dir = _dir ? *_dir : 1;

	int linear_average;

	if (dir == 0)
		return 0;

	ASSERT(dir == 1 || dir == -1);

	linear_average = dst && dumb_resampling_quality >= 2 && dt > 65536;

	if (dir < 0) dt = -dt;

	if (linear_average)
		volume /= delta;

	while (s < dst_size) {

		long todo;

		/* Process pick-ups first, just in case. */

		if (linear_average) {

			/* For linear average, the pick-up point could split a sum into
			 * two parts. We handle this by putting the pick-up code inside
			 * the summing loop. Note that this code is only executed when we
			 * know that a pick-up is necessary somewhere during this sum
			 * (although it is always executed once for the first sample).
			 * We use a separate loop further down when we know we won't have
			 * to do a pick-up, so the condition does not need testing inside
			 * the loop.
			 */

			float sum;
			long i;
			int advance;
			int x[3];

			advance = src_subpos + dt;

			/* Make these negative. Then they stay within the necessary
			 * range for integer multiplication, -32768 to 32767 ;)
			 */
			x[0] = ~(src_subpos >> 1); /* = -1 - (src_subpos >> 1) */
			x[2] = x[0] ^ 0x7FFF; /* = -32768 + (src_subpos >> 1) */

			sum = (float)(-((src[src_pos] * (x+1)[dir]) >> 15));

			i = src_pos + (advance >> 16);
			src_pos += dir;
			src_subpos = (dir >> 1) & 65535; /* changes 1,-1 to 0,65535 */

			advance &= 65535;

			/* i is the index of the first sample NOT to sum fully,
			 * regardless of the direction of resampling.
			 */

			while (dir < 0 ? (i < src_start) : (i >= src_end)) {
				if (dir < 0) {
					while (src_pos >= src_start)
						sum += src[src_pos--];
				} else {
					while (src_pos < src_end)
						sum += src[src_pos++];
				}

				i -= src_pos;
				/* i is now the number of samples left to sum fully, except
				 * it's negative if we're going backwards.
				 */

				if (!pickup) {
					dir = 0;
					goto endsum;
				}

				dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);

				if (dir == 0)
					goto endsum;

				ASSERT(dir == 1 || dir == -1);

				if ((dir ^ dt) < 0) {
					dt = -dt;
					advance ^= 65535;
					i = -i;
				}

				i += src_pos;
				/* There, i is back to normal. */
			}

			for (; src_pos != i; src_pos += dir)
				sum += src[src_pos];

			src_subpos = advance;

			x[2] = src_subpos >> 1;
			x[0] = x[2] ^ 0x7FFF; /* = 32767 - (src_subpos >> 1) */

			sum += (src[src_pos] * (x+1)[dir]) >> 15;

			endsum:

			sum *= volume;
			dst[s] += (int)sum;

			s++;

			if (dir == 0)
				break;

		} else if (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end)) {

			/* If there's no pick-up function, we stop. */
			if (!pickup) {
				dir = 0;
				break;
			}

			/* Process the pick-up. It may need invoking more than once. */
			do {
				dir = (*pickup)(src, &src_pos, &src_subpos, &src_start, &src_end, dir, pickup_data);

				if (dir == 0)
					goto end;

				ASSERT(dir == 1 || dir == -1);
			} while (dir < 0 ? (src_pos < src_start) : (src_pos >= src_end));

			/* Update sign of dt to match that of dir. */
			if ((dir ^ dt) < 0)
				dt = -dt;
		}

		/* Work out how many contiguous samples we can now render. */
		if (dir < 0)
			todo = (long)((((LONG_LONG)(src_pos - src_start) << 16) + src_subpos) / -dt);
		else
			todo = (long)((((LONG_LONG)(src_end - src_pos) << 16) - src_subpos - 1) / dt);

		/* The above equations work out how many complete dt-sized
		 * intervals there are between the current position and the loop
		 * point (provided there is a little fractional extra). The linear
		 * average function needs complete intervals - but the other
		 * resamplers only read a sample from the beginning of each interval,
		 * so they can process one extra sample in their main loops (so we
		 * increment todo in a moment).
		 *
		 * The linear average function makes up the extra sample using the
		 * specialised pick-up code above.
		 *
		 * Note that our above pick-up process should have absolutely ensured
		 * that the result of this function will be nonnegative.
		 */

		ASSERT(todo >= 0);

		if (!linear_average)
			todo++;

		/* Of course we don't want to overrun the output buffer! */
		if (todo > dst_size - s)
			todo = dst_size - s;

		if (!dst) {

			LONG_LONG t = src_subpos + (LONG_LONG)dt * todo;
			src_pos += (long)(t >> 16);
			src_subpos = (int)t & 0xFFFFl;

			s += todo;

		} else if (linear_average) {

			float sum;
			long i;
			int advance;
			int x[3];

			while (todo) {

				advance = src_subpos + dt;

				/* Make these negative. Then they stay within the necessary
				 * range for integer multiplication, -32768 to 32767 ;)
				 */
				x[0] = ~(src_subpos >> 1); /* = -1 - (src_subpos >> 1) */
				x[2] = x[0] ^ 0x7FFF; /* = -32768 + (src_subpos >> 1) */

				sum = (float)(-((src[src_pos] * (x+1)[dir]) >> 15));

				i = src_pos + (advance >> 16);
				src_pos += dir;
				src_subpos = (dir >> 1) & 65535; /* changes 1,-1 to 0,65535 */

				advance &= 65535;

				/* i is the index of the first sample NOT to sum fully,
				 * regardless of the direction of resampling.
				 */

				HEAVYASSERT(dir < 0 ? (i >= src_start) : (i < src_end));

				for (; src_pos != i; src_pos += dir)
					sum += src[src_pos];

				src_subpos = advance;

				x[2] = src_subpos >> 1;
				x[0] = x[2] ^ 0x7FFF; /* = 32767 - (src_subpos >> 1) */

				sum += (src[src_pos] * (x+1)[dir]) >> 15;

				sum *= volume;
				dst[s] += (int)sum;

				s++;
				todo--;
			}

		} else if (dumb_resampling_quality == 0 || (dumb_resampling_quality == 1 && delta >= 1.0)) {

			/* Aliasing (coarse) */
			int volume_fact = (int)(volume * 16384.0);

			do {
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end));
				dst[s] += ((src[src_pos] * volume_fact) >> 14);
				src_subpos += dt;
				src_pos += src_subpos >> 16;
				src_subpos &= 0xFFFFl;
				s++;
			} while (--todo);

		} else if (dumb_resampling_quality <= 2) {

			/* Linear interpolation */
			int volume_fact = (int)(volume * 16384.0);

			#define DO_RESAMPLE(inc)		 \
			{								 \
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
											 \
				dst[s] += (((src[src_pos] + ((((src[src_pos + 1] - src[src_pos]) >> 1) * src_subpos) >> 14)) * volume_fact) >> 14); \
											 \
				src_subpos += dt + inc;		 \
				src_pos += src_subpos >> 15; \
				src_subpos &= 0x7FFFl;		 \
				s++;						 \
			}

			MAKE_RESAMPLER();

			#undef DO_RESAMPLE

		} else if (dumb_resampling_quality == 3) {

			/* Quadratic interpolation */

			int volume_fact = (int)(volume * 16384.0);
			int a = 0, b = 0;
			sample_t *x = NULL;
			int last_src_pos = -1;

			/* AIM: no integer multiplicands must transcend the range -32768 to 32767.
			 * This limitation is imposed by most compilers, including gcc and MSVC.
			 *
			 * a = 0.5 * (s0 + s2) - s1
			 * b = -1.5 * s0 + 2 * s1 - 0.5 * s2
			 * c = s0
			 *
			 * s = (a * t + b) * t + c
			 *
			 * In fixed-point:
			 *
			 * a = ((s0 + s2) >> 1) - s1
			 * b = ((-3 * s0 - s2) >> 1) + (s1 << 1)
			 *
			 * s = (((((a * t) >> 16) + b) * t) >> 16) + s0
			 *
			 * With t halved (since t can reach 65535):
			 *
			 * s = (((((a * t) >> 15) + b) * t) >> 15) + s0
			 *
			 * a currently reaches 65536
			 * b currently reaches 131072
			 *
			 * So we must use aon2
			 *
			 * s = (((((aon2 * t) >> 14) + b) * t) >> 15) + s0
			 *
			 * ((aon2 * t) >> 14) + b is 5 times too big
			 * so we must divide by 8
			 *
			 * s = (((((aon2 * t) >> 17) + bon8) * t) >> 12) + s0
			 *
			 * aon2 = ((s0 + s2) >> 2) - (s1 >> 1)
			 * bon8 = ((-3 * s0 - s2) >> 4) + (s1 >> 2)
			 * or:
			 * bon8 = ((s2 - s0) >> 4) - (aon2 >> 1)
			 */

			/* Unh4x0r3d version:
			#define DO_RESAMPLE(inc)						\
			{												\
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
															\
				if (src_pos != last_src_pos) {				\
					last_src_pos = src_pos;					\
					x = &src[src_pos];						\
					a = ((x[0] + x[2]) >> 2) - (x[1] >> 1); \
					b = ((x[2] - x[0]) >> 4) - (a >> 1);	\
				}											\
															\
				dst[s] += ((((((((a * src_subpos) >> 17) + b) * src_subpos) >> 12) + x[0]) * volume_fact) >> 14); \
															\
				src_subpos += dt + inc;						\
				src_pos += src_subpos >> 15;				\
				src_subpos &= 0x7FFFl;						\
				s++;										\
			}
			*/

			/* H4x0r3d version: */
			#define DO_RESAMPLE(inc)						\
			{												\
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
															\
				if (src_pos != last_src_pos) {				\
					last_src_pos = src_pos;					\
					x = &src[src_pos];						\
					a = ((x[0] + x[2]) >> 1) - x[1];		\
					b = ((x[2] - x[0]) >> 1) - (a << 1);	\
				}											\
															\
				dst[s] += ((((((((a * src_subpos) >> 15) + b) * src_subpos) >> 15) + x[0]) * volume_fact) >> 14); \
															\
				src_subpos += dt + inc;						\
				src_pos += src_subpos >> 15;				\
				src_subpos &= 0x7FFFl;						\
				s++;										\
			}

			MAKE_RESAMPLER();

			#undef DO_RESAMPLE

		} else {

			/* Cubic interpolation */

			int volume_fact = (int)(volume * 16384.0);
			int a = 0, b = 0, c = 0;
			sample_t *x = NULL;
			int last_src_pos = -1;

			/* AIM: never multiply integers outside the range -32768 to 32767.
			 *
			 * a = 1.5f * (x[1] - x[2]) + (x[3] - x[0]) * 0.5f;
			 * b = 2.0f * x[2] + x[0] - 2.5f * x[1] - x[3] * 0.5f;
			 * c = (x[2] - x[0]) * 0.5f;
			 *
			 * s = ((a * t + b) * t + c) * t + x[1];
			 *
			 * Fixed-point version:
			 *
			 * a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1;
			 * b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1);
			 * c = (x[2] - x[0]) >> 1;
			 *
			 * s = ((((((((a * t) >> 15) + b) * t) >> 15) + c) * t) >> 15) + x[1];
			 *   (with t already halved, maximum 32767)
			 *
			 * a is in (((1+1)*2)+(1+1)+(1+1))/2 = 8 times the required range
			 * b is in (1*2)+1+((5*1+1)/2) = 6 times
			 * c is in the required range
			 *
			 * We must use aon8
			 *
			 * s = ((((((((aon8 * t) >> 12) + b) * t) >> 15) + c) * t) >> 15) + x[1];
			 *
			 * But ((aon8 * t) >> 12) is in 2^(15+15-12) = 2^18 = 8 times
			 * b is in 6 times
			 * so we divide both ((aon8 * t) >> 12) and b by 16
			 *
			 * s = ((((((((aon8 * t) >> 16) + bon16) * t) >> 11) + c) * t) >> 15) + x[1];
			 *
			 * ((... + bon16) * t) >> 11 is 16 times too big
			 * c is in the correct range
			 * we must divide both by 32
			 *
			 * s = ((((((((aon8 * t) >> 16) + bon16) * t) >> 16) + con32) * t) >> 10) + x[1];
			 *
			 * aon8  = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 4;
			 * bon16 = ((x[2] << 2) + (x[0] << 1) - (5 * x[1] + x[3])) >> 5;
			 * con32 = (x[2] - x[0]) >> 6;
			 *
			 * A lot of accuracy is lost here. It is quite likely that some
			 * of the above would cancel anyway, so the scaling down wouldn't
			 * have to be so severe. However, I'm not in the mood to work it
			 * out now :P
			 *
			 * It may also be worth investigating whether doing this stuff
			 * in floats would be faster.
			 */

			/* Unh4x0r3d version:
			#define DO_RESAMPLE(inc)		   \
			{								   \
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
											   \
				if (src_pos != last_src_pos) { \
					last_src_pos = src_pos;	   \
					x = &src[src_pos];		   \
					a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 4; \
					b = ((x[2] << 2) + (x[0] << 1) - (5 * x[1] + x[3])) >> 5; \
					c = (x[2] - x[0]) >> 6;	   \
				}							   \
											   \
				dst[s] += ((((((((((a * src_subpos) >> 16) + b) * src_subpos) >> 16) + c) * src_subpos) >> 10) + x[1]) * volume_fact) >> 14; \
											   \
				src_subpos += dt + inc;		   \
				src_pos += src_subpos >> 15;   \
				src_subpos &= 0x7FFFl;		   \
				s++;						   \
			}
			*/

			/* H4x0r3d version: */
			#define DO_RESAMPLE(inc)		   \
			{								   \
				HEAVYASSERT(dir < 0 ? (src_pos >= src_start) : (src_pos < src_end)); \
											   \
				if (src_pos != last_src_pos) { \
					last_src_pos = src_pos;	   \
					x = &src[src_pos];		   \
					a = (((x[1] - x[2]) << 1) + (x[1] - x[2]) + (x[3] - x[0])) >> 1; \
					b = (x[2] << 1) + x[0] - ((5 * x[1] + x[3]) >> 1); \
					c = (x[2] - x[0]) >> 1;	   \
				}							   \
											   \
				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; \
											   \
				src_subpos += dt + inc;		   \
				src_pos += src_subpos >> 15;   \
				src_subpos &= 0x7FFFl;		   \
				s++;						   \
			}

			MAKE_RESAMPLER();

			#undef DO_RESAMPLE

		}

	}

	end:

	ASSERT(s <= dst_size);

	*_src_pos = src_pos;
	*_src_subpos = src_subpos;
	if (_dir) *_dir = dir;

	return s;
}
#endif