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