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