/* 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 voidfir_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 voidfir_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 intfiltering_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 voidsample_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 voidint_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 */intcst_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 */intcst_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;} intcst_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 doublesinc(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 doubleinterpol_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 intcoefficient(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 voidmake_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;} voiddelete_rateconv(cst_rateconv *filt){	cst_free(filt->coep);	cst_free(filt->sin);	cst_free(filt->sout);	cst_free(filt);}