1 /* THIS IS A MODIFIED VERSION. It was modified by David
2 Huggins-Daines <dhd@cepstral.com> in December 2001 for use in the
3 Flite speech synthesis systems. */
4
5 /*
6 *
7 * RATECONV.C
8 *
9 * Convert sampling rate stdin to stdout
10 *
11 * Copyright (c) 1992, 1995 by Markus Mummert
12 *
13 * Redistribution and use of this software, modifcation and inclusion
14 * into other forms of software are permitted provided that the following
15 * conditions are met:
16 *
17 * 1. Redistributions of this software must retain the above copyright
18 * notice, this list of conditions and the following disclaimer.
19 * 2. If this software is redistributed in a modified condition
20 * it must reveal clearly that it has been modified.
21 *
22 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS''
23 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
24 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
25 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR
26 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
28 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
29 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
30 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
31 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
32 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
33 * DAMAGE.
34 *
35 *
36 * history: 2.9.92 begin coding
37 * 5.9.92 fully operational
38 * 14.2.95 provide BIG_ENDIAN, SWAPPED_BYTES_DEFAULT
39 * switches, Copyright note and References
40 * 25.11.95 changed XXX_ENDIAN to I_AM_XXX_ENDIAN
41 * default gain set to 0.8
42 * 3.12.95 stereo implementation
43 * SWAPPED_BYTES_DEFAULT now HBYTE1ST_DEFAULT
44 * changed [L/2] to (L-1)/2 for exact symmetry
45 *
46 *
47 * IMPLEMENTATION NOTES
48 *
49 * Converting is achieved by interpolating the input samples in
50 * order to evaluate the represented continuous input slope at
51 * sample instances of the new rate (resampling). It is implemented
52 * as a polyphase FIR-filtering process (see reference). The rate
53 * conversion factor is determined by a rational factor. Its
54 * nominator and denominator are integers of almost arbitrary
55 * value, limited only by coefficient memory size.
56 *
57 * General rate conversion formula:
58 *
59 * out(n*Tout) = SUM in(m*Tin) * g((n*d/u-m)*Tin) * Tin
60 * over all m
61 *
62 * FIR-based rate conversion formula for polyphase processing:
63 *
64 * L-1
65 * out(n*Tout) = SUM in(A(i,n)*Tin) * g(B(i,n)*Tin) * Tin
66 * i=0
67 *
68 * A(i,n) = i - (L-1)/2 + [n*d/u]
69 * = i - (L-1)/2 + [(n%u)*d/u] + [n/u]*d
70 * B(i,n) = n*d/u - [n*d/u] + (L-1)/2 - i
71 * = ((n%u)*d/u)%1 + (L-1)/2 - i
72 * Tout = Tin * d/u
73 *
74 * where:
75 * n,i running integers
76 * out(t) output signal sampled at t=n*Tout
77 * in(t) input signal sampled in intervalls Tin
78 * u,d up- and downsampling factor, integers
79 * g(t) interpolating function
80 * L FIR-length of realized g(t), integer
81 * / float-division-operator
82 * % float-modulo-operator
83 * [] integer-operator
84 *
85 * note:
86 * (L-1)/2 in A(i,n) can be omitted as pure time shift yielding
87 * a causal design with a delay of ((L-1)/2)*Tin.
88 * n%u is a cyclic modulo-u counter clocked by out-rate
89 * [n/u]*d is a d-increment counter, advanced when n%u resets
90 * B(i,n)*Tin can take on L*u differnt values, at which g(t)
91 * has to be sampled as a coefficient array
92 *
93 * Interpolation function design:
94 *
95 * The interpolation function design is based on a sinc-function
96 * windowed by a gaussian function. The former determines the
97 * cutoff frequency, the latter limits the necessary FIR-length by
98 * pushing the outer skirts of the resulting impulse response below
99 * a certain threshold fast enough. The drawback is a smoothed
100 * cutoff inducing some aliasing. Due to the symmetry of g(t) the
101 * group delay of the filtering process is contant (linear phase).
102 *
103 * g(t) = 2*fgK*sinc(pi*2*fgK*t) * exp(-pi*(2*fgG*t)**2)
104 *
105 * where:
106 * fgK cutoff frequency of sinc function in f-domain
107 * fgG key frequency of gaussian window in f-domain
108 * reflecting the 6.82dB-down point
109 *
110 * note:
111 * Taking fsin=1/Tin as the input sampling frequncy, it turns out
112 * that in conjunction with L, u and d only the ratios fgK/(fsin/2)
113 * and fgG/(fsin/2) specify the whole proces. Requiring fsin, fgK
114 * and fgG as input purposely keeps the notion of absolute
115 * frequencies.
116 *
117 * Numerical design:
118 *
119 * Samples are expected to be 16bit-signed integers, alternating
120 * left and right channel in case of stereo mode- The byte order
121 * per sample is selectable. FIR-filtering is implemented using
122 * 32bit-signed integer arithmetic. Coefficients are scaled to
123 * find the output sample in the high word of accumulated FIR-sum.
124 *
125 * Interpolation can lead to sample magnitudes exceeding the
126 * input maximum. Worst case is a full scale step function on the
127 * input. In this case the sinc-function exhibits an overshoot of
128 * 2*9=18percent (Gibb's phaenomenon). In any case sample overflow
129 * can be avoided by a gain of 0.8.
130 *
131 * If u=d=1 and if the input stream contains only a single sample,
132 * the whole length of the FIR-filter will be written to the output.
133 * In general the resulting output signal will be (L-1)*fsout/fsin
134 * samples longer than the input signal. The effect is that a
135 * finite input sequence is viewed as padded with zeros before the
136 * `beginning' and after the `end'.
137 *
138 * The output lags ((L-1)/2)*Tin behind to implement g(t) as a
139 * causal system corresponding to a causal relationship of the
140 * discrete-time sequences in(m/fsin) and out(n/fsout) with
141 * resepect to a sequence time origin at t=n*Tin=m*Tout=0.
142 *
143 *
144 * REFERENCES
145 *
146 * Crochiere, R. E., Rabiner, L. R.: "Multirate Digital Signal
147 * Processing", Prentice-Hall, Englewood Cliffs, New Jersey, 1983
148 *
149 * Zwicker, E., Fastl, H.: "Psychoacoustics - Facts and Models",
150 * Springer-Verlag, Berlin, Heidelberg, New-York, Tokyo, 1990 */
151
152 #include "cst_string.h"
153 #include "cst_math.h"
154 #include "cst_alloc.h"
155 #include "cst_error.h"
156 #include "cst_wave.h"
157
158 /*
159 * adaptable defines and globals
160 */
161
162 #define DPRINTF(l, x)
163
164 #define FIXSHIFT 16
165 #define FIXMUL (1<<FIXSHIFT)
166
167 #ifndef M_PI
168 #define M_PI 3.1415926535
169 #endif
170
171 #define sqr(a) ((a)*(a))
172
173 /*
174 * FIR-routines, mono and stereo
175 * this is where we need all the MIPS
176 */
177 static void
fir_mono(int * inp,int * coep,int firlen,int * outp)178 fir_mono(int *inp, int *coep, int firlen, int *outp)
179 {
180 register int akku = 0, *endp;
181 int n1 = (firlen / 8) * 8, n0 = firlen % 8;
182
183 endp = coep + n1;
184 while (coep != endp) {
185 akku += *inp++ * *coep++;
186 akku += *inp++ * *coep++;
187 akku += *inp++ * *coep++;
188 akku += *inp++ * *coep++;
189 akku += *inp++ * *coep++;
190 akku += *inp++ * *coep++;
191 akku += *inp++ * *coep++;
192 akku += *inp++ * *coep++;
193 }
194
195 endp = coep + n0;
196 while (coep != endp) {
197 akku += *inp++ * *coep++;
198 }
199 *outp = akku;
200 }
201
202 static void
fir_stereo(int * inp,int * coep,int firlen,int * out1p,int * out2p)203 fir_stereo(int *inp, int *coep,
204 int firlen, int *out1p, int *out2p)
205 {
206 register int akku1 = 0, akku2 = 0, *endp;
207 int n1 = (firlen / 8) * 8, n0 = firlen % 8;
208
209 endp = coep + n1;
210 while (coep != endp) {
211 akku1 += *inp++ * *coep;
212 akku2 += *inp++ * *coep++;
213 akku1 += *inp++ * *coep;
214 akku2 += *inp++ * *coep++;
215 akku1 += *inp++ * *coep;
216 akku2 += *inp++ * *coep++;
217 akku1 += *inp++ * *coep;
218 akku2 += *inp++ * *coep++;
219 akku1 += *inp++ * *coep;
220 akku2 += *inp++ * *coep++;
221 akku1 += *inp++ * *coep;
222 akku2 += *inp++ * *coep++;
223 akku1 += *inp++ * *coep;
224 akku2 += *inp++ * *coep++;
225 akku1 += *inp++ * *coep;
226 akku2 += *inp++ * *coep++;
227 }
228
229 endp = coep + n0;
230 while (coep != endp) {
231 akku1 += *inp++ * *coep;
232 akku2 += *inp++ * *coep++;
233 }
234 *out1p = akku1;
235 *out2p = akku2;
236 }
237
238 /*
239 * filtering from input buffer to output buffer;
240 * returns number of processed samples in output buffer:
241 * if it is not equal to output buffer size,
242 * the input buffer is expected to be refilled upon entry, so that
243 * the last firlen numbers of the old input buffer are
244 * the first firlen numbers of the new input buffer;
245 * if it is equal to output buffer size, the output buffer
246 * is full and is expected to be stowed away;
247 *
248 */
249 static int
filtering_on_buffers(cst_rateconv * filt)250 filtering_on_buffers(cst_rateconv *filt)
251 {
252 int insize;
253
254 DPRINTF(0, ("filtering_on_buffers(%d)\n", filt->incount));
255 insize = filt->incount + filt->lag;
256 if (filt->channels == 1) {
257 while (1) {
258 filt->inoffset = (filt->cycctr * filt->down)/filt->up;
259 if ((filt->inbaseidx + filt->inoffset + filt->len) > insize) {
260 filt->inbaseidx -= insize - filt->len + 1;
261 memcpy(filt->sin, filt->sin + insize - filt->lag,
262 filt->lag * sizeof(int));
263 /* Prevent people from re-filtering the same stuff. */
264 filt->incount = 0;
265 return 0;
266 }
267 fir_mono(filt->sin + filt->inoffset + filt->inbaseidx,
268 filt->coep + filt->cycctr * filt->len,
269 filt->len, filt->sout + filt->outidx);
270 DPRINTF(1, ("in(%d + %d) = %d cycctr %d out(%d) = %d\n",
271 filt->inoffset, filt->inbaseidx,
272 filt->sin[filt->inoffset + filt->inbaseidx],
273 filt->cycctr, filt->outidx,
274 filt->sout[filt->outidx] >> FIXSHIFT));
275 ++filt->outidx;
276 ++filt->cycctr;
277 if (!(filt->cycctr %= filt->up))
278 filt->inbaseidx += filt->down;
279 if (!(filt->outidx %= filt->outsize))
280 return filt->outsize;
281 }
282 } else if (filt->channels == 2) {
283 /*
284 * rule how to convert mono routine to stereo routine:
285 * firlen, up, down and cycctr relate to samples in general,
286 * wether mono or stereo; inbaseidx, inoffset and outidx as
287 * well as insize and outsize still account for mono samples.
288 */
289 while (1) {
290 filt->inoffset = 2*((filt->cycctr * filt->down)/filt->up);
291 if ((filt->inbaseidx + filt->inoffset + 2*filt->len) > insize) {
292 filt->inbaseidx -= insize - 2*filt->len + 2;
293 return filt->outidx;
294 }
295 fir_stereo(filt->sin + filt->inoffset + filt->inbaseidx,
296 filt->coep + filt->cycctr * filt->len,
297 filt->len,
298 filt->sout + filt->outidx,
299 filt->sout + filt->outidx + 1);
300 filt->outidx += 2;
301 ++filt->cycctr;
302 if (!(filt->cycctr %= filt->up))
303 filt->inbaseidx += 2*filt->down;
304 if (!(filt->outidx %= filt->outsize))
305 return filt->outsize;
306 }
307 } else {
308 cst_errmsg("filtering_on_buffers: only 1 or 2 channels supported!\n");
309 cst_error();
310 }
311 return 0;
312 }
313
314 /*
315 * convert buffer of n samples to ints
316 */
317 static void
sample_to_int(short * buff,int n)318 sample_to_int(short *buff, int n)
319 {
320 short *s, *e;
321 int *d;
322
323 if (n < 1)
324 return;
325 s = buff + n;
326 d = (int*)buff + n;
327 e = buff;
328 while (s != e) {
329 *--d = (int)(*--s);
330 }
331 }
332
333 /*
334 * convert buffer of n ints to samples
335 */
336 static void
int_to_sample(short * buff,int n)337 int_to_sample(short *buff, int n)
338 {
339 int *s;
340 short *e, *d;
341
342 if (n < 1)
343 return;
344 s = (int *)buff;
345 d = buff;
346 e = buff + n;
347 while (d != e)
348 *d++ = (*s++ >> FIXSHIFT);
349 }
350
351 /*
352 * read and convert input sample format to integer
353 */
354 int
cst_rateconv_in(cst_rateconv * filt,const short * inptr,int max)355 cst_rateconv_in(cst_rateconv *filt, const short *inptr, int max)
356 {
357 if (max > filt->insize - filt->lag)
358 max = filt->insize - filt->lag;
359 if (max > 0) {
360 memcpy(filt->sin + filt->lag, inptr, max * sizeof(short));
361 sample_to_int((short *)(filt->sin + filt->lag), max);
362 }
363 filt->incount = max;
364 return max;
365 }
366
367 /*
368 * do some conversion jobs and write
369 */
370 int
cst_rateconv_out(cst_rateconv * filt,short * outptr,int max)371 cst_rateconv_out(cst_rateconv *filt, short *outptr, int max)
372 {
373 int outsize;
374
375 if ((outsize = filtering_on_buffers(filt)) == 0)
376 return 0;
377 if (max > outsize)
378 max = outsize;
379 int_to_sample((short *)filt->sout, max);
380 memcpy(outptr, filt->sout, max * sizeof(short));
381 return max;
382 }
383
384 int
cst_rateconv_leadout(cst_rateconv * filt)385 cst_rateconv_leadout(cst_rateconv *filt)
386 {
387 memset(filt->sin + filt->lag, 0, filt->lag * sizeof(int));
388 filt->incount = filt->lag;
389 return filt->lag;
390 }
391
392 /*
393 * evaluate sinc(x) = sin(x)/x safely
394 */
395 static double
sinc(double x)396 sinc(double x)
397 {
398 return(fabs(x) < 1E-50 ? 1.0 : sin(fmod(x,2*M_PI))/x);
399 }
400
401 /*
402 * evaluate interpolation function g(t) at t
403 * integral of g(t) over all times is expected to be one
404 */
405 static double
interpol_func(double t,double fgk,double fgg)406 interpol_func(double t, double fgk, double fgg)
407 {
408 return (2*fgk*sinc(M_PI*2*fgk*t)*exp(-M_PI*sqr(2*fgg*t)));
409 }
410
411 /*
412 * evaluate coefficient from i, q=n%u by sampling interpolation function
413 * and scale it for integer multiplication used by FIR-filtering
414 */
415 static int
coefficient(int i,int q,cst_rateconv * filt)416 coefficient(int i, int q, cst_rateconv *filt)
417 {
418 return (int)(FIXMUL * filt->gain *
419 interpol_func((fmod(q* filt->down/(double)filt->up,1.0)
420 + (filt->len-1)/2.0 - i) / filt->fsin,
421 filt->fgk, filt->fgg) / filt->fsin);
422 }
423
424 /*
425 * set up coefficient array
426 */
427 static void
make_coe(cst_rateconv * filt)428 make_coe(cst_rateconv *filt)
429 {
430 int i, q;
431
432 filt->coep = cst_alloc(int, filt->len * filt->up);
433 for (i = 0; i < filt->len; i++) {
434 for (q = 0; q < filt->up; q++) {
435 filt->coep[q * filt->len + i]
436 = coefficient(i, q, filt);
437 }
438 }
439 }
440
441 cst_rateconv *
new_rateconv(int up,int down,int channels)442 new_rateconv(int up, int down, int channels)
443 {
444 cst_rateconv *filt;
445
446 if (!(channels == 1 || channels == 2)) {
447 cst_errmsg("new_rateconv: channels must be 1 or 2\n");
448 cst_error();
449 }
450
451 filt = cst_alloc(cst_rateconv, 1);
452 filt->fsin = 1.0;
453 filt->gain = 0.8;
454 filt->fgg = 0.0116;
455 filt->fgk = 0.461;
456 filt->len = 162;
457 filt->down = down;
458 filt->up = up;
459 filt->channels = channels;
460
461 if (down > up) {
462 filt->fgg *= (double) up / down;
463 filt->fgk *= (double) up / down;
464 filt->len = filt->len * down / up;
465 }
466
467 make_coe(filt);
468 filt->lag = (filt->len - 1) * channels;
469 filt->insize = channels*filt->len + filt->lag;
470 filt->outsize = channels*filt->len;
471 filt->sin = cst_alloc(int, filt->insize);
472 filt->sout = cst_alloc(int, filt->outsize);
473
474 return filt;
475 }
476
477 void
delete_rateconv(cst_rateconv * filt)478 delete_rateconv(cst_rateconv *filt)
479 {
480 cst_free(filt->coep);
481 cst_free(filt->sin);
482 cst_free(filt->sout);
483 cst_free(filt);
484 }
485