ref: 82ccfd17ad7ad52e64651a9d63092a3e8f282567
dir: /src/speech/rateconv.c/
/* THIS IS A MODIFIED VERSION. It was modified by David
Huggins-Daines <dhd@cepstral.com> in December 2001 for use in the
Flite speech synthesis systems. */
/*
*
* RATECONV.C
*
* Convert sampling rate stdin to stdout
*
* Copyright (c) 1992, 1995 by Markus Mummert
*
* Redistribution and use of this software, modifcation and inclusion
* into other forms of software are permitted provided that the following
* conditions are met:
*
* 1. Redistributions of this software must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. If this software is redistributed in a modified condition
* it must reveal clearly that it has been modified.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS''
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
* PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR
* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
* PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
* PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
* OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
* USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*
* history: 2.9.92 begin coding
* 5.9.92 fully operational
* 14.2.95 provide BIG_ENDIAN, SWAPPED_BYTES_DEFAULT
* switches, Copyright note and References
* 25.11.95 changed XXX_ENDIAN to I_AM_XXX_ENDIAN
* default gain set to 0.8
* 3.12.95 stereo implementation
* SWAPPED_BYTES_DEFAULT now HBYTE1ST_DEFAULT
* changed [L/2] to (L-1)/2 for exact symmetry
*
*
* IMPLEMENTATION NOTES
*
* Converting is achieved by interpolating the input samples in
* order to evaluate the represented continuous input slope at
* sample instances of the new rate (resampling). It is implemented
* as a polyphase FIR-filtering process (see reference). The rate
* conversion factor is determined by a rational factor. Its
* nominator and denominator are integers of almost arbitrary
* value, limited only by coefficient memory size.
*
* General rate conversion formula:
*
* out(n*Tout) = SUM in(m*Tin) * g((n*d/u-m)*Tin) * Tin
* over all m
*
* FIR-based rate conversion formula for polyphase processing:
*
* L-1
* out(n*Tout) = SUM in(A(i,n)*Tin) * g(B(i,n)*Tin) * Tin
* i=0
*
* A(i,n) = i - (L-1)/2 + [n*d/u]
* = i - (L-1)/2 + [(n%u)*d/u] + [n/u]*d
* B(i,n) = n*d/u - [n*d/u] + (L-1)/2 - i
* = ((n%u)*d/u)%1 + (L-1)/2 - i
* Tout = Tin * d/u
*
* where:
* n,i running integers
* out(t) output signal sampled at t=n*Tout
* in(t) input signal sampled in intervalls Tin
* u,d up- and downsampling factor, integers
* g(t) interpolating function
* L FIR-length of realized g(t), integer
* / float-division-operator
* % float-modulo-operator
* [] integer-operator
*
* note:
* (L-1)/2 in A(i,n) can be omitted as pure time shift yielding
* a causal design with a delay of ((L-1)/2)*Tin.
* n%u is a cyclic modulo-u counter clocked by out-rate
* [n/u]*d is a d-increment counter, advanced when n%u resets
* B(i,n)*Tin can take on L*u differnt values, at which g(t)
* has to be sampled as a coefficient array
*
* Interpolation function design:
*
* The interpolation function design is based on a sinc-function
* windowed by a gaussian function. The former determines the
* cutoff frequency, the latter limits the necessary FIR-length by
* pushing the outer skirts of the resulting impulse response below
* a certain threshold fast enough. The drawback is a smoothed
* cutoff inducing some aliasing. Due to the symmetry of g(t) the
* group delay of the filtering process is contant (linear phase).
*
* g(t) = 2*fgK*sinc(pi*2*fgK*t) * exp(-pi*(2*fgG*t)**2)
*
* where:
* fgK cutoff frequency of sinc function in f-domain
* fgG key frequency of gaussian window in f-domain
* reflecting the 6.82dB-down point
*
* note:
* Taking fsin=1/Tin as the input sampling frequncy, it turns out
* that in conjunction with L, u and d only the ratios fgK/(fsin/2)
* and fgG/(fsin/2) specify the whole proces. Requiring fsin, fgK
* and fgG as input purposely keeps the notion of absolute
* frequencies.
*
* Numerical design:
*
* Samples are expected to be 16bit-signed integers, alternating
* left and right channel in case of stereo mode- The byte order
* per sample is selectable. FIR-filtering is implemented using
* 32bit-signed integer arithmetic. Coefficients are scaled to
* find the output sample in the high word of accumulated FIR-sum.
*
* Interpolation can lead to sample magnitudes exceeding the
* input maximum. Worst case is a full scale step function on the
* input. In this case the sinc-function exhibits an overshoot of
* 2*9=18percent (Gibb's phaenomenon). In any case sample overflow
* can be avoided by a gain of 0.8.
*
* If u=d=1 and if the input stream contains only a single sample,
* the whole length of the FIR-filter will be written to the output.
* In general the resulting output signal will be (L-1)*fsout/fsin
* samples longer than the input signal. The effect is that a
* finite input sequence is viewed as padded with zeros before the
* `beginning' and after the `end'.
*
* The output lags ((L-1)/2)*Tin behind to implement g(t) as a
* causal system corresponding to a causal relationship of the
* discrete-time sequences in(m/fsin) and out(n/fsout) with
* resepect to a sequence time origin at t=n*Tin=m*Tout=0.
*
*
* REFERENCES
*
* Crochiere, R. E., Rabiner, L. R.: "Multirate Digital Signal
* Processing", Prentice-Hall, Englewood Cliffs, New Jersey, 1983
*
* Zwicker, E., Fastl, H.: "Psychoacoustics - Facts and Models",
* Springer-Verlag, Berlin, Heidelberg, New-York, Tokyo, 1990 */
#include "cst_string.h"
#include "cst_math.h"
#include "cst_alloc.h"
#include "cst_error.h"
#include "cst_wave.h"
/*
* adaptable defines and globals
*/
#define DPRINTF(l, x)
#define FIXSHIFT 16
#define FIXMUL (1<<FIXSHIFT)
#ifndef M_PI
#define M_PI 3.1415926535
#endif
#define sqr(a) ((a)*(a))
/*
* FIR-routines, mono and stereo
* this is where we need all the MIPS
*/
static void
fir_mono(int *inp, int *coep, int firlen, int *outp)
{
register int akku = 0, *endp;
int n1 = (firlen / 8) * 8, n0 = firlen % 8;
endp = coep + n1;
while (coep != endp) {
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
akku += *inp++ * *coep++;
}
endp = coep + n0;
while (coep != endp) {
akku += *inp++ * *coep++;
}
*outp = akku;
}
static void
fir_stereo(int *inp, int *coep,
int firlen, int *out1p, int *out2p)
{
register int akku1 = 0, akku2 = 0, *endp;
int n1 = (firlen / 8) * 8, n0 = firlen % 8;
endp = coep + n1;
while (coep != endp) {
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
}
endp = coep + n0;
while (coep != endp) {
akku1 += *inp++ * *coep;
akku2 += *inp++ * *coep++;
}
*out1p = akku1;
*out2p = akku2;
}
/*
* filtering from input buffer to output buffer;
* returns number of processed samples in output buffer:
* if it is not equal to output buffer size,
* the input buffer is expected to be refilled upon entry, so that
* the last firlen numbers of the old input buffer are
* the first firlen numbers of the new input buffer;
* if it is equal to output buffer size, the output buffer
* is full and is expected to be stowed away;
*
*/
static int
filtering_on_buffers(cst_rateconv *filt)
{
int insize;
DPRINTF(0, ("filtering_on_buffers(%d)\n", filt->incount));
insize = filt->incount + filt->lag;
if (filt->channels == 1) {
while (1) {
filt->inoffset = (filt->cycctr * filt->down)/filt->up;
if ((filt->inbaseidx + filt->inoffset + filt->len) > insize) {
filt->inbaseidx -= insize - filt->len + 1;
memcpy(filt->sin, filt->sin + insize - filt->lag,
filt->lag * sizeof(int));
/* Prevent people from re-filtering the same stuff. */
filt->incount = 0;
return 0;
}
fir_mono(filt->sin + filt->inoffset + filt->inbaseidx,
filt->coep + filt->cycctr * filt->len,
filt->len, filt->sout + filt->outidx);
DPRINTF(1, ("in(%d + %d) = %d cycctr %d out(%d) = %d\n",
filt->inoffset, filt->inbaseidx,
filt->sin[filt->inoffset + filt->inbaseidx],
filt->cycctr, filt->outidx,
filt->sout[filt->outidx] >> FIXSHIFT));
++filt->outidx;
++filt->cycctr;
if (!(filt->cycctr %= filt->up))
filt->inbaseidx += filt->down;
if (!(filt->outidx %= filt->outsize))
return filt->outsize;
}
} else if (filt->channels == 2) {
/*
* rule how to convert mono routine to stereo routine:
* firlen, up, down and cycctr relate to samples in general,
* wether mono or stereo; inbaseidx, inoffset and outidx as
* well as insize and outsize still account for mono samples.
*/
while (1) {
filt->inoffset = 2*((filt->cycctr * filt->down)/filt->up);
if ((filt->inbaseidx + filt->inoffset + 2*filt->len) > insize) {
filt->inbaseidx -= insize - 2*filt->len + 2;
return filt->outidx;
}
fir_stereo(filt->sin + filt->inoffset + filt->inbaseidx,
filt->coep + filt->cycctr * filt->len,
filt->len,
filt->sout + filt->outidx,
filt->sout + filt->outidx + 1);
filt->outidx += 2;
++filt->cycctr;
if (!(filt->cycctr %= filt->up))
filt->inbaseidx += 2*filt->down;
if (!(filt->outidx %= filt->outsize))
return filt->outsize;
}
} else {
cst_errmsg("filtering_on_buffers: only 1 or 2 channels supported!\n");
cst_error();
}
return 0;
}
/*
* convert buffer of n samples to ints
*/
static void
sample_to_int(short *buff, int n)
{
short *s, *e;
int *d;
if (n < 1)
return;
s = buff + n;
d = (int*)buff + n;
e = buff;
while (s != e) {
*--d = (int)(*--s);
}
}
/*
* convert buffer of n ints to samples
*/
static void
int_to_sample(short *buff, int n)
{
int *s;
short *e, *d;
if (n < 1)
return;
s = (int *)buff;
d = buff;
e = buff + n;
while (d != e)
*d++ = (*s++ >> FIXSHIFT);
}
/*
* read and convert input sample format to integer
*/
int
cst_rateconv_in(cst_rateconv *filt, const short *inptr, int max)
{
if (max > filt->insize - filt->lag)
max = filt->insize - filt->lag;
if (max > 0) {
memcpy(filt->sin + filt->lag, inptr, max * sizeof(short));
sample_to_int((short *)(filt->sin + filt->lag), max);
}
filt->incount = max;
return max;
}
/*
* do some conversion jobs and write
*/
int
cst_rateconv_out(cst_rateconv *filt, short *outptr, int max)
{
int outsize;
if ((outsize = filtering_on_buffers(filt)) == 0)
return 0;
if (max > outsize)
max = outsize;
int_to_sample((short *)filt->sout, max);
memcpy(outptr, filt->sout, max * sizeof(short));
return max;
}
int
cst_rateconv_leadout(cst_rateconv *filt)
{
memset(filt->sin + filt->lag, 0, filt->lag * sizeof(int));
filt->incount = filt->lag;
return filt->lag;
}
/*
* evaluate sinc(x) = sin(x)/x safely
*/
static double
sinc(double x)
{
return(fabs(x) < 1E-50 ? 1.0 : sin(fmod(x,2*M_PI))/x);
}
/*
* evaluate interpolation function g(t) at t
* integral of g(t) over all times is expected to be one
*/
static double
interpol_func(double t, double fgk, double fgg)
{
return (2*fgk*sinc(M_PI*2*fgk*t)*exp(-M_PI*sqr(2*fgg*t)));
}
/*
* evaluate coefficient from i, q=n%u by sampling interpolation function
* and scale it for integer multiplication used by FIR-filtering
*/
static int
coefficient(int i, int q, cst_rateconv *filt)
{
return (int)(FIXMUL * filt->gain *
interpol_func((fmod(q* filt->down/(double)filt->up,1.0)
+ (filt->len-1)/2.0 - i) / filt->fsin,
filt->fgk, filt->fgg) / filt->fsin);
}
/*
* set up coefficient array
*/
static void
make_coe(cst_rateconv *filt)
{
int i, q;
filt->coep = cst_alloc(int, filt->len * filt->up);
for (i = 0; i < filt->len; i++) {
for (q = 0; q < filt->up; q++) {
filt->coep[q * filt->len + i]
= coefficient(i, q, filt);
}
}
}
cst_rateconv *
new_rateconv(int up, int down, int channels)
{
cst_rateconv *filt;
if (!(channels == 1 || channels == 2)) {
cst_errmsg("new_rateconv: channels must be 1 or 2\n");
cst_error();
}
filt = cst_alloc(cst_rateconv, 1);
filt->fsin = 1.0;
filt->gain = 0.8;
filt->fgg = 0.0116;
filt->fgk = 0.461;
filt->len = 162;
filt->down = down;
filt->up = up;
filt->channels = channels;
if (down > up) {
filt->fgg *= (double) up / down;
filt->fgk *= (double) up / down;
filt->len = filt->len * down / up;
}
make_coe(filt);
filt->lag = (filt->len - 1) * channels;
filt->insize = channels*filt->len + filt->lag;
filt->outsize = channels*filt->len;
filt->sin = cst_alloc(int, filt->insize);
filt->sout = cst_alloc(int, filt->outsize);
return filt;
}
void
delete_rateconv(cst_rateconv *filt)
{
cst_free(filt->coep);
cst_free(filt->sin);
cst_free(filt->sout);
cst_free(filt);
}