Fixes a crash on destroy with the pull API
[libopusenc.git] / src / resample.c
1 /* Copyright (C) 2007-2008 Jean-Marc Valin
2    Copyright (C) 2008      Thorvald Natvig
3
4    File: resample.c
5    Arbitrary resampling code
6
7    Redistribution and use in source and binary forms, with or without
8    modification, are permitted provided that the following conditions are
9    met:
10
11    1. Redistributions of source code must retain the above copyright notice,
12    this list of conditions and the following disclaimer.
13
14    2. Redistributions in binary form must reproduce the above copyright
15    notice, this list of conditions and the following disclaimer in the
16    documentation and/or other materials provided with the distribution.
17
18    3. The name of the author may not be used to endorse or promote products
19    derived from this software without specific prior written permission.
20
21    THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
22    IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
23    OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24    DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
25    INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
26    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27    SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28    HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
29    STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30    ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31    POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 /*
35    The design goals of this code are:
36       - Very fast algorithm
37       - SIMD-friendly algorithm
38       - Low memory requirement
39       - Good *perceptual* quality (and not best SNR)
40
41    Warning: This resampler is relatively new. Although I think I got rid of
42    all the major bugs and I don't expect the API to change anymore, there
43    may be something I've missed. So use with caution.
44
45    This algorithm is based on this original resampling algorithm:
46    Smith, Julius O. Digital Audio Resampling Home Page
47    Center for Computer Research in Music and Acoustics (CCRMA),
48    Stanford University, 2007.
49    Web published at http://www-ccrma.stanford.edu/~jos/resample/.
50
51    There is one main difference, though. This resampler uses cubic
52    interpolation instead of linear interpolation in the above paper. This
53    makes the table much smaller and makes it possible to compute that table
54    on a per-stream basis. In turn, being able to tweak the table for each
55    stream makes it possible to both reduce complexity on simple ratios
56    (e.g. 2/3), and get rid of the rounding operations in the inner loop.
57    The latter both reduces CPU time and makes the algorithm more SIMD-friendly.
58 */
59
60 #ifdef HAVE_CONFIG_H
61 #include "config.h"
62 #endif
63
64 #ifdef OUTSIDE_SPEEX
65 #include <stdlib.h>
66 static void *speex_alloc (int size) {return calloc(size,1);}
67 static void *speex_realloc (void *ptr, int size) {return realloc(ptr, size);}
68 static void speex_free (void *ptr) {free(ptr);}
69 #define speex_assert(x)
70 #define EXPORT
71 #include "speex_resampler.h"
72 #include "arch.h"
73 #else /* OUTSIDE_SPEEX */
74
75 #include "speex/speex_resampler.h"
76 #include "arch.h"
77 #include "os_support.h"
78 #endif /* OUTSIDE_SPEEX */
79
80 #include "stack_alloc.h"
81 #include <math.h>
82 #include <limits.h>
83
84 #ifndef M_PI
85 #define M_PI 3.14159265358979323846
86 #endif
87
88 #define IMAX(a,b) ((a) > (b) ? (a) : (b))
89 #define IMIN(a,b) ((a) < (b) ? (a) : (b))
90
91 #ifndef NULL
92 #define NULL 0
93 #endif
94
95 #ifndef UINT32_MAX
96 #define UINT32_MAX 4294967296U
97 #endif
98
99 #if defined(__SSE__) && !defined(FIXED_POINT)
100 #include "resample_sse.h"
101 #endif
102
103 #ifdef _USE_NEON
104 #include "resample_neon.h"
105 #endif
106
107 /* Numer of elements to allocate on the stack */
108 #ifdef VAR_ARRAYS
109 #define FIXED_STACK_ALLOC 8192
110 #else
111 #define FIXED_STACK_ALLOC 1024
112 #endif
113
114 typedef int (*resampler_basic_func)(SpeexResamplerState *, spx_uint32_t , const spx_word16_t *, spx_uint32_t *, spx_word16_t *, spx_uint32_t *);
115
116 struct SpeexResamplerState_ {
117    spx_uint32_t in_rate;
118    spx_uint32_t out_rate;
119    spx_uint32_t num_rate;
120    spx_uint32_t den_rate;
121
122    int    quality;
123    spx_uint32_t nb_channels;
124    spx_uint32_t filt_len;
125    spx_uint32_t mem_alloc_size;
126    spx_uint32_t buffer_size;
127    int          int_advance;
128    int          frac_advance;
129    float  cutoff;
130    spx_uint32_t oversample;
131    int          initialised;
132    int          started;
133
134    /* These are per-channel */
135    spx_int32_t  *last_sample;
136    spx_uint32_t *samp_frac_num;
137    spx_uint32_t *magic_samples;
138
139    spx_word16_t *mem;
140    spx_word16_t *sinc_table;
141    spx_uint32_t sinc_table_length;
142    resampler_basic_func resampler_ptr;
143
144    int    in_stride;
145    int    out_stride;
146 } ;
147
148 static const double kaiser12_table[68] = {
149    0.99859849, 1.00000000, 0.99859849, 0.99440475, 0.98745105, 0.97779076,
150    0.96549770, 0.95066529, 0.93340547, 0.91384741, 0.89213598, 0.86843014,
151    0.84290116, 0.81573067, 0.78710866, 0.75723148, 0.72629970, 0.69451601,
152    0.66208321, 0.62920216, 0.59606986, 0.56287762, 0.52980938, 0.49704014,
153    0.46473455, 0.43304576, 0.40211431, 0.37206735, 0.34301800, 0.31506490,
154    0.28829195, 0.26276832, 0.23854851, 0.21567274, 0.19416736, 0.17404546,
155    0.15530766, 0.13794294, 0.12192957, 0.10723616, 0.09382272, 0.08164178,
156    0.07063950, 0.06075685, 0.05193064, 0.04409466, 0.03718069, 0.03111947,
157    0.02584161, 0.02127838, 0.01736250, 0.01402878, 0.01121463, 0.00886058,
158    0.00691064, 0.00531256, 0.00401805, 0.00298291, 0.00216702, 0.00153438,
159    0.00105297, 0.00069463, 0.00043489, 0.00025272, 0.00013031, 0.0000527734,
160    0.00001000, 0.00000000};
161 /*
162 static const double kaiser12_table[36] = {
163    0.99440475, 1.00000000, 0.99440475, 0.97779076, 0.95066529, 0.91384741,
164    0.86843014, 0.81573067, 0.75723148, 0.69451601, 0.62920216, 0.56287762,
165    0.49704014, 0.43304576, 0.37206735, 0.31506490, 0.26276832, 0.21567274,
166    0.17404546, 0.13794294, 0.10723616, 0.08164178, 0.06075685, 0.04409466,
167    0.03111947, 0.02127838, 0.01402878, 0.00886058, 0.00531256, 0.00298291,
168    0.00153438, 0.00069463, 0.00025272, 0.0000527734, 0.00000500, 0.00000000};
169 */
170 static const double kaiser10_table[36] = {
171    0.99537781, 1.00000000, 0.99537781, 0.98162644, 0.95908712, 0.92831446,
172    0.89005583, 0.84522401, 0.79486424, 0.74011713, 0.68217934, 0.62226347,
173    0.56155915, 0.50119680, 0.44221549, 0.38553619, 0.33194107, 0.28205962,
174    0.23636152, 0.19515633, 0.15859932, 0.12670280, 0.09935205, 0.07632451,
175    0.05731132, 0.04193980, 0.02979584, 0.02044510, 0.01345224, 0.00839739,
176    0.00488951, 0.00257636, 0.00115101, 0.00035515, 0.00000000, 0.00000000};
177
178 static const double kaiser8_table[36] = {
179    0.99635258, 1.00000000, 0.99635258, 0.98548012, 0.96759014, 0.94302200,
180    0.91223751, 0.87580811, 0.83439927, 0.78875245, 0.73966538, 0.68797126,
181    0.63451750, 0.58014482, 0.52566725, 0.47185369, 0.41941150, 0.36897272,
182    0.32108304, 0.27619388, 0.23465776, 0.19672670, 0.16255380, 0.13219758,
183    0.10562887, 0.08273982, 0.06335451, 0.04724088, 0.03412321, 0.02369490,
184    0.01563093, 0.00959968, 0.00527363, 0.00233883, 0.00050000, 0.00000000};
185
186 static const double kaiser6_table[36] = {
187    0.99733006, 1.00000000, 0.99733006, 0.98935595, 0.97618418, 0.95799003,
188    0.93501423, 0.90755855, 0.87598009, 0.84068475, 0.80211977, 0.76076565,
189    0.71712752, 0.67172623, 0.62508937, 0.57774224, 0.53019925, 0.48295561,
190    0.43647969, 0.39120616, 0.34752997, 0.30580127, 0.26632152, 0.22934058,
191    0.19505503, 0.16360756, 0.13508755, 0.10953262, 0.08693120, 0.06722600,
192    0.05031820, 0.03607231, 0.02432151, 0.01487334, 0.00752000, 0.00000000};
193
194 struct FuncDef {
195    const double *table;
196    int oversample;
197 };
198
199 static const struct FuncDef _KAISER12 = {kaiser12_table, 64};
200 #define KAISER12 (&_KAISER12)
201 /*static struct FuncDef _KAISER12 = {kaiser12_table, 32};
202 #define KAISER12 (&_KAISER12)*/
203 static const struct FuncDef _KAISER10 = {kaiser10_table, 32};
204 #define KAISER10 (&_KAISER10)
205 static const struct FuncDef _KAISER8 = {kaiser8_table, 32};
206 #define KAISER8 (&_KAISER8)
207 static const struct FuncDef _KAISER6 = {kaiser6_table, 32};
208 #define KAISER6 (&_KAISER6)
209
210 struct QualityMapping {
211    int base_length;
212    int oversample;
213    float downsample_bandwidth;
214    float upsample_bandwidth;
215    const struct FuncDef *window_func;
216 };
217
218
219 /* This table maps conversion quality to internal parameters. There are two
220    reasons that explain why the up-sampling bandwidth is larger than the
221    down-sampling bandwidth:
222    1) When up-sampling, we can assume that the spectrum is already attenuated
223       close to the Nyquist rate (from an A/D or a previous resampling filter)
224    2) Any aliasing that occurs very close to the Nyquist rate will be masked
225       by the sinusoids/noise just below the Nyquist rate (guaranteed only for
226       up-sampling).
227 */
228 static const struct QualityMapping quality_map[11] = {
229    {  8,  4, 0.830f, 0.860f, KAISER6 }, /* Q0 */
230    { 16,  4, 0.850f, 0.880f, KAISER6 }, /* Q1 */
231    { 32,  4, 0.882f, 0.910f, KAISER6 }, /* Q2 */  /* 82.3% cutoff ( ~60 dB stop) 6  */
232    { 48,  8, 0.895f, 0.917f, KAISER8 }, /* Q3 */  /* 84.9% cutoff ( ~80 dB stop) 8  */
233    { 64,  8, 0.921f, 0.940f, KAISER8 }, /* Q4 */  /* 88.7% cutoff ( ~80 dB stop) 8  */
234    { 80, 16, 0.922f, 0.940f, KAISER10}, /* Q5 */  /* 89.1% cutoff (~100 dB stop) 10 */
235    { 96, 16, 0.940f, 0.945f, KAISER10}, /* Q6 */  /* 91.5% cutoff (~100 dB stop) 10 */
236    {128, 16, 0.950f, 0.950f, KAISER10}, /* Q7 */  /* 93.1% cutoff (~100 dB stop) 10 */
237    {160, 16, 0.960f, 0.960f, KAISER10}, /* Q8 */  /* 94.5% cutoff (~100 dB stop) 10 */
238    {192, 32, 0.968f, 0.968f, KAISER12}, /* Q9 */  /* 95.5% cutoff (~100 dB stop) 10 */
239    {256, 32, 0.975f, 0.975f, KAISER12}, /* Q10 */ /* 96.6% cutoff (~100 dB stop) 10 */
240 };
241 /*8,24,40,56,80,104,128,160,200,256,320*/
242 static double compute_func(float x, const struct FuncDef *func)
243 {
244    float y, frac;
245    double interp[4];
246    int ind;
247    y = x*func->oversample;
248    ind = (int)floor(y);
249    frac = (y-ind);
250    /* CSE with handle the repeated powers */
251    interp[3] =  -0.1666666667*frac + 0.1666666667*(frac*frac*frac);
252    interp[2] = frac + 0.5*(frac*frac) - 0.5*(frac*frac*frac);
253    /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac;*/
254    interp[0] = -0.3333333333*frac + 0.5*(frac*frac) - 0.1666666667*(frac*frac*frac);
255    /* Just to make sure we don't have rounding problems */
256    interp[1] = 1.f-interp[3]-interp[2]-interp[0];
257
258    /*sum = frac*accum[1] + (1-frac)*accum[2];*/
259    return interp[0]*func->table[ind] + interp[1]*func->table[ind+1] + interp[2]*func->table[ind+2] + interp[3]*func->table[ind+3];
260 }
261
262 #if 0
263 #include <stdio.h>
264 int main(int argc, char **argv)
265 {
266    int i;
267    for (i=0;i<256;i++)
268    {
269       printf ("%f\n", compute_func(i/256., KAISER12));
270    }
271    return 0;
272 }
273 #endif
274
275 #ifdef FIXED_POINT
276 /* The slow way of computing a sinc for the table. Should improve that some day */
277 static spx_word16_t sinc(float cutoff, float x, int N, const struct FuncDef *window_func)
278 {
279    /*fprintf (stderr, "%f ", x);*/
280    float xx = x * cutoff;
281    if (fabs(x)<1e-6f)
282       return WORD2INT(32768.*cutoff);
283    else if (fabs(x) > .5f*N)
284       return 0;
285    /*FIXME: Can it really be any slower than this? */
286    return WORD2INT(32768.*cutoff*sin(M_PI*xx)/(M_PI*xx) * compute_func(fabs(2.*x/N), window_func));
287 }
288 #else
289 /* The slow way of computing a sinc for the table. Should improve that some day */
290 static spx_word16_t sinc(float cutoff, float x, int N, const struct FuncDef *window_func)
291 {
292    /*fprintf (stderr, "%f ", x);*/
293    float xx = x * cutoff;
294    if (fabs(x)<1e-6)
295       return cutoff;
296    else if (fabs(x) > .5*N)
297       return 0;
298    /*FIXME: Can it really be any slower than this? */
299    return cutoff*sin(M_PI*xx)/(M_PI*xx) * compute_func(fabs(2.*x/N), window_func);
300 }
301 #endif
302
303 #ifdef FIXED_POINT
304 static void cubic_coef(spx_word16_t x, spx_word16_t interp[4])
305 {
306    /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
307    but I know it's MMSE-optimal on a sinc */
308    spx_word16_t x2, x3;
309    x2 = MULT16_16_P15(x, x);
310    x3 = MULT16_16_P15(x, x2);
311    interp[0] = PSHR32(MULT16_16(QCONST16(-0.16667f, 15),x) + MULT16_16(QCONST16(0.16667f, 15),x3),15);
312    interp[1] = EXTRACT16(EXTEND32(x) + SHR32(SUB32(EXTEND32(x2),EXTEND32(x3)),1));
313    interp[3] = PSHR32(MULT16_16(QCONST16(-0.33333f, 15),x) + MULT16_16(QCONST16(.5f,15),x2) - MULT16_16(QCONST16(0.16667f, 15),x3),15);
314    /* Just to make sure we don't have rounding problems */
315    interp[2] = Q15_ONE-interp[0]-interp[1]-interp[3];
316    if (interp[2]<32767)
317       interp[2]+=1;
318 }
319 #else
320 static void cubic_coef(spx_word16_t frac, spx_word16_t interp[4])
321 {
322    /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
323    but I know it's MMSE-optimal on a sinc */
324    interp[0] =  -0.16667f*frac + 0.16667f*frac*frac*frac;
325    interp[1] = frac + 0.5f*frac*frac - 0.5f*frac*frac*frac;
326    /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac;*/
327    interp[3] = -0.33333f*frac + 0.5f*frac*frac - 0.16667f*frac*frac*frac;
328    /* Just to make sure we don't have rounding problems */
329    interp[2] = 1.-interp[0]-interp[1]-interp[3];
330 }
331 #endif
332
333 static int resampler_basic_direct_single(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
334 {
335    const int N = st->filt_len;
336    int out_sample = 0;
337    int last_sample = st->last_sample[channel_index];
338    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
339    const spx_word16_t *sinc_table = st->sinc_table;
340    const int out_stride = st->out_stride;
341    const int int_advance = st->int_advance;
342    const int frac_advance = st->frac_advance;
343    const spx_uint32_t den_rate = st->den_rate;
344    spx_word32_t sum;
345
346    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
347    {
348       const spx_word16_t *sinct = & sinc_table[samp_frac_num*N];
349       const spx_word16_t *iptr = & in[last_sample];
350
351 #ifndef OVERRIDE_INNER_PRODUCT_SINGLE
352       int j;
353       sum = 0;
354       for(j=0;j<N;j++) sum += MULT16_16(sinct[j], iptr[j]);
355
356 /*    This code is slower on most DSPs which have only 2 accumulators.
357       Plus this this forces truncation to 32 bits and you lose the HW guard bits.
358       I think we can trust the compiler and let it vectorize and/or unroll itself.
359       spx_word32_t accum[4] = {0,0,0,0};
360       for(j=0;j<N;j+=4) {
361         accum[0] += MULT16_16(sinct[j], iptr[j]);
362         accum[1] += MULT16_16(sinct[j+1], iptr[j+1]);
363         accum[2] += MULT16_16(sinct[j+2], iptr[j+2]);
364         accum[3] += MULT16_16(sinct[j+3], iptr[j+3]);
365       }
366       sum = accum[0] + accum[1] + accum[2] + accum[3];
367 */
368       sum = SATURATE32PSHR(sum, 15, 32767);
369 #else
370       sum = inner_product_single(sinct, iptr, N);
371 #endif
372
373       out[out_stride * out_sample++] = sum;
374       last_sample += int_advance;
375       samp_frac_num += frac_advance;
376       if (samp_frac_num >= den_rate)
377       {
378          samp_frac_num -= den_rate;
379          last_sample++;
380       }
381    }
382
383    st->last_sample[channel_index] = last_sample;
384    st->samp_frac_num[channel_index] = samp_frac_num;
385    return out_sample;
386 }
387
388 #ifdef FIXED_POINT
389 #else
390 /* This is the same as the previous function, except with a double-precision accumulator */
391 static int resampler_basic_direct_double(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
392 {
393    const int N = st->filt_len;
394    int out_sample = 0;
395    int last_sample = st->last_sample[channel_index];
396    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
397    const spx_word16_t *sinc_table = st->sinc_table;
398    const int out_stride = st->out_stride;
399    const int int_advance = st->int_advance;
400    const int frac_advance = st->frac_advance;
401    const spx_uint32_t den_rate = st->den_rate;
402    double sum;
403
404    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
405    {
406       const spx_word16_t *sinct = & sinc_table[samp_frac_num*N];
407       const spx_word16_t *iptr = & in[last_sample];
408
409 #ifndef OVERRIDE_INNER_PRODUCT_DOUBLE
410       int j;
411       double accum[4] = {0,0,0,0};
412
413       for(j=0;j<N;j+=4) {
414         accum[0] += sinct[j]*iptr[j];
415         accum[1] += sinct[j+1]*iptr[j+1];
416         accum[2] += sinct[j+2]*iptr[j+2];
417         accum[3] += sinct[j+3]*iptr[j+3];
418       }
419       sum = accum[0] + accum[1] + accum[2] + accum[3];
420 #else
421       sum = inner_product_double(sinct, iptr, N);
422 #endif
423
424       out[out_stride * out_sample++] = PSHR32(sum, 15);
425       last_sample += int_advance;
426       samp_frac_num += frac_advance;
427       if (samp_frac_num >= den_rate)
428       {
429          samp_frac_num -= den_rate;
430          last_sample++;
431       }
432    }
433
434    st->last_sample[channel_index] = last_sample;
435    st->samp_frac_num[channel_index] = samp_frac_num;
436    return out_sample;
437 }
438 #endif
439
440 static int resampler_basic_interpolate_single(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
441 {
442    const int N = st->filt_len;
443    int out_sample = 0;
444    int last_sample = st->last_sample[channel_index];
445    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
446    const int out_stride = st->out_stride;
447    const int int_advance = st->int_advance;
448    const int frac_advance = st->frac_advance;
449    const spx_uint32_t den_rate = st->den_rate;
450    spx_word32_t sum;
451
452    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
453    {
454       const spx_word16_t *iptr = & in[last_sample];
455
456       const int offset = samp_frac_num*st->oversample/st->den_rate;
457 #ifdef FIXED_POINT
458       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
459 #else
460       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
461 #endif
462       spx_word16_t interp[4];
463
464
465 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_SINGLE
466       int j;
467       spx_word32_t accum[4] = {0,0,0,0};
468
469       for(j=0;j<N;j++) {
470         const spx_word16_t curr_in=iptr[j];
471         accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
472         accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
473         accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
474         accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
475       }
476
477       cubic_coef(frac, interp);
478       sum = MULT16_32_Q15(interp[0],SHR32(accum[0], 1)) + MULT16_32_Q15(interp[1],SHR32(accum[1], 1)) + MULT16_32_Q15(interp[2],SHR32(accum[2], 1)) + MULT16_32_Q15(interp[3],SHR32(accum[3], 1));
479       sum = SATURATE32PSHR(sum, 15, 32767);
480 #else
481       cubic_coef(frac, interp);
482       sum = interpolate_product_single(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
483 #endif
484
485       out[out_stride * out_sample++] = sum;
486       last_sample += int_advance;
487       samp_frac_num += frac_advance;
488       if (samp_frac_num >= den_rate)
489       {
490          samp_frac_num -= den_rate;
491          last_sample++;
492       }
493    }
494
495    st->last_sample[channel_index] = last_sample;
496    st->samp_frac_num[channel_index] = samp_frac_num;
497    return out_sample;
498 }
499
500 #ifdef FIXED_POINT
501 #else
502 /* This is the same as the previous function, except with a double-precision accumulator */
503 static int resampler_basic_interpolate_double(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
504 {
505    const int N = st->filt_len;
506    int out_sample = 0;
507    int last_sample = st->last_sample[channel_index];
508    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
509    const int out_stride = st->out_stride;
510    const int int_advance = st->int_advance;
511    const int frac_advance = st->frac_advance;
512    const spx_uint32_t den_rate = st->den_rate;
513    spx_word32_t sum;
514
515    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
516    {
517       const spx_word16_t *iptr = & in[last_sample];
518
519       const int offset = samp_frac_num*st->oversample/st->den_rate;
520 #ifdef FIXED_POINT
521       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
522 #else
523       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
524 #endif
525       spx_word16_t interp[4];
526
527
528 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE
529       int j;
530       double accum[4] = {0,0,0,0};
531
532       for(j=0;j<N;j++) {
533         const double curr_in=iptr[j];
534         accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
535         accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
536         accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
537         accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
538       }
539
540       cubic_coef(frac, interp);
541       sum = MULT16_32_Q15(interp[0],accum[0]) + MULT16_32_Q15(interp[1],accum[1]) + MULT16_32_Q15(interp[2],accum[2]) + MULT16_32_Q15(interp[3],accum[3]);
542 #else
543       cubic_coef(frac, interp);
544       sum = interpolate_product_double(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
545 #endif
546
547       out[out_stride * out_sample++] = PSHR32(sum,15);
548       last_sample += int_advance;
549       samp_frac_num += frac_advance;
550       if (samp_frac_num >= den_rate)
551       {
552          samp_frac_num -= den_rate;
553          last_sample++;
554       }
555    }
556
557    st->last_sample[channel_index] = last_sample;
558    st->samp_frac_num[channel_index] = samp_frac_num;
559    return out_sample;
560 }
561 #endif
562
563 /* This resampler is used to produce zero output in situations where memory
564    for the filter could not be allocated.  The expected numbers of input and
565    output samples are still processed so that callers failing to check error
566    codes are not surprised, possibly getting into infinite loops. */
567 static int resampler_basic_zero(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_word16_t *in, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
568 {
569    int out_sample = 0;
570    int last_sample = st->last_sample[channel_index];
571    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
572    const int out_stride = st->out_stride;
573    const int int_advance = st->int_advance;
574    const int frac_advance = st->frac_advance;
575    const spx_uint32_t den_rate = st->den_rate;
576
577    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
578    {
579       out[out_stride * out_sample++] = 0;
580       last_sample += int_advance;
581       samp_frac_num += frac_advance;
582       if (samp_frac_num >= den_rate)
583       {
584          samp_frac_num -= den_rate;
585          last_sample++;
586       }
587    }
588
589    st->last_sample[channel_index] = last_sample;
590    st->samp_frac_num[channel_index] = samp_frac_num;
591    return out_sample;
592 }
593
594 static int _muldiv(spx_uint32_t *result, spx_uint32_t value, spx_uint32_t mul, spx_uint32_t div)
595 {
596    speex_assert(result);
597    spx_uint32_t major = value / div;
598    spx_uint32_t remainder = value % div;
599    /* TODO: Could use 64 bits operation to check for overflow. But only guaranteed in C99+ */
600    if (remainder > UINT32_MAX / mul || major > UINT32_MAX / mul
601        || major * mul > UINT32_MAX - remainder * mul / div)
602       return RESAMPLER_ERR_OVERFLOW;
603    *result = remainder * mul / div + major * mul;
604    return RESAMPLER_ERR_SUCCESS;
605 }
606
607 static int update_filter(SpeexResamplerState *st)
608 {
609    spx_uint32_t old_length = st->filt_len;
610    spx_uint32_t old_alloc_size = st->mem_alloc_size;
611    int use_direct;
612    spx_uint32_t min_sinc_table_length;
613    spx_uint32_t min_alloc_size;
614
615    st->int_advance = st->num_rate/st->den_rate;
616    st->frac_advance = st->num_rate%st->den_rate;
617    st->oversample = quality_map[st->quality].oversample;
618    st->filt_len = quality_map[st->quality].base_length;
619
620    if (st->num_rate > st->den_rate)
621    {
622       /* down-sampling */
623       st->cutoff = quality_map[st->quality].downsample_bandwidth * st->den_rate / st->num_rate;
624       if (_muldiv(&st->filt_len,st->filt_len,st->num_rate,st->den_rate) != RESAMPLER_ERR_SUCCESS)
625          goto fail;
626       /* Round up to make sure we have a multiple of 8 for SSE */
627       st->filt_len = ((st->filt_len-1)&(~0x7))+8;
628       if (2*st->den_rate < st->num_rate)
629          st->oversample >>= 1;
630       if (4*st->den_rate < st->num_rate)
631          st->oversample >>= 1;
632       if (8*st->den_rate < st->num_rate)
633          st->oversample >>= 1;
634       if (16*st->den_rate < st->num_rate)
635          st->oversample >>= 1;
636       if (st->oversample < 1)
637          st->oversample = 1;
638    } else {
639       /* up-sampling */
640       st->cutoff = quality_map[st->quality].upsample_bandwidth;
641    }
642
643    /* Choose the resampling type that requires the least amount of memory */
644 #ifdef RESAMPLE_FULL_SINC_TABLE
645    use_direct = 1;
646    if (INT_MAX/sizeof(spx_word16_t)/st->den_rate < st->filt_len)
647       goto fail;
648 #else
649    use_direct = st->filt_len*st->den_rate <= st->filt_len*st->oversample+8
650                 && INT_MAX/sizeof(spx_word16_t)/st->den_rate >= st->filt_len;
651 #endif
652    if (use_direct)
653    {
654       min_sinc_table_length = st->filt_len*st->den_rate;
655    } else {
656       if ((INT_MAX/sizeof(spx_word16_t)-8)/st->oversample < st->filt_len)
657          goto fail;
658
659       min_sinc_table_length = st->filt_len*st->oversample+8;
660    }
661    if (st->sinc_table_length < min_sinc_table_length)
662    {
663       spx_word16_t *sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,min_sinc_table_length*sizeof(spx_word16_t));
664       if (!sinc_table)
665          goto fail;
666
667       st->sinc_table = sinc_table;
668       st->sinc_table_length = min_sinc_table_length;
669    }
670    if (use_direct)
671    {
672       spx_uint32_t i;
673       for (i=0;i<st->den_rate;i++)
674       {
675          spx_int32_t j;
676          for (j=0;j<st->filt_len;j++)
677          {
678             st->sinc_table[i*st->filt_len+j] = sinc(st->cutoff,((j-(spx_int32_t)st->filt_len/2+1)-((float)i)/st->den_rate), st->filt_len, quality_map[st->quality].window_func);
679          }
680       }
681 #ifdef FIXED_POINT
682       st->resampler_ptr = resampler_basic_direct_single;
683 #else
684       if (st->quality>8)
685          st->resampler_ptr = resampler_basic_direct_double;
686       else
687          st->resampler_ptr = resampler_basic_direct_single;
688 #endif
689       /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff);*/
690    } else {
691       spx_int32_t i;
692       for (i=-4;i<(spx_int32_t)(st->oversample*st->filt_len+4);i++)
693          st->sinc_table[i+4] = sinc(st->cutoff,(i/(float)st->oversample - st->filt_len/2), st->filt_len, quality_map[st->quality].window_func);
694 #ifdef FIXED_POINT
695       st->resampler_ptr = resampler_basic_interpolate_single;
696 #else
697       if (st->quality>8)
698          st->resampler_ptr = resampler_basic_interpolate_double;
699       else
700          st->resampler_ptr = resampler_basic_interpolate_single;
701 #endif
702       /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff);*/
703    }
704
705    /* Here's the place where we update the filter memory to take into account
706       the change in filter length. It's probably the messiest part of the code
707       due to handling of lots of corner cases. */
708
709    /* Adding buffer_size to filt_len won't overflow here because filt_len
710       could be multiplied by sizeof(spx_word16_t) above. */
711    min_alloc_size = st->filt_len-1 + st->buffer_size;
712    if (min_alloc_size > st->mem_alloc_size)
713    {
714       spx_word16_t *mem;
715       if (INT_MAX/sizeof(spx_word16_t)/st->nb_channels < min_alloc_size)
716           goto fail;
717       else if (!(mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*min_alloc_size * sizeof(*mem))))
718           goto fail;
719
720       st->mem = mem;
721       st->mem_alloc_size = min_alloc_size;
722    }
723    if (!st->started)
724    {
725       spx_uint32_t i;
726       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
727          st->mem[i] = 0;
728       /*speex_warning("reinit filter");*/
729    } else if (st->filt_len > old_length)
730    {
731       spx_uint32_t i;
732       /* Increase the filter length */
733       /*speex_warning("increase filter size");*/
734       for (i=st->nb_channels;i--;)
735       {
736          spx_uint32_t j;
737          spx_uint32_t olen = old_length;
738          /*if (st->magic_samples[i])*/
739          {
740             /* Try and remove the magic samples as if nothing had happened */
741
742             /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
743             olen = old_length + 2*st->magic_samples[i];
744             for (j=old_length-1+st->magic_samples[i];j--;)
745                st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]] = st->mem[i*old_alloc_size+j];
746             for (j=0;j<st->magic_samples[i];j++)
747                st->mem[i*st->mem_alloc_size+j] = 0;
748             st->magic_samples[i] = 0;
749          }
750          if (st->filt_len > olen)
751          {
752             /* If the new filter length is still bigger than the "augmented" length */
753             /* Copy data going backward */
754             for (j=0;j<olen-1;j++)
755                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = st->mem[i*st->mem_alloc_size+(olen-2-j)];
756             /* Then put zeros for lack of anything better */
757             for (;j<st->filt_len-1;j++)
758                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = 0;
759             /* Adjust last_sample */
760             st->last_sample[i] += (st->filt_len - olen)/2;
761          } else {
762             /* Put back some of the magic! */
763             st->magic_samples[i] = (olen - st->filt_len)/2;
764             for (j=0;j<st->filt_len-1+st->magic_samples[i];j++)
765                st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
766          }
767       }
768    } else if (st->filt_len < old_length)
769    {
770       spx_uint32_t i;
771       /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"
772          samples so they can be used directly as input the next time(s) */
773       for (i=0;i<st->nb_channels;i++)
774       {
775          spx_uint32_t j;
776          spx_uint32_t old_magic = st->magic_samples[i];
777          st->magic_samples[i] = (old_length - st->filt_len)/2;
778          /* We must copy some of the memory that's no longer used */
779          /* Copy data going backward */
780          for (j=0;j<st->filt_len-1+st->magic_samples[i]+old_magic;j++)
781             st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
782          st->magic_samples[i] += old_magic;
783       }
784    }
785    return RESAMPLER_ERR_SUCCESS;
786
787 fail:
788    st->resampler_ptr = resampler_basic_zero;
789    /* st->mem may still contain consumed input samples for the filter.
790       Restore filt_len so that filt_len - 1 still points to the position after
791       the last of these samples. */
792    st->filt_len = old_length;
793    return RESAMPLER_ERR_ALLOC_FAILED;
794 }
795
796 EXPORT SpeexResamplerState *speex_resampler_init(spx_uint32_t nb_channels, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
797 {
798    return speex_resampler_init_frac(nb_channels, in_rate, out_rate, in_rate, out_rate, quality, err);
799 }
800
801 EXPORT SpeexResamplerState *speex_resampler_init_frac(spx_uint32_t nb_channels, spx_uint32_t ratio_num, spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
802 {
803    spx_uint32_t i;
804    SpeexResamplerState *st;
805    int filter_err;
806
807    if (quality > 10 || quality < 0)
808    {
809       if (err)
810          *err = RESAMPLER_ERR_INVALID_ARG;
811       return NULL;
812    }
813    st = (SpeexResamplerState *)speex_alloc(sizeof(SpeexResamplerState));
814    if (!st)
815    {
816       if (err)
817          *err = RESAMPLER_ERR_ALLOC_FAILED;
818       return NULL;
819    }
820    st->initialised = 0;
821    st->started = 0;
822    st->in_rate = 0;
823    st->out_rate = 0;
824    st->num_rate = 0;
825    st->den_rate = 0;
826    st->quality = -1;
827    st->sinc_table_length = 0;
828    st->mem_alloc_size = 0;
829    st->filt_len = 0;
830    st->mem = 0;
831    st->resampler_ptr = 0;
832
833    st->cutoff = 1.f;
834    st->nb_channels = nb_channels;
835    st->in_stride = 1;
836    st->out_stride = 1;
837
838    st->buffer_size = 160;
839
840    /* Per channel data */
841    if (!(st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(spx_int32_t))))
842       goto fail;
843    if (!(st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t))))
844       goto fail;
845    if (!(st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t))))
846       goto fail;
847
848    speex_resampler_set_quality(st, quality);
849    speex_resampler_set_rate_frac(st, ratio_num, ratio_den, in_rate, out_rate);
850
851    filter_err = update_filter(st);
852    if (filter_err == RESAMPLER_ERR_SUCCESS)
853    {
854       st->initialised = 1;
855    } else {
856       speex_resampler_destroy(st);
857       st = NULL;
858    }
859    if (err)
860       *err = filter_err;
861
862    return st;
863
864 fail:
865    if (err)
866       *err = RESAMPLER_ERR_ALLOC_FAILED;
867    speex_resampler_destroy(st);
868    return NULL;
869 }
870
871 EXPORT void speex_resampler_destroy(SpeexResamplerState *st)
872 {
873    speex_free(st->mem);
874    speex_free(st->sinc_table);
875    speex_free(st->last_sample);
876    speex_free(st->magic_samples);
877    speex_free(st->samp_frac_num);
878    speex_free(st);
879 }
880
881 static int speex_resampler_process_native(SpeexResamplerState *st, spx_uint32_t channel_index, spx_uint32_t *in_len, spx_word16_t *out, spx_uint32_t *out_len)
882 {
883    int j=0;
884    const int N = st->filt_len;
885    int out_sample = 0;
886    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
887    spx_uint32_t ilen;
888
889    st->started = 1;
890
891    /* Call the right resampler through the function ptr */
892    out_sample = st->resampler_ptr(st, channel_index, mem, in_len, out, out_len);
893
894    if (st->last_sample[channel_index] < (spx_int32_t)*in_len)
895       *in_len = st->last_sample[channel_index];
896    *out_len = out_sample;
897    st->last_sample[channel_index] -= *in_len;
898
899    ilen = *in_len;
900
901    for(j=0;j<N-1;++j)
902      mem[j] = mem[j+ilen];
903
904    return RESAMPLER_ERR_SUCCESS;
905 }
906
907 static int speex_resampler_magic(SpeexResamplerState *st, spx_uint32_t channel_index, spx_word16_t **out, spx_uint32_t out_len) {
908    spx_uint32_t tmp_in_len = st->magic_samples[channel_index];
909    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
910    const int N = st->filt_len;
911
912    speex_resampler_process_native(st, channel_index, &tmp_in_len, *out, &out_len);
913
914    st->magic_samples[channel_index] -= tmp_in_len;
915
916    /* If we couldn't process all "magic" input samples, save the rest for next time */
917    if (st->magic_samples[channel_index])
918    {
919       spx_uint32_t i;
920       for (i=0;i<st->magic_samples[channel_index];i++)
921          mem[N-1+i]=mem[N-1+i+tmp_in_len];
922    }
923    *out += out_len*st->out_stride;
924    return out_len;
925 }
926
927 #ifdef FIXED_POINT
928 EXPORT int speex_resampler_process_int(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_int16_t *in, spx_uint32_t *in_len, spx_int16_t *out, spx_uint32_t *out_len)
929 #else
930 EXPORT int speex_resampler_process_float(SpeexResamplerState *st, spx_uint32_t channel_index, const float *in, spx_uint32_t *in_len, float *out, spx_uint32_t *out_len)
931 #endif
932 {
933    int j;
934    spx_uint32_t ilen = *in_len;
935    spx_uint32_t olen = *out_len;
936    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
937    const int filt_offs = st->filt_len - 1;
938    const spx_uint32_t xlen = st->mem_alloc_size - filt_offs;
939    const int istride = st->in_stride;
940
941    if (st->magic_samples[channel_index])
942       olen -= speex_resampler_magic(st, channel_index, &out, olen);
943    if (! st->magic_samples[channel_index]) {
944       while (ilen && olen) {
945         spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
946         spx_uint32_t ochunk = olen;
947
948         if (in) {
949            for(j=0;j<ichunk;++j)
950               x[j+filt_offs]=in[j*istride];
951         } else {
952           for(j=0;j<ichunk;++j)
953             x[j+filt_offs]=0;
954         }
955         speex_resampler_process_native(st, channel_index, &ichunk, out, &ochunk);
956         ilen -= ichunk;
957         olen -= ochunk;
958         out += ochunk * st->out_stride;
959         if (in)
960            in += ichunk * istride;
961       }
962    }
963    *in_len -= ilen;
964    *out_len -= olen;
965    return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
966 }
967
968 #ifdef FIXED_POINT
969 EXPORT int speex_resampler_process_float(SpeexResamplerState *st, spx_uint32_t channel_index, const float *in, spx_uint32_t *in_len, float *out, spx_uint32_t *out_len)
970 #else
971 EXPORT int speex_resampler_process_int(SpeexResamplerState *st, spx_uint32_t channel_index, const spx_int16_t *in, spx_uint32_t *in_len, spx_int16_t *out, spx_uint32_t *out_len)
972 #endif
973 {
974    int j;
975    const int istride_save = st->in_stride;
976    const int ostride_save = st->out_stride;
977    spx_uint32_t ilen = *in_len;
978    spx_uint32_t olen = *out_len;
979    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
980    const spx_uint32_t xlen = st->mem_alloc_size - (st->filt_len - 1);
981 #ifdef VAR_ARRAYS
982    const unsigned int ylen = (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;
983    VARDECL(spx_word16_t *ystack);
984    ALLOC(ystack, ylen, spx_word16_t);
985 #else
986    const unsigned int ylen = FIXED_STACK_ALLOC;
987    spx_word16_t ystack[FIXED_STACK_ALLOC];
988 #endif
989
990    st->out_stride = 1;
991
992    while (ilen && olen) {
993      spx_word16_t *y = ystack;
994      spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
995      spx_uint32_t ochunk = (olen > ylen) ? ylen : olen;
996      spx_uint32_t omagic = 0;
997
998      if (st->magic_samples[channel_index]) {
999        omagic = speex_resampler_magic(st, channel_index, &y, ochunk);
1000        ochunk -= omagic;
1001        olen -= omagic;
1002      }
1003      if (! st->magic_samples[channel_index]) {
1004        if (in) {
1005          for(j=0;j<ichunk;++j)
1006 #ifdef FIXED_POINT
1007            x[j+st->filt_len-1]=WORD2INT(in[j*istride_save]);
1008 #else
1009            x[j+st->filt_len-1]=in[j*istride_save];
1010 #endif
1011        } else {
1012          for(j=0;j<ichunk;++j)
1013            x[j+st->filt_len-1]=0;
1014        }
1015
1016        speex_resampler_process_native(st, channel_index, &ichunk, y, &ochunk);
1017      } else {
1018        ichunk = 0;
1019        ochunk = 0;
1020      }
1021
1022      for (j=0;j<ochunk+omagic;++j)
1023 #ifdef FIXED_POINT
1024         out[j*ostride_save] = ystack[j];
1025 #else
1026         out[j*ostride_save] = WORD2INT(ystack[j]);
1027 #endif
1028
1029      ilen -= ichunk;
1030      olen -= ochunk;
1031      out += (ochunk+omagic) * ostride_save;
1032      if (in)
1033        in += ichunk * istride_save;
1034    }
1035    st->out_stride = ostride_save;
1036    *in_len -= ilen;
1037    *out_len -= olen;
1038
1039    return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
1040 }
1041
1042 EXPORT int speex_resampler_process_interleaved_float(SpeexResamplerState *st, const float *in, spx_uint32_t *in_len, float *out, spx_uint32_t *out_len)
1043 {
1044    spx_uint32_t i;
1045    int istride_save, ostride_save;
1046    spx_uint32_t bak_out_len = *out_len;
1047    spx_uint32_t bak_in_len = *in_len;
1048    istride_save = st->in_stride;
1049    ostride_save = st->out_stride;
1050    st->in_stride = st->out_stride = st->nb_channels;
1051    for (i=0;i<st->nb_channels;i++)
1052    {
1053       *out_len = bak_out_len;
1054       *in_len = bak_in_len;
1055       if (in != NULL)
1056          speex_resampler_process_float(st, i, in+i, in_len, out+i, out_len);
1057       else
1058          speex_resampler_process_float(st, i, NULL, in_len, out+i, out_len);
1059    }
1060    st->in_stride = istride_save;
1061    st->out_stride = ostride_save;
1062    return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
1063 }
1064
1065 EXPORT int speex_resampler_process_interleaved_int(SpeexResamplerState *st, const spx_int16_t *in, spx_uint32_t *in_len, spx_int16_t *out, spx_uint32_t *out_len)
1066 {
1067    spx_uint32_t i;
1068    int istride_save, ostride_save;
1069    spx_uint32_t bak_out_len = *out_len;
1070    spx_uint32_t bak_in_len = *in_len;
1071    istride_save = st->in_stride;
1072    ostride_save = st->out_stride;
1073    st->in_stride = st->out_stride = st->nb_channels;
1074    for (i=0;i<st->nb_channels;i++)
1075    {
1076       *out_len = bak_out_len;
1077       *in_len = bak_in_len;
1078       if (in != NULL)
1079          speex_resampler_process_int(st, i, in+i, in_len, out+i, out_len);
1080       else
1081          speex_resampler_process_int(st, i, NULL, in_len, out+i, out_len);
1082    }
1083    st->in_stride = istride_save;
1084    st->out_stride = ostride_save;
1085    return st->resampler_ptr == resampler_basic_zero ? RESAMPLER_ERR_ALLOC_FAILED : RESAMPLER_ERR_SUCCESS;
1086 }
1087
1088 EXPORT int speex_resampler_set_rate(SpeexResamplerState *st, spx_uint32_t in_rate, spx_uint32_t out_rate)
1089 {
1090    return speex_resampler_set_rate_frac(st, in_rate, out_rate, in_rate, out_rate);
1091 }
1092
1093 EXPORT void speex_resampler_get_rate(SpeexResamplerState *st, spx_uint32_t *in_rate, spx_uint32_t *out_rate)
1094 {
1095    *in_rate = st->in_rate;
1096    *out_rate = st->out_rate;
1097 }
1098
1099 static inline spx_uint32_t _gcd(spx_uint32_t a, spx_uint32_t b)
1100 {
1101    while (b != 0)
1102    {
1103       spx_uint32_t temp = a;
1104
1105       a = b;
1106       b = temp % b;
1107    }
1108    return a;
1109 }
1110
1111 EXPORT int speex_resampler_set_rate_frac(SpeexResamplerState *st, spx_uint32_t ratio_num, spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate)
1112 {
1113    spx_uint32_t fact;
1114    spx_uint32_t old_den;
1115    spx_uint32_t i;
1116    if (st->in_rate == in_rate && st->out_rate == out_rate && st->num_rate == ratio_num && st->den_rate == ratio_den)
1117       return RESAMPLER_ERR_SUCCESS;
1118
1119    old_den = st->den_rate;
1120    st->in_rate = in_rate;
1121    st->out_rate = out_rate;
1122    st->num_rate = ratio_num;
1123    st->den_rate = ratio_den;
1124
1125    fact = _gcd (st->num_rate, st->den_rate);
1126
1127    st->num_rate /= fact;
1128    st->den_rate /= fact;
1129
1130    if (old_den > 0)
1131    {
1132       for (i=0;i<st->nb_channels;i++)
1133       {
1134          if (_muldiv(&st->samp_frac_num[i],st->samp_frac_num[i],st->den_rate,old_den) != RESAMPLER_ERR_SUCCESS)
1135             return RESAMPLER_ERR_OVERFLOW;
1136          /* Safety net */
1137          if (st->samp_frac_num[i] >= st->den_rate)
1138             st->samp_frac_num[i] = st->den_rate-1;
1139       }
1140    }
1141
1142    if (st->initialised)
1143       return update_filter(st);
1144    return RESAMPLER_ERR_SUCCESS;
1145 }
1146
1147 EXPORT void speex_resampler_get_ratio(SpeexResamplerState *st, spx_uint32_t *ratio_num, spx_uint32_t *ratio_den)
1148 {
1149    *ratio_num = st->num_rate;
1150    *ratio_den = st->den_rate;
1151 }
1152
1153 EXPORT int speex_resampler_set_quality(SpeexResamplerState *st, int quality)
1154 {
1155    if (quality > 10 || quality < 0)
1156       return RESAMPLER_ERR_INVALID_ARG;
1157    if (st->quality == quality)
1158       return RESAMPLER_ERR_SUCCESS;
1159    st->quality = quality;
1160    if (st->initialised)
1161       return update_filter(st);
1162    return RESAMPLER_ERR_SUCCESS;
1163 }
1164
1165 EXPORT void speex_resampler_get_quality(SpeexResamplerState *st, int *quality)
1166 {
1167    *quality = st->quality;
1168 }
1169
1170 EXPORT void speex_resampler_set_input_stride(SpeexResamplerState *st, spx_uint32_t stride)
1171 {
1172    st->in_stride = stride;
1173 }
1174
1175 EXPORT void speex_resampler_get_input_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1176 {
1177    *stride = st->in_stride;
1178 }
1179
1180 EXPORT void speex_resampler_set_output_stride(SpeexResamplerState *st, spx_uint32_t stride)
1181 {
1182    st->out_stride = stride;
1183 }
1184
1185 EXPORT void speex_resampler_get_output_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1186 {
1187    *stride = st->out_stride;
1188 }
1189
1190 EXPORT int speex_resampler_get_input_latency(SpeexResamplerState *st)
1191 {
1192   return st->filt_len / 2;
1193 }
1194
1195 EXPORT int speex_resampler_get_output_latency(SpeexResamplerState *st)
1196 {
1197   return ((st->filt_len / 2) * st->den_rate + (st->num_rate >> 1)) / st->num_rate;
1198 }
1199
1200 EXPORT int speex_resampler_skip_zeros(SpeexResamplerState *st)
1201 {
1202    spx_uint32_t i;
1203    for (i=0;i<st->nb_channels;i++)
1204       st->last_sample[i] = st->filt_len/2;
1205    return RESAMPLER_ERR_SUCCESS;
1206 }
1207
1208 EXPORT int speex_resampler_reset_mem(SpeexResamplerState *st)
1209 {
1210    spx_uint32_t i;
1211    for (i=0;i<st->nb_channels;i++)
1212    {
1213       st->last_sample[i] = 0;
1214       st->magic_samples[i] = 0;
1215       st->samp_frac_num[i] = 0;
1216    }
1217    for (i=0;i<st->nb_channels*(st->filt_len-1);i++)
1218       st->mem[i] = 0;
1219    return RESAMPLER_ERR_SUCCESS;
1220 }
1221
1222 EXPORT const char *speex_resampler_strerror(int err)
1223 {
1224    switch (err)
1225    {
1226       case RESAMPLER_ERR_SUCCESS:
1227          return "Success.";
1228       case RESAMPLER_ERR_ALLOC_FAILED:
1229          return "Memory allocation failed.";
1230       case RESAMPLER_ERR_BAD_STATE:
1231          return "Bad resampler state.";
1232       case RESAMPLER_ERR_INVALID_ARG:
1233          return "Invalid argument.";
1234       case RESAMPLER_ERR_PTR_OVERLAP:
1235          return "Input and output buffers overlap.";
1236       default:
1237          return "Unknown error. Bad error code or strange version mismatch.";
1238    }
1239 }