resample: Calculate full sinc table (e.g. no sinc interpolation)
[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 #ifdef RESAMPLE_FULL_SINC_TABLE
589    if (1)
590 #else
591    if (st->filt_len*st->den_rate <= st->filt_len*st->oversample+8)
592 #endif
593    {
594       spx_uint32_t i;
595       if (!st->sinc_table)
596          st->sinc_table = (spx_word16_t *)speex_alloc(st->filt_len*st->den_rate*sizeof(spx_word16_t));
597       else if (st->sinc_table_length < st->filt_len*st->den_rate)
598       {
599          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,st->filt_len*st->den_rate*sizeof(spx_word16_t));
600          st->sinc_table_length = st->filt_len*st->den_rate;
601       }
602       for (i=0;i<st->den_rate;i++)
603       {
604          spx_int32_t j;
605          for (j=0;j<st->filt_len;j++)
606          {
607             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);
608          }
609       }
610 #ifdef FIXED_POINT
611       st->resampler_ptr = resampler_basic_direct_single;
612 #else
613       if (st->quality>8)
614          st->resampler_ptr = resampler_basic_direct_double;
615       else
616          st->resampler_ptr = resampler_basic_direct_single;
617 #endif
618       /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff);*/
619    } else {
620       spx_int32_t i;
621       if (!st->sinc_table)
622          st->sinc_table = (spx_word16_t *)speex_alloc((st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
623       else if (st->sinc_table_length < st->filt_len*st->oversample+8)
624       {
625          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,(st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
626          st->sinc_table_length = st->filt_len*st->oversample+8;
627       }
628       for (i=-4;i<(spx_int32_t)(st->oversample*st->filt_len+4);i++)
629          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);
630 #ifdef FIXED_POINT
631       st->resampler_ptr = resampler_basic_interpolate_single;
632 #else
633       if (st->quality>8)
634          st->resampler_ptr = resampler_basic_interpolate_double;
635       else
636          st->resampler_ptr = resampler_basic_interpolate_single;
637 #endif
638       /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff);*/
639    }
640    st->int_advance = st->num_rate/st->den_rate;
641    st->frac_advance = st->num_rate%st->den_rate;
642
643    
644    /* Here's the place where we update the filter memory to take into account
645       the change in filter length. It's probably the messiest part of the code
646       due to handling of lots of corner cases. */
647    if (!st->mem)
648    {
649       spx_uint32_t i;
650       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
651       st->mem = (spx_word16_t*)speex_alloc(st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
652       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
653          st->mem[i] = 0;
654       /*speex_warning("init filter");*/
655    } else if (!st->started)
656    {
657       spx_uint32_t i;
658       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
659       st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
660       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
661          st->mem[i] = 0;
662       /*speex_warning("reinit filter");*/
663    } else if (st->filt_len > old_length)
664    {
665       spx_int32_t i;
666       /* Increase the filter length */
667       /*speex_warning("increase filter size");*/
668       int old_alloc_size = st->mem_alloc_size;
669       if ((st->filt_len-1 + st->buffer_size) > st->mem_alloc_size)
670       {
671          st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
672          st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
673       }
674       for (i=st->nb_channels-1;i>=0;i--)
675       {
676          spx_int32_t j;
677          spx_uint32_t olen = old_length;
678          /*if (st->magic_samples[i])*/
679          {
680             /* Try and remove the magic samples as if nothing had happened */
681             
682             /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
683             olen = old_length + 2*st->magic_samples[i];
684             for (j=old_length-2+st->magic_samples[i];j>=0;j--)
685                st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]] = st->mem[i*old_alloc_size+j];
686             for (j=0;j<st->magic_samples[i];j++)
687                st->mem[i*st->mem_alloc_size+j] = 0;
688             st->magic_samples[i] = 0;
689          }
690          if (st->filt_len > olen)
691          {
692             /* If the new filter length is still bigger than the "augmented" length */
693             /* Copy data going backward */
694             for (j=0;j<olen-1;j++)
695                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = st->mem[i*st->mem_alloc_size+(olen-2-j)];
696             /* Then put zeros for lack of anything better */
697             for (;j<st->filt_len-1;j++)
698                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = 0;
699             /* Adjust last_sample */
700             st->last_sample[i] += (st->filt_len - olen)/2;
701          } else {
702             /* Put back some of the magic! */
703             st->magic_samples[i] = (olen - st->filt_len)/2;
704             for (j=0;j<st->filt_len-1+st->magic_samples[i];j++)
705                st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
706          }
707       }
708    } else if (st->filt_len < old_length)
709    {
710       spx_uint32_t i;
711       /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"
712          samples so they can be used directly as input the next time(s) */
713       for (i=0;i<st->nb_channels;i++)
714       {
715          spx_uint32_t j;
716          spx_uint32_t old_magic = st->magic_samples[i];
717          st->magic_samples[i] = (old_length - st->filt_len)/2;
718          /* We must copy some of the memory that's no longer used */
719          /* Copy data going backward */
720          for (j=0;j<st->filt_len-1+st->magic_samples[i]+old_magic;j++)
721             st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
722          st->magic_samples[i] += old_magic;
723       }
724    }
725
726 }
727
728 EXPORT SpeexResamplerState *speex_resampler_init(spx_uint32_t nb_channels, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
729 {
730    return speex_resampler_init_frac(nb_channels, in_rate, out_rate, in_rate, out_rate, quality, err);
731 }
732
733 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)
734 {
735    spx_uint32_t i;
736    SpeexResamplerState *st;
737    if (quality > 10 || quality < 0)
738    {
739       if (err)
740          *err = RESAMPLER_ERR_INVALID_ARG;
741       return NULL;
742    }
743    st = (SpeexResamplerState *)speex_alloc(sizeof(SpeexResamplerState));
744    st->initialised = 0;
745    st->started = 0;
746    st->in_rate = 0;
747    st->out_rate = 0;
748    st->num_rate = 0;
749    st->den_rate = 0;
750    st->quality = -1;
751    st->sinc_table_length = 0;
752    st->mem_alloc_size = 0;
753    st->filt_len = 0;
754    st->mem = 0;
755    st->resampler_ptr = 0;
756          
757    st->cutoff = 1.f;
758    st->nb_channels = nb_channels;
759    st->in_stride = 1;
760    st->out_stride = 1;
761    
762    st->buffer_size = 160;
763    
764    /* Per channel data */
765    st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(spx_int32_t));
766    st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t));
767    st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(spx_uint32_t));
768    for (i=0;i<nb_channels;i++)
769    {
770       st->last_sample[i] = 0;
771       st->magic_samples[i] = 0;
772       st->samp_frac_num[i] = 0;
773    }
774
775    speex_resampler_set_quality(st, quality);
776    speex_resampler_set_rate_frac(st, ratio_num, ratio_den, in_rate, out_rate);
777
778    
779    update_filter(st);
780    
781    st->initialised = 1;
782    if (err)
783       *err = RESAMPLER_ERR_SUCCESS;
784
785    return st;
786 }
787
788 EXPORT void speex_resampler_destroy(SpeexResamplerState *st)
789 {
790    speex_free(st->mem);
791    speex_free(st->sinc_table);
792    speex_free(st->last_sample);
793    speex_free(st->magic_samples);
794    speex_free(st->samp_frac_num);
795    speex_free(st);
796 }
797
798 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)
799 {
800    int j=0;
801    const int N = st->filt_len;
802    int out_sample = 0;
803    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
804    spx_uint32_t ilen;
805    
806    st->started = 1;
807    
808    /* Call the right resampler through the function ptr */
809    out_sample = st->resampler_ptr(st, channel_index, mem, in_len, out, out_len);
810    
811    if (st->last_sample[channel_index] < (spx_int32_t)*in_len)
812       *in_len = st->last_sample[channel_index];
813    *out_len = out_sample;
814    st->last_sample[channel_index] -= *in_len;
815    
816    ilen = *in_len;
817
818    for(j=0;j<N-1;++j)
819      mem[j] = mem[j+ilen];
820
821    return RESAMPLER_ERR_SUCCESS;
822 }
823
824 static int speex_resampler_magic(SpeexResamplerState *st, spx_uint32_t channel_index, spx_word16_t **out, spx_uint32_t out_len) {
825    spx_uint32_t tmp_in_len = st->magic_samples[channel_index];
826    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
827    const int N = st->filt_len;
828    
829    speex_resampler_process_native(st, channel_index, &tmp_in_len, *out, &out_len);
830
831    st->magic_samples[channel_index] -= tmp_in_len;
832    
833    /* If we couldn't process all "magic" input samples, save the rest for next time */
834    if (st->magic_samples[channel_index])
835    {
836       spx_uint32_t i;
837       for (i=0;i<st->magic_samples[channel_index];i++)
838          mem[N-1+i]=mem[N-1+i+tmp_in_len];
839    }
840    *out += out_len*st->out_stride;
841    return out_len;
842 }
843
844 #ifdef FIXED_POINT
845 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)
846 #else
847 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)
848 #endif
849 {
850    int j;
851    spx_uint32_t ilen = *in_len;
852    spx_uint32_t olen = *out_len;
853    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
854    const int filt_offs = st->filt_len - 1;
855    const spx_uint32_t xlen = st->mem_alloc_size - filt_offs;
856    const int istride = st->in_stride;
857
858    if (st->magic_samples[channel_index]) 
859       olen -= speex_resampler_magic(st, channel_index, &out, olen);
860    if (! st->magic_samples[channel_index]) {
861       while (ilen && olen) {
862         spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
863         spx_uint32_t ochunk = olen;
864  
865         if (in) {
866            for(j=0;j<ichunk;++j)
867               x[j+filt_offs]=in[j*istride];
868         } else {
869           for(j=0;j<ichunk;++j)
870             x[j+filt_offs]=0;
871         }
872         speex_resampler_process_native(st, channel_index, &ichunk, out, &ochunk);
873         ilen -= ichunk;
874         olen -= ochunk;
875         out += ochunk * st->out_stride;
876         if (in)
877            in += ichunk * istride;
878       }
879    }
880    *in_len -= ilen;
881    *out_len -= olen;
882    return RESAMPLER_ERR_SUCCESS;
883 }
884
885 #ifdef FIXED_POINT
886 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)
887 #else
888 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)
889 #endif
890 {
891    int j;
892    const int istride_save = st->in_stride;
893    const int ostride_save = st->out_stride;
894    spx_uint32_t ilen = *in_len;
895    spx_uint32_t olen = *out_len;
896    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
897    const spx_uint32_t xlen = st->mem_alloc_size - (st->filt_len - 1);
898 #ifdef VAR_ARRAYS
899    const unsigned int ylen = (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;
900    VARDECL(spx_word16_t *ystack);
901    ALLOC(ystack, ylen, spx_word16_t);
902 #else
903    const unsigned int ylen = FIXED_STACK_ALLOC;
904    spx_word16_t ystack[FIXED_STACK_ALLOC];
905 #endif
906
907    st->out_stride = 1;
908    
909    while (ilen && olen) {
910      spx_word16_t *y = ystack;
911      spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
912      spx_uint32_t ochunk = (olen > ylen) ? ylen : olen;
913      spx_uint32_t omagic = 0;
914
915      if (st->magic_samples[channel_index]) {
916        omagic = speex_resampler_magic(st, channel_index, &y, ochunk);
917        ochunk -= omagic;
918        olen -= omagic;
919      }
920      if (! st->magic_samples[channel_index]) {
921        if (in) {
922          for(j=0;j<ichunk;++j)
923 #ifdef FIXED_POINT
924            x[j+st->filt_len-1]=WORD2INT(in[j*istride_save]);
925 #else
926            x[j+st->filt_len-1]=in[j*istride_save];
927 #endif
928        } else {
929          for(j=0;j<ichunk;++j)
930            x[j+st->filt_len-1]=0;
931        }
932
933        speex_resampler_process_native(st, channel_index, &ichunk, y, &ochunk);
934      } else {
935        ichunk = 0;
936        ochunk = 0;
937      }
938
939      for (j=0;j<ochunk+omagic;++j)
940 #ifdef FIXED_POINT
941         out[j*ostride_save] = ystack[j];
942 #else
943         out[j*ostride_save] = WORD2INT(ystack[j]);
944 #endif
945      
946      ilen -= ichunk;
947      olen -= ochunk;
948      out += (ochunk+omagic) * ostride_save;
949      if (in)
950        in += ichunk * istride_save;
951    }
952    st->out_stride = ostride_save;
953    *in_len -= ilen;
954    *out_len -= olen;
955
956    return RESAMPLER_ERR_SUCCESS;
957 }
958
959 EXPORT int speex_resampler_process_interleaved_float(SpeexResamplerState *st, const float *in, spx_uint32_t *in_len, float *out, spx_uint32_t *out_len)
960 {
961    spx_uint32_t i;
962    int istride_save, ostride_save;
963    spx_uint32_t bak_out_len = *out_len;
964    spx_uint32_t bak_in_len = *in_len;
965    istride_save = st->in_stride;
966    ostride_save = st->out_stride;
967    st->in_stride = st->out_stride = st->nb_channels;
968    for (i=0;i<st->nb_channels;i++)
969    {
970       *out_len = bak_out_len;
971       *in_len = bak_in_len;
972       if (in != NULL)
973          speex_resampler_process_float(st, i, in+i, in_len, out+i, out_len);
974       else
975          speex_resampler_process_float(st, i, NULL, in_len, out+i, out_len);
976    }
977    st->in_stride = istride_save;
978    st->out_stride = ostride_save;
979    return RESAMPLER_ERR_SUCCESS;
980 }
981                
982 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)
983 {
984    spx_uint32_t i;
985    int istride_save, ostride_save;
986    spx_uint32_t bak_out_len = *out_len;
987    spx_uint32_t bak_in_len = *in_len;
988    istride_save = st->in_stride;
989    ostride_save = st->out_stride;
990    st->in_stride = st->out_stride = st->nb_channels;
991    for (i=0;i<st->nb_channels;i++)
992    {
993       *out_len = bak_out_len;
994       *in_len = bak_in_len;
995       if (in != NULL)
996          speex_resampler_process_int(st, i, in+i, in_len, out+i, out_len);
997       else
998          speex_resampler_process_int(st, i, NULL, in_len, out+i, out_len);
999    }
1000    st->in_stride = istride_save;
1001    st->out_stride = ostride_save;
1002    return RESAMPLER_ERR_SUCCESS;
1003 }
1004
1005 EXPORT int speex_resampler_set_rate(SpeexResamplerState *st, spx_uint32_t in_rate, spx_uint32_t out_rate)
1006 {
1007    return speex_resampler_set_rate_frac(st, in_rate, out_rate, in_rate, out_rate);
1008 }
1009
1010 EXPORT void speex_resampler_get_rate(SpeexResamplerState *st, spx_uint32_t *in_rate, spx_uint32_t *out_rate)
1011 {
1012    *in_rate = st->in_rate;
1013    *out_rate = st->out_rate;
1014 }
1015
1016 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)
1017 {
1018    spx_uint32_t fact;
1019    spx_uint32_t old_den;
1020    spx_uint32_t i;
1021    if (st->in_rate == in_rate && st->out_rate == out_rate && st->num_rate == ratio_num && st->den_rate == ratio_den)
1022       return RESAMPLER_ERR_SUCCESS;
1023    
1024    old_den = st->den_rate;
1025    st->in_rate = in_rate;
1026    st->out_rate = out_rate;
1027    st->num_rate = ratio_num;
1028    st->den_rate = ratio_den;
1029    /* FIXME: This is terribly inefficient, but who cares (at least for now)? */
1030    for (fact=2;fact<=IMIN(st->num_rate, st->den_rate);fact++)
1031    {
1032       while ((st->num_rate % fact == 0) && (st->den_rate % fact == 0))
1033       {
1034          st->num_rate /= fact;
1035          st->den_rate /= fact;
1036       }
1037    }
1038       
1039    if (old_den > 0)
1040    {
1041       for (i=0;i<st->nb_channels;i++)
1042       {
1043          st->samp_frac_num[i]=st->samp_frac_num[i]*st->den_rate/old_den;
1044          /* Safety net */
1045          if (st->samp_frac_num[i] >= st->den_rate)
1046             st->samp_frac_num[i] = st->den_rate-1;
1047       }
1048    }
1049    
1050    if (st->initialised)
1051       update_filter(st);
1052    return RESAMPLER_ERR_SUCCESS;
1053 }
1054
1055 EXPORT void speex_resampler_get_ratio(SpeexResamplerState *st, spx_uint32_t *ratio_num, spx_uint32_t *ratio_den)
1056 {
1057    *ratio_num = st->num_rate;
1058    *ratio_den = st->den_rate;
1059 }
1060
1061 EXPORT int speex_resampler_set_quality(SpeexResamplerState *st, int quality)
1062 {
1063    if (quality > 10 || quality < 0)
1064       return RESAMPLER_ERR_INVALID_ARG;
1065    if (st->quality == quality)
1066       return RESAMPLER_ERR_SUCCESS;
1067    st->quality = quality;
1068    if (st->initialised)
1069       update_filter(st);
1070    return RESAMPLER_ERR_SUCCESS;
1071 }
1072
1073 EXPORT void speex_resampler_get_quality(SpeexResamplerState *st, int *quality)
1074 {
1075    *quality = st->quality;
1076 }
1077
1078 EXPORT void speex_resampler_set_input_stride(SpeexResamplerState *st, spx_uint32_t stride)
1079 {
1080    st->in_stride = stride;
1081 }
1082
1083 EXPORT void speex_resampler_get_input_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1084 {
1085    *stride = st->in_stride;
1086 }
1087
1088 EXPORT void speex_resampler_set_output_stride(SpeexResamplerState *st, spx_uint32_t stride)
1089 {
1090    st->out_stride = stride;
1091 }
1092
1093 EXPORT void speex_resampler_get_output_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1094 {
1095    *stride = st->out_stride;
1096 }
1097
1098 EXPORT int speex_resampler_get_input_latency(SpeexResamplerState *st)
1099 {
1100   return st->filt_len / 2;
1101 }
1102
1103 EXPORT int speex_resampler_get_output_latency(SpeexResamplerState *st)
1104 {
1105   return ((st->filt_len / 2) * st->den_rate + (st->num_rate >> 1)) / st->num_rate;
1106 }
1107
1108 EXPORT int speex_resampler_skip_zeros(SpeexResamplerState *st)
1109 {
1110    spx_uint32_t i;
1111    for (i=0;i<st->nb_channels;i++)
1112       st->last_sample[i] = st->filt_len/2;
1113    return RESAMPLER_ERR_SUCCESS;
1114 }
1115
1116 EXPORT int speex_resampler_reset_mem(SpeexResamplerState *st)
1117 {
1118    spx_uint32_t i;
1119    for (i=0;i<st->nb_channels;i++)
1120    {
1121       st->last_sample[i] = 0;
1122       st->magic_samples[i] = 0;
1123       st->samp_frac_num[i] = 0;
1124    }
1125    for (i=0;i<st->nb_channels*(st->filt_len-1);i++)
1126       st->mem[i] = 0;
1127    return RESAMPLER_ERR_SUCCESS;
1128 }
1129
1130 EXPORT const char *speex_resampler_strerror(int err)
1131 {
1132    switch (err)
1133    {
1134       case RESAMPLER_ERR_SUCCESS:
1135          return "Success.";
1136       case RESAMPLER_ERR_ALLOC_FAILED:
1137          return "Memory allocation failed.";
1138       case RESAMPLER_ERR_BAD_STATE:
1139          return "Bad resampler state.";
1140       case RESAMPLER_ERR_INVALID_ARG:
1141          return "Invalid argument.";
1142       case RESAMPLER_ERR_PTR_OVERLAP:
1143          return "Input and output buffers overlap.";
1144       default:
1145          return "Unknown error. Bad error code or strange version mismatch.";
1146    }
1147 }