thirdparty/scummvm
2
0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

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

Browse source 
svn-id: r9176
This commit is contained in:
Max Horn 2003-07-24 17:46:38 +00:00
parent becd70d243
commit b9d380bba4
5 changed files with 1101 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

91
sound/audiostream.h Normal file
View file

@ -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

164
sound/rate.cpp Normal file
View file

@ -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);
}

72
sound/rate.h Normal file
View file

@ -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

713
sound/resample.cpp Normal file
View file

@ -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);
}
}
}

61
sound/resample.h Normal file
View file

@ -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.
*
*/