From b9d380bba43b25e1654b320a42b05499c060656c Mon Sep 17 00:00:00 2001 From: Max Horn Date: Thu, 24 Jul 2003 17:46:38 +0000 Subject: [PATCH] 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 --- sound/audiostream.h | 91 ++++++ sound/rate.cpp | 164 ++++++++++ sound/rate.h | 72 +++++ sound/resample.cpp | 713 ++++++++++++++++++++++++++++++++++++++++++++ sound/resample.h | 61 ++++ 5 files changed, 1101 insertions(+) create mode 100644 sound/audiostream.h create mode 100644 sound/rate.cpp create mode 100644 sound/rate.h create mode 100644 sound/resample.cpp create mode 100644 sound/resample.h diff --git a/sound/audiostream.h b/sound/audiostream.h new file mode 100644 index 00000000000..c3d43d6b2a3 --- /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 +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 +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 +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 +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 00000000000..bde477647af --- /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 +/* + * 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 00000000000..aa8edea1b00 --- /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 +#include +#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 00000000000..8dc63af6fb5 --- /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 +#include +#include +#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.01rolloff); + 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 00000000000..f66ad75a5d9 --- /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<