new files, based on SoX (http://sox.sf.net): better resampling code. Note that my mixer.cpp changes are on purpose not yet in CVS since they are not complete. Only reasons I checkin these files is that it's much more comfortable to have CVS, since I need to rewrite parts of resample.cpp now (I already have lots of modifications in). Also expect more OO in the future

svn-id: r9176

5 files changed, 1101 insertions, 0 deletions

diff --git a/sound/audiostream.h b/sound/audiostream.h
new file mode 100644
index 0000000000..c3d43d6b2a
--- /dev/null
+++ b/sound/audiostream.h
@@ -0,0 +1,91 @@
+/* ScummVM - Scumm Interpreter
+ * Copyright (C) 2001-2003 The ScummVM project
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * $Header$
+ *
+ */
+
+#ifndef AUDIOSTREAM_H
+#define AUDIOSTREAM_H
+
+/**
+ * Generic input stream for the resampling code.
+ */
+class InputStream {
+public:
+	byte *_ptr;
+	byte *_end;
+	InputStream(byte *ptr, uint len) : _ptr(ptr), _end(ptr+len) { }
+	virtual int16 readIntern() = 0;
+	virtual void advance() = 0;
+public:
+	int16 read() { assert(size() > 0); int16 val = readIntern(); advance(); return val; }
+	int16 peek() { assert(size() > 0); return readIntern(); }
+	virtual int size() = 0;
+	bool eof() { return size() <= 0; }
+};
+
+class ZeroInputStream : public InputStream {
+protected:
+	int16 readIntern() { return 0; }
+	void advance() { _ptr++; }
+public:
+	ZeroInputStream(uint len) : InputStream(0, len) { }
+	virtual int size() { return _end - _ptr; }
+};
+
+template<int channels>
+class Input8bitSignedStream : public InputStream {
+protected:
+	int16 readIntern() { int8 v = (int8)*_ptr; return v << 8; }
+	void advance() { _ptr += channels; }
+public:
+	Input8bitSignedStream(byte *ptr, int len) : InputStream(ptr, len) { }
+	virtual int size() { return (_end - _ptr) / channels; }
+};
+
+template<int channels>
+class Input8bitUnsignedStream : public InputStream {
+protected:
+	int16 readIntern() { int8 v = (int8)(*_ptr ^ 0x80); return v << 8; }
+	void advance() { _ptr += channels; }
+public:
+	Input8bitUnsignedStream(byte *ptr, int len) : InputStream(ptr, len) { }
+	virtual int size() { return (_end - _ptr) / channels; }
+};
+
+template<int channels>
+class Input16bitSignedStream : public InputStream {
+protected:
+	int16 readIntern() { return (int16)READ_BE_UINT16(_ptr); }
+	void advance() { _ptr += 2*channels; }
+public:
+	Input16bitSignedStream(byte *ptr, int len) : InputStream(ptr, len) { }
+	virtual int size() { return (_end - _ptr) / (2 * channels); }
+};
+
+template<int channels>
+class Input16bitUnsignedStream : public InputStream {
+protected:
+	int16 readIntern() { return (int16)(READ_BE_UINT16(_ptr) ^ 0x8000); }
+	void advance() { _ptr += 2*channels; }
+public:
+	Input16bitUnsignedStream(byte *ptr, int len) : InputStream(ptr, len) { }
+	virtual int size() { return (_end - _ptr) / (2 * channels); }
+};
+
+#endif
diff --git a/sound/rate.cpp b/sound/rate.cpp
new file mode 100644
index 0000000000..bde477647a
--- /dev/null
+++ b/sound/rate.cpp
@@ -0,0 +1,164 @@
+/*
+* August 21, 1998
+* Copyright 1998 Fabrice Bellard.
+*
+* [Rewrote completly the code of Lance Norskog And Sundry
+* Contributors with a more efficient algorithm.]
+*
+* This source code is freely redistributable and may be used for
+* any purpose.	 This copyright notice must be maintained. 
+* Lance Norskog And Sundry Contributors are not responsible for 
+* the consequences of using this software.  
+*/
+
+/*
+ * Sound Tools rate change effect file.
+ */
+
+#include "rate.h"
+
+#include <math.h>
+/*
+ * Linear Interpolation.
+ *
+ * The use of fractional increment allows us to use no buffer. It
+ * avoid the problems at the end of the buffer we had with the old
+ * method which stored a possibly big buffer of size
+ * lcm(in_rate,out_rate).
+ *
+ * Limited to 16 bit samples and sampling frequency <= 65535 Hz. If
+ * the input & output frequencies are equal, a delay of one sample is
+ * introduced.	Limited to processing 32-bit count worth of samples.
+ *
+ * 1 << FRAC_BITS evaluating to zero in several places.	 Changed with
+ * an (unsigned long) cast to make it safe.  MarkMLl 2/1/99
+ *
+ * Replaced all uses of floating point arithmetic by fixed point
+ * calculations (Max Horn 2003-07-18).
+ */
+
+#define FRAC_BITS 16
+
+/* Private data */
+
+typedef struct ratestuff
+{
+	unsigned long opos_frac;  /* fractional position of the output stream in input stream unit */
+	unsigned long opos;
+
+	unsigned long opos_inc_frac;  /* fractional position increment in the output stream */
+	unsigned long opos_inc;
+
+	unsigned long ipos;	 /* position in the input stream (integer) */
+
+	st_sample_t ilast; /* last sample in the input stream */
+} *rate_t;
+
+/*
+ * Process options
+ */
+int st_rate_getopts(eff_t effp, int n, char **argv)
+{
+	if (n) {
+		st_fail("Rate effect takes no options.");
+		return (ST_EOF);
+	}
+
+	return (ST_SUCCESS);
+}
+
+/*
+ * Prepare processing.
+ */
+int st_rate_start(eff_t effp, st_rate_t inrate, st_rate_t outrate)
+{
+	rate_t rate = (rate_t) effp->priv;
+	unsigned long incr;
+
+	if (inrate == outrate) {
+		st_fail("Input and Output rates must be different to use rate effect");
+		return (ST_EOF);
+	}
+
+	if (inrate >= 65536 || outrate >= 65536) {
+		st_fail("rate effect can only handle rates < 65536");
+		return (ST_EOF);
+	}
+
+	rate->opos_frac = 0;
+	rate->opos = 0;
+
+	/* increment */
+	incr = (inrate << FRAC_BITS) / outrate;
+
+	rate->opos_inc_frac = incr & ((1UL << FRAC_BITS) - 1);
+	rate->opos_inc = incr >> FRAC_BITS;
+
+	rate->ipos = 0;
+
+	rate->ilast = 0;
+	return (ST_SUCCESS);
+}
+
+/*
+ * Processed signed long samples from ibuf to obuf.
+ * Return number of samples processed.
+ */
+int st_rate_flow(eff_t effp, InputStream &input, st_sample_t *obuf, st_size_t *osamp)
+{
+	rate_t rate = (rate_t) effp->priv;
+	st_sample_t *ostart, *oend;
+	st_sample_t ilast, icur, out;
+	unsigned long tmp;
+
+	ilast = rate->ilast;
+
+	ostart = obuf;
+	oend = obuf + *osamp;
+
+	while (obuf < oend && !input.eof()) {
+
+		/* read as many input samples so that ipos > opos */
+		while (rate->ipos <= rate->opos) {
+			ilast = input.read();
+			rate->ipos++;
+			/* See if we finished the input buffer yet */
+
+			if (input.eof())
+				goto the_end;
+		}
+
+		icur = input.peek();
+
+		/* interpolate */
+		out = ilast + (((icur - ilast) * rate->opos_frac) >> FRAC_BITS);
+
+		/* output sample & increment position */
+
+		clampedAdd(*obuf++, out);
+		#if 1	// FIXME: Hack to generate stereo output
+
+		clampedAdd(*obuf++, out);
+		#endif
+
+		tmp = rate->opos_frac + rate->opos_inc_frac;
+		rate->opos = rate->opos + rate->opos_inc + (tmp >> FRAC_BITS);
+		rate->opos_frac = tmp & ((1UL << FRAC_BITS) - 1);
+	}
+
+the_end:
+	*osamp = obuf - ostart;
+	rate->ilast = ilast;
+	return (ST_SUCCESS);
+}
+
+/*
+ * Do anything required when you stop reading samples.	
+ * Don't close input file! 
+ */
+int st_rate_stop(eff_t effp)
+{
+	/* nothing to do */
+	return (ST_SUCCESS);
+}
+
diff --git a/sound/rate.h b/sound/rate.h
new file mode 100644
index 0000000000..aa8edea1b0
--- /dev/null
+++ b/sound/rate.h
@@ -0,0 +1,72 @@
+// HACK: Instead of using the full st_i.h (and then st.h and stconfig.h etc.)
+// from SoX, we use this minimal variant which is just sufficient to make
+// resample.c and rate.c compile.
+
+#ifndef RATE_H
+#define RATE_H
+
+#include <stdio.h>
+#include <assert.h>
+#include "scummsys.h"
+#include "common/engine.h"
+#include "common/util.h"
+
+#include "audiostream.h"
+
+typedef int16 st_sample_t;
+typedef uint32 st_size_t;
+typedef uint32 st_rate_t;
+
+typedef struct {
+	bool used;
+	byte priv[1024];
+} eff_struct;
+typedef eff_struct *eff_t;
+
+/* Minimum and maximum values a sample can hold. */
+#define ST_SAMPLE_MAX 0x7fffL
+#define ST_SAMPLE_MIN (-ST_SAMPLE_MAX - 1L)
+
+#define ST_EOF (-1)
+#define ST_SUCCESS (0)
+
+/* here for linear interp.  might be useful for other things */
+static st_rate_t st_gcd(st_rate_t a, st_rate_t b)
+{
+	if (b == 0)
+		return a;
+	else
+		return st_gcd(b, a % b);
+}
+
+static inline void clampedAdd(int16& a, int b) {
+	int val = a + b;
+
+	if (val > ST_SAMPLE_MAX)
+		a = ST_SAMPLE_MAX;
+	else if (val < ST_SAMPLE_MIN)
+		a = ST_SAMPLE_MIN;
+	else
+		a = val;
+}
+
+// Q&D hack to get this SOX stuff to work
+#define st_report warning
+#define st_warn warning
+#define st_fail error
+
+
+// Resample (high quality)
+int st_resample_getopts(eff_t effp, int n, char **argv);
+int st_resample_start(eff_t effp, st_rate_t inrate, st_rate_t outrate);
+int st_resample_flow(eff_t effp, InputStream &input, st_sample_t *obuf, st_size_t *osamp);
+int st_resample_drain(eff_t effp, st_sample_t *obuf, st_size_t *osamp);
+int st_resample_stop(eff_t effp);
+
+// Rate (linear filter, low quality)
+int st_rate_getopts(eff_t effp, int n, char **argv);
+int st_rate_start(eff_t effp, st_rate_t inrate, st_rate_t outrate);
+int st_rate_flow(eff_t effp, InputStream &input, st_sample_t *obuf, st_size_t *osamp);
+int st_rate_stop(eff_t effp);
+
+#endif
diff --git a/sound/resample.cpp b/sound/resample.cpp
new file mode 100644
index 0000000000..8dc63af6fb
--- /dev/null
+++ b/sound/resample.cpp
@@ -0,0 +1,713 @@
+/*
+* July 5, 1991
+* Copyright 1991 Lance Norskog And Sundry Contributors
+* This source code is freely redistributable and may be used for
+* any purpose.	 This copyright notice must be maintained. 
+* Lance Norskog And Sundry Contributors are not responsible for 
+* the consequences of using this software.
+*/
+
+/*
+ * Sound Tools rate change effect file.
+ * Spiffy rate changer using Smith & Wesson Bandwidth-Limited Interpolation.
+ * The algorithm is described in "Bandlimited Interpolation -
+ * Introduction and Algorithm" by Julian O. Smith III.
+ * Available on ccrma-ftp.stanford.edu as
+ * pub/BandlimitedInterpolation.eps.Z or similar.
+ *
+ * The latest stand alone version of this algorithm can be found
+ * at ftp://ccrma-ftp.stanford.edu/pub/NeXT/
+ * under the name of resample-version.number.tar.Z
+ *
+ * NOTE: There is a newer version of the resample routine then what
+ * this file was originally based on.  Those adventurous might be
+ * interested in reviewing its improvesments and porting it to this
+ * version.
+ */
+
+/* Fixed bug: roll off frequency was wrong, too high by 2 when upsampling,
+ * too low by 2 when downsampling.
+ * Andreas Wilde, 12. Feb. 1999, andreas@eakaw2.et.tu-dresden.de
+*/
+
+/*
+ * October 29, 1999
+ * Various changes, bugfixes(?), increased precision, by Stan Brooks.
+ *
+ * This source code is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ */ 
+/*
+ * SJB: [11/25/99]
+ * TODO: another idea for improvement...
+ * note that upsampling usually doesn't require interpolation,
+ * therefore is faster and more accurate than downsampling.
+ * Downsampling by an integer factor is also simple, since
+ * it just involves decimation if the input is already 
+ * lowpass-filtered to the output Nyquist freqency.
+ * Get the idea? :)
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "rate.h"
+
+/* resample includes */
+#include "resample.h"
+
+/* this Float MUST match that in filter.c */
+#define Float double/*float*/
+
+/* largest factor for which exact-coefficients upsampling will be used */
+#define NQMAX 511
+
+#define BUFFSIZE 8192 /*16384*/	 /* Total I/O buffer size */
+
+/* Private data for Lerp via LCM file */
+typedef struct resamplestuff {
+	double Factor;     /* Factor = Fout/Fin sample rates */
+	double rolloff;    /* roll-off frequency */
+	double beta;	      /* passband/stopband tuning magic */
+	int quadr;	      /* non-zero to use qprodUD quadratic interpolation */
+	long Nmult;
+	long Nwing;
+	long Nq;
+	Float *Imp;	      /* impulse [Nwing+1] Filter coefficients */
+
+	double Time;	      /* Current time/pos in input sample */
+	long dhb;
+
+	long a, b;	      /* gcd-reduced input,output rates	  */
+	long t;	      /* Current time/pos for exact-coeff's method */
+
+	long Xh;	      /* number of past/future samples needed by filter	 */
+	long Xoff;	      /* Xh plus some room for creep  */
+	long Xread;	      /* X[Xread] is start-position to enter new samples */
+	long Xp;	      /* X[Xp] is position to start filter application	 */
+	long Xsize, Ysize;  /* size (Floats) of X[],Y[]	  */
+	long Yposition;		/* FIXME: offset into Y buffer */
+	Float *X, *Y;      /* I/O buffers */
+} *resample_t;
+
+static void LpFilter(double c[],
+                     long N,
+                     double frq,
+                     double Beta,
+                     long Num);
+
+/* makeFilter is used by filter.c */
+int makeFilter(Float Imp[],
+               long Nwing,
+               double Froll,
+               double Beta,
+               long Num,
+               int Normalize);
+
+static long SrcUD(resample_t r, long Nx);
+static long SrcEX(resample_t r, long Nx);
+
+
+/*
+ * Process options
+ */
+int st_resample_getopts(eff_t effp, int n, char **argv) {
+	resample_t r = (resample_t) effp->priv;
+
+	/* These defaults are conservative with respect to aliasing. */
+	r->rolloff = 0.80;
+	r->beta = 16; /* anything <=2 means Nutall window */
+	r->quadr = 0;
+	r->Nmult = 45;
+
+	/* This used to fail, but with sox-12.15 it works. AW */
+	if ((n >= 1)) {
+		if (!strcmp(argv[0], "-qs")) {
+			r->quadr = 1;
+			n--;
+			argv++;
+		} else if (!strcmp(argv[0], "-q")) {
+			r->rolloff = 0.875;
+			r->quadr = 1;
+			r->Nmult = 75;
+			n--;
+			argv++;
+		} else if (!strcmp(argv[0], "-ql")) {
+			r->rolloff = 0.94;
+			r->quadr = 1;
+			r->Nmult = 149;
+			n--;
+			argv++;
+		}
+	}
+
+	if ((n >= 1) && (sscanf(argv[0], "%lf", &r->rolloff) != 1)) {
+		st_fail("Usage: resample [ rolloff [ beta ] ]");
+		return (ST_EOF);
+	} else if ((r->rolloff <= 0.01) || (r->rolloff >= 1.0)) {
+		st_fail("resample: rolloff factor (%f) no good, should be 0.01<x<1.0", r->rolloff);
+		return (ST_EOF);
+	}
+
+	if ((n >= 2) && !sscanf(argv[1], "%lf", &r->beta)) {
+		st_fail("Usage: resample [ rolloff [ beta ] ]");
+		return (ST_EOF);
+	} else if (r->beta <= 2.0) {
+		r->beta = 0;
+		st_report("resample opts: Nuttall window, cutoff %f\n", r->rolloff);
+	} else {
+		st_report("resample opts: Kaiser window, cutoff %f, beta %f\n", r->rolloff, r->beta);
+	}
+	return (ST_SUCCESS);
+}
+
+/*
+ * Prepare processing.
+ */
+int st_resample_start(eff_t effp, st_rate_t inrate, st_rate_t outrate) {
+	resample_t r = (resample_t) effp->priv;
+	long Xoff, gcdrate;
+	int i;
+
+	if (inrate == outrate) {
+		st_fail("Input and Output rates must be different to use resample effect");
+		return (ST_EOF);
+	}
+
+	r->Factor = (double)outrate / (double)inrate;
+
+	gcdrate = st_gcd(inrate, outrate);
+	r->a = inrate / gcdrate;
+	r->b = outrate / gcdrate;
+
+	if (r->a <= r->b && r->b <= NQMAX) {
+		r->quadr = -1; /* exact coeff's	  */
+		r->Nq = r->b;  /* MAX(r->a,r->b);	*/
+	} else {
+		r->Nq = Nc; /* for now */
+	}
+
+	/* Check for illegal constants */
+# if 0
+	if (Lp >= 16)
+		st_fail("Error: Lp>=16");
+	if (Nb + Nhg + NLpScl >= 32)
+		st_fail("Error: Nb+Nhg+NLpScl>=32");
+	if (Nh + Nb > 32)
+		st_fail("Error: Nh+Nb>32");
+# endif
+
+	/* Nwing: # of filter coeffs in right wing */
+	r->Nwing = r->Nq * (r->Nmult / 2 + 1) + 1;
+
+	r->Imp = (Float *)malloc(sizeof(Float) * (r->Nwing + 2)) + 1;
+	/* need Imp[-1] and Imp[Nwing] for quadratic interpolation */
+	/* returns error # <=0, or adjusted wing-len > 0 */
+	i = makeFilter(r->Imp, r->Nwing, r->rolloff, r->beta, r->Nq, 1);
+	if (i <= 0) {
+		st_fail("resample: Unable to make filter\n");
+		return (ST_EOF);
+	}
+
+	st_report("Nmult: %ld, Nwing: %ld, Nq: %ld\n",r->Nmult,r->Nwing,r->Nq);	// FIXME
+
+	if (r->quadr < 0) { /* exact coeff's method */
+		r->Xh = r->Nwing / r->b;
+		st_report("resample: rate ratio %ld:%ld, coeff interpolation not needed\n", r->a, r->b);
+	} else {
+		r->dhb = Np;	/* Fixed-point Filter sampling-time-increment */
+		if (r->Factor < 1.0)
+			r->dhb = (long)(r->Factor * Np + 0.5);
+		r->Xh = (r->Nwing << La) / r->dhb;
+		/* (Xh * dhb)>>La is max index into Imp[] */
+	}
+
+	/* reach of LP filter wings + some creeping room */
+	Xoff = r->Xh + 10;
+	r->Xoff = Xoff;
+
+	/* Current "now"-sample pointer for input to filter */
+	r->Xp = Xoff;
+	/* Position in input array to read into */
+	r->Xread = Xoff;
+	/* Current-time pointer for converter */
+	r->Time = Xoff;
+	if (r->quadr < 0) { /* exact coeff's method */
+		r->t = Xoff * r->Nq;
+	}
+	i = BUFFSIZE - 2 * Xoff;
+	if (i < r->Factor + 1.0 / r->Factor)	/* Check input buffer size */
+	{
+		st_fail("Factor is too small or large for BUFFSIZE");
+		return (ST_EOF);
+	}
+
+	r->Xsize = (long)(2 * Xoff + i / (1.0 + r->Factor));
+	r->Ysize = BUFFSIZE - r->Xsize;
+	st_report("Xsize %ld, Ysize %ld, Xoff %ld",r->Xsize,r->Ysize,r->Xoff);	// FIXME
+
+	r->X = (Float *) malloc(sizeof(Float) * (BUFFSIZE));
+	r->Y = r->X + r->Xsize;
+	r->Yposition = 0;
+
+	/* Need Xoff zeros at beginning of sample */
+	for (i = 0; i < Xoff; i++)
+		r->X[i] = 0;
+	return (ST_SUCCESS);
+}
+
+/*
+ * Processed signed long samples from ibuf to obuf.
+ * Return number of samples processed.
+ */
+int st_resample_flow(eff_t effp, InputStream &input, st_sample_t *obuf, st_size_t *osamp) {
+	resample_t r = (resample_t) effp->priv;
+	long i, k, last;
+	long Nout;		// The number of bytes we effectively output
+	long Nx;		// The number of bytes we will read from input
+	long Nproc;		// The number of bytes we process to generate Nout output bytes
+
+#if 1	// FIXME: Hack to generate stereo output
+	*osamp >>= 1;
+#endif
+
+	/* constrain amount we actually process */
+	fprintf(stderr,"Xp %d, Xread %d\n",r->Xp, r->Xread);
+
+	// Initially assume we process the full X buffer starting at the filter
+	// start position.
+	Nproc = r->Xsize - r->Xp;
+printf("FOO(1) Nproc %ld\n", Nproc);
+
+	// Nproc is bounded indirectly by the size of output buffer, and also by
+	// the remaining size of the Y buffer (whichever is smaller).
+	i = MIN(r->Ysize - r->Yposition, (long)*osamp);
+	if (Nproc * r->Factor >= i)
+		Nproc = (long)(i / r->Factor);
+
+printf("FOO(2) Nproc %ld\n", Nproc);
+
+	// Now that we know how many bytes we want to process, we determine
+	// how many bytes to read. We already have Xread bytes in our input
+	// buffer, so we need Nproc - r->Xread more bytes.
+	Nx = Nproc - r->Xread + r->Xoff + r->Xp; // FIXME: Fingolfing thinks this is the correct thing, not what's in the next line!
+//	Nx = Nproc - r->Xread; /* space for right-wing future-data */
+	if (Nx <= 0) {
+		st_fail("resample: Can not handle this sample rate change. Nx not positive: %d", Nx);
+		return (ST_EOF);
+	}
+	
+	// Nx is the number of bytes we'd like to read, but of course that is limited
+	// by the number of bytes actually available...
+	if (Nx > (long)input.size())
+		Nx = (long)input.size();
+	fprintf(stderr,"Nx %d\n",Nx);
+
+	// Read in Nx bytes
+	for (i = r->Xread; i < Nx + r->Xread ; i++)
+		r->X[i] = (Float)input.read();
+
+	last = Nx + r->Xread;	// 'last' is the idx after the last valid byte in X (i.e. number of bytes are in buffer X right now)
+	
+	// Finally compute the effective number of bytes to process
+	Nproc = last - r->Xoff - r->Xp;
+printf("FOO(3) Nproc %ld\n", Nproc);
+
+	if (Nproc <= 0) {
+		/* fill in starting here next time */
+		r->Xread = last;
+		/* leave *isamp alone, we consumed it */
+		*osamp = 0;
+		return (ST_SUCCESS);
+	}
+	if (r->quadr < 0) { /* exact coeff's method */
+		long creep;
+		Nout = SrcEX(r, Nproc);
+		fprintf(stderr,"Nproc %d --> %d\n",Nproc,Nout);
+		/* Move converter Nproc samples back in time */
+		r->t -= Nproc * r->b;
+		/* Advance by number of samples processed */
+		r->Xp += Nproc;
+		/* Calc time accumulation in Time */
+		creep = r->t / r->b - r->Xoff;
+		if (creep) {
+			r->t -= creep * r->b;	 /* Remove time accumulation   */
+			r->Xp += creep;	 /* and add it to read pointer */
+			fprintf(stderr,"Nproc %ld, creep %ld\n",Nproc,creep);
+		}
+	} else { /* approx coeff's method */
+		long creep;
+		Nout = SrcUD(r, Nproc);
+		fprintf(stderr,"Nproc %d --> %d\n",Nproc,Nout);
+		/* Move converter Nproc samples back in time */
+		r->Time -= Nproc;
+		/* Advance by number of samples processed */
+		r->Xp += Nproc;
+		/* Calc time accumulation in Time */
+		creep = (long)(r->Time - r->Xoff);
+		if (creep) {
+			r->Time -= creep;   /* Remove time accumulation   */
+			r->Xp += creep;     /* and add it to read pointer */
+			fprintf(stderr,"Nproc %ld, creep %ld\n",Nproc,creep);
+		}
+	}
+
+	/* Copy back portion of input signal that must be re-used */
+	k = r->Xp - r->Xoff;
+	fprintf(stderr,"k %d, last %d\n",k,last);
+	for (i = 0; i < last - k; i++)
+		r->X[i] = r->X[i + k];
+
+	/* Pos in input buff to read new data into */
+	r->Xread = i;
+	r->Xp = r->Xoff;
+
+printf("osamp = %d, Nout = %ld\n", *osamp, Nout);
+	for (i = 0; i < Nout; i++) {
+		clampedAdd(*obuf++, (int)r->Y[i]);
+#if 1	// FIXME: Hack to generate stereo output
+		clampedAdd(*obuf++, (int)r->Y[i]);
+#endif
+	}
+
+	// At this point, we used *osamp bytes out of Nout available bytes in the Y buffer.
+	// If there are any bytes remaining, shift them to the start of the buffer,
+	// and update Yposition
+
+
+	*osamp = Nout;
+
+	return (ST_SUCCESS);
+}
+
+/*
+ * Process tail of input samples.
+ */
+int st_resample_drain(eff_t effp, st_sample_t *obuf, st_size_t *osamp) {
+	resample_t r = (resample_t) effp->priv;
+	long isamp_res, osamp_res;
+	st_sample_t *Obuf;
+	int rc;
+
+	/*fprintf(stderr,"Xoff %d, Xt %d  <--- DRAIN\n",r->Xoff, r->Xt);*/
+
+	/* stuff end with Xoff zeros */
+	isamp_res = r->Xoff;
+	osamp_res = *osamp;
+	Obuf = obuf;
+	while (isamp_res > 0 && osamp_res > 0) {
+		st_sample_t Osamp;
+		Osamp = osamp_res;
+		ZeroInputStream zero(isamp_res);
+		rc = st_resample_flow(effp, zero, Obuf, (st_size_t *) & Osamp);
+		if (rc)
+			return rc;
+		/*fprintf(stderr,"DRAIN isamp,osamp	(%d,%d) -> (%d,%d)\n",
+		    isamp_res,osamp_res,Isamp,Osamp);*/
+		Obuf += Osamp;
+		osamp_res -= Osamp;
+		isamp_res = zero.size();
+	}
+	*osamp -= osamp_res;
+	fprintf(stderr,"DRAIN osamp %d\n", *osamp);
+	if (isamp_res)
+		st_warn("drain overran obuf by %d\n", isamp_res);
+	fflush(stderr);
+	return (ST_SUCCESS);
+}
+
+/*
+ * Do anything required when you stop reading samples.	
+ * Don't close input file! 
+ */
+int st_resample_stop(eff_t effp) {
+	resample_t r = (resample_t) effp->priv;
+
+	free(r->Imp - 1);
+	free(r->X);
+	/* free(r->Y); Y is in same block starting at X */
+	return (ST_SUCCESS);
+}
+
+/* over 90% of CPU time spent in this iprodUD() function */
+/* quadratic interpolation */
+static double qprodUD(const Float Imp[], const Float *Xp, long Inc, double T0,
+                      long dhb, long ct) {
+	const double f = 1.0 / (1 << La);
+	double v;
+	long Ho;
+
+	Ho = (long)(T0 * dhb);
+	Ho += (ct - 1) * dhb; /* so Float sum starts with smallest coef's */
+	Xp += (ct - 1) * Inc;
+	v = 0;
+	do {
+		Float coef;
+		long Hoh;
+		Hoh = Ho >> La;
+		coef = Imp[Hoh];
+		{
+			Float dm, dp, t;
+			dm = coef - Imp[Hoh - 1];
+			dp = Imp[Hoh + 1] - coef;
+			t = (Ho & Amask) * f;
+			coef += ((dp - dm) * t + (dp + dm)) * t * 0.5;
+		}
+		/* filter coef, lower La bits by quadratic interpolation */
+		v += coef * *Xp;   /* sum coeff * input sample */
+		Xp -= Inc;	   /* Input signal step. NO CHECK ON ARRAY BOUNDS */
+		Ho -= dhb;	   /* IR step */
+	} while (--ct);
+	return v;
+}
+
+/* linear interpolation */
+static double iprodUD(const Float Imp[], const Float *Xp, long Inc,
+                      double T0, long dhb, long ct) {
+	const double f = 1.0 / (1 << La);
+	double v;
+	long Ho;
+
+	Ho = (long)(T0 * dhb);
+	Ho += (ct - 1) * dhb; /* so Float sum starts with smallest coef's */
+	Xp += (ct - 1) * Inc;
+	v = 0;
+	do {
+		Float coef;
+		long Hoh;
+		Hoh = Ho >> La;
+		/* if (Hoh >= End) break; */
+		coef = Imp[Hoh] + (Imp[Hoh + 1] - Imp[Hoh]) * (Ho & Amask) * f;
+		/* filter coef, lower La bits by linear interpolation */
+		v += coef * *Xp;   /* sum coeff * input sample */
+		Xp -= Inc;	   /* Input signal step. NO CHECK ON ARRAY BOUNDS */
+		Ho -= dhb;	   /* IR step */
+	} while (--ct);
+	return v;
+}
+
+/* From resample:filters.c */
+/* Sampling rate conversion subroutine */
+
+static long SrcUD(resample_t r, long Nx) {
+	Float *Ystart, *Y;
+	double Factor;
+	double dt;		       /* Step through input signal */
+	double time;
+	double (*prodUD)(const Float Imp[], const Float *Xp, long Inc, double T0, long dhb, long ct);
+	int n;
+
+	prodUD = (r->quadr) ? qprodUD : iprodUD; /* quadratic or linear interp */
+	Factor = r->Factor;
+	time = r->Time;
+	dt = 1.0 / Factor;	   /* Output sampling period */
+	fprintf(stderr,"Factor %f, dt %f, ",Factor,dt);
+	fprintf(stderr,"Time %f, ",r->Time);
+	/* (Xh * dhb)>>La is max index into Imp[] */
+	/*fprintf(stderr,"ct=%d\n",ct);*/
+	fprintf(stderr,"ct=%.2f %d\n",(double)r->Nwing*Na/r->dhb, r->Xh);
+	fprintf(stderr,"ct=%ld, T=%.6f, dhb=%6f, dt=%.6f\n", r->Xh, time-floor(time),(double)r->dhb/Na,dt);
+	Ystart = Y = r->Y + r->Yposition;
+	n = (int)ceil((double)Nx / dt);
+	while (n--) {
+		Float *Xp;
+		double v;
+		double T;
+		T = time - floor(time);	   /* fractional part of Time */
+		Xp = r->X + (long)time;	   /* Ptr to current input sample */
+
+		/* Past  inner product: */
+		v = (*prodUD)(r->Imp, Xp, -1, T, r->dhb, r->Xh); /* needs Np*Nmult in 31 bits */
+		/* Future inner product: */
+		v += (*prodUD)(r->Imp, Xp + 1, 1, (1.0 - T), r->dhb, r->Xh); /* prefer even total */
+
+		if (Factor < 1)
+			v *= Factor;
+		*Y++ = v;		     /* Deposit output */
+		time += dt;	     /* Move to next sample by time increment */
+	}
+	r->Time = time;
+	fprintf(stderr,"Time %f\n",r->Time);
+	return (Y - Ystart);	       /* Return the number of output samples */
+}
+
+/* exact coeff's */
+static double prodEX(const Float Imp[], const Float *Xp,
+                     long Inc, long T0, long dhb, long ct) {
+	double v;
+	const Float *Cp;
+
+	Cp = Imp + (ct - 1) * dhb + T0; /* so Float sum starts with smallest coef's */
+	Xp += (ct - 1) * Inc;
+	v = 0;
+	do {
+		v += *Cp * *Xp;   /* sum coeff * input sample */
+		Cp -= dhb;	   /* IR step */
+		Xp -= Inc;	   /* Input signal step. */
+	} while (--ct);
+	return v;
+}
+
+static long SrcEX(resample_t r, long Nx) {
+	Float *Ystart, *Y;
+	double Factor;
+	long a, b;
+	long time;
+	int n;
+
+	Factor = r->Factor;
+	time = r->t;
+	a = r->a;
+	b = r->b;
+	Ystart = Y = r->Y + r->Yposition;
+	n = (Nx * b + (a - 1)) / a;
+	while (n--) {
+		Float *Xp;
+		double v;
+		long T;
+		T = time % b;		   /* fractional part of Time */
+		Xp = r->X + (time / b);	   /* Ptr to current input sample */
+
+		/* Past	 inner product: */
+		v = prodEX(r->Imp, Xp, -1, T, b, r->Xh);
+		/* Future inner product: */
+		v += prodEX(r->Imp, Xp + 1, 1, b - T, b, r->Xh);
+
+		if (Factor < 1)
+			v *= Factor;
+		*Y++ = v;	      /* Deposit output */
+		time += a;	      /* Move to next sample by time increment */
+	}
+	r->t = time;
+	return (Y - Ystart);	       /* Return the number of output samples */
+}
+
+int makeFilter(Float Imp[], long Nwing, double Froll, double Beta,
+               long Num, int Normalize) {
+	double *ImpR;
+	long Mwing, i;
+
+	if (Nwing > MAXNWING)		      /* Check for valid parameters */
+		return ( -1);
+	if ((Froll <= 0) || (Froll > 1))
+		return ( -2);
+
+	/* it does help accuracy a bit to have the window stop at
+	 * a zero-crossing of the sinc function */
+	Mwing = (long)(floor((double)Nwing / (Num / Froll)) * (Num / Froll) + 0.5);
+	if (Mwing == 0)
+		return ( -4);
+
+	ImpR = (double *) malloc(sizeof(double) * Mwing);
+
+	/* Design a Nuttall or Kaiser windowed Sinc low-pass filter */
+	LpFilter(ImpR, Mwing, Froll, Beta, Num);
+
+	if (Normalize) { /* 'correct' the DC gain of the lowpass filter */
+		long Dh;
+		double DCgain;
+		DCgain = 0;
+		Dh = Num;			 /* Filter sampling period for factors>=1 */
+		for (i = Dh; i < Mwing; i += Dh)
+			DCgain += ImpR[i];
+		DCgain = 2 * DCgain + ImpR[0];    /* DC gain of real coefficients */
+		st_report("DCgain err=%.12f",DCgain-1.0);	// FIXME
+
+		DCgain = 1.0 / DCgain;
+		for (i = 0; i < Mwing; i++)
+			Imp[i] = ImpR[i] * DCgain;
+
+	} else {
+		for (i = 0; i < Mwing; i++)
+			Imp[i] = ImpR[i];
+	}
+	free(ImpR);
+	for (i = Mwing; i <= Nwing; i++)
+		Imp[i] = 0;
+	/* Imp[Mwing] and Imp[-1] needed for quadratic interpolation */
+	Imp[ -1] = Imp[1];
+
+	return (Mwing);
+}
+
+/* LpFilter()
+ *
+ * reference: "Digital Filters, 2nd edition"
+ *	      R.W. Hamming, pp. 178-179
+ *
+ * Izero() computes the 0th order modified bessel function of the first kind.
+ *    (Needed to compute Kaiser window).
+ *
+ * LpFilter() computes the coeffs of a Kaiser-windowed low pass filter with
+ *    the following characteristics:
+ *
+ *	 c[]  = array in which to store computed coeffs
+ *	 frq  = roll-off frequency of filter
+ *	 N    = Half the window length in number of coeffs
+ *	 Beta = parameter of Kaiser window
+ *	 Num  = number of coeffs before 1/frq
+ *
+ * Beta trades the rejection of the lowpass filter against the transition
+ *    width from passband to stopband.	Larger Beta means a slower
+ *    transition and greater stopband rejection.  See Rabiner and Gold
+ *    (Theory and Application of DSP) under Kaiser windows for more about
+ *    Beta.  The following table from Rabiner and Gold gives some feel
+ *    for the effect of Beta:
+ *
+ * All ripples in dB, width of transition band = D*N where N = window length
+ *
+ *		 BETA	 D	 PB RIP	  SB RIP
+ *		 2.120	 1.50  +-0.27	   -30
+ *		 3.384	 2.23	 0.0864	   -40
+ *		 4.538	 2.93	 0.0274	   -50
+ *		 5.658	 3.62	 0.00868   -60
+ *		 6.764	 4.32	 0.00275   -70
+ *		 7.865	 5.0	 0.000868  -80
+ *		 8.960	 5.7	 0.000275  -90
+ *		 10.056	 6.4	 0.000087  -100
+ */
+
+
+#define IzeroEPSILON 1E-21		 /* Max error acceptable in Izero */
+
+static double Izero(double x) {
+	double sum, u, halfx, temp;
+	long n;
+
+	sum = u = n = 1;
+	halfx = x / 2.0;
+	do {
+		temp = halfx / (double)n;
+		n += 1;
+		temp *= temp;
+		u *= temp;
+		sum += u;
+	} while (u >= IzeroEPSILON*sum);
+	return (sum);
+}
+
+static void LpFilter(double *c, long N, double frq, double Beta, long Num) {
+	long i;
+
+	/* Calculate filter coeffs: */
+	c[0] = frq;
+	for (i = 1; i < N; i++) {
+		double x = M_PI * (double)i / (double)(Num);
+		c[i] = sin(x * frq) / x;
+	}
+
+	if (Beta > 2) { /* Apply Kaiser window to filter coeffs: */
+		double IBeta = 1.0 / Izero(Beta);
+		for (i = 1; i < N; i++) {
+			double x = (double)i / (double)(N);
+			c[i] *= Izero(Beta * sqrt(1.0 - x * x)) * IBeta;
+		}
+	} else { /* Apply Nuttall window: */
+		for (i = 0; i < N; i++) {
+			double x = M_PI * i / N;
+			c[i] *= 0.36335819 + 0.4891775 * cos(x) + 0.1365995 * cos(2 * x) + 0.0106411 * cos(3 * x);
+		}
+	}
+}
diff --git a/sound/resample.h b/sound/resample.h
new file mode 100644
index 0000000000..f66ad75a5d
--- /dev/null
+++ b/sound/resample.h
@@ -0,0 +1,61 @@
+/*
+ * FILE: resample.h
+ *   BY: Julius Smith (at CCRMA, Stanford U)
+ * C BY: translated from SAIL to C by Christopher Lee Fraley
+ *          (cf0v@andrew.cmu.edu)
+ * DATE: 7-JUN-88
+ * VERS: 2.0  (17-JUN-88, 3:00pm)
+ */
+
+/*
+ * October 29, 1999
+ * Various changes, bugfixes(?), increased precision, by Stan Brooks.
+ *
+ * This source code is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ */
+
+/* Conversion constants */
+#define Lc        7
+#define Nc       (1<<Lc)
+#define La        16
+#define Na       (1<<La)
+#define Lp       (Lc+La)
+#define Np       (1<<Lp)
+#define Amask    (Na-1)
+#define Pmask    (Np-1)
+
+#define MAXNWING  (80<<Lc)
+/* Description of constants:
+ *
+ * Nc - is the number of look-up values available for the lowpass filter
+ *    between the beginning of its impulse response and the "cutoff time"
+ *    of the filter.  The cutoff time is defined as the reciprocal of the
+ *    lowpass-filter cut off frequence in Hz.  For example, if the
+ *    lowpass filter were a sinc function, Nc would be the index of the
+ *    impulse-response lookup-table corresponding to the first zero-
+ *    crossing of the sinc function.  (The inverse first zero-crossing
+ *    time of a sinc function equals its nominal cutoff frequency in Hz.)
+ *    Nc must be a power of 2 due to the details of the current
+ *    implementation. The default value of 128 is sufficiently high that
+ *    using linear interpolation to fill in between the table entries
+ *    gives approximately 16-bit precision, and quadratic interpolation
+ *    gives about 23-bit (float) precision in filter coefficients.
+ *
+ * Lc - is log base 2 of Nc.
+ *
+ * La - is the number of bits devoted to linear interpolation of the
+ *    filter coefficients.
+ *
+ * Lp - is La + Lc, the number of bits to the right of the binary point
+ *    in the integer "time" variable. To the left of the point, it indexes
+ *    the input array (X), and to the right, it is interpreted as a number
+ *    between 0 and 1 sample of the input X.  The default value of 23 is
+ *    about right.  There is a constraint that the filter window must be
+ *    "addressable" in a int32_t, more precisely, if Nmult is the number
+ *    of sinc zero-crossings in the right wing of the filter window, then
+ *    (Nwing<<Lp) must be expressible in 31 bits.
+ *
+ */