| /* 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); |
| } |