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