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