libspeex/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.14159263
  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    int j;
 341
 342    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 343    {
 344       const spx_word16_t *sinc = & sinc_table[samp_frac_num*N];
 345       const spx_word16_t *iptr = & in[last_sample];
 346
 347 #ifndef OVERRIDE_INNER_PRODUCT_SINGLE
 348       sum = 0;
 349       for(j=0;j<N;j++) sum += MULT16_16(sinc[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(sinc[j], iptr[j]);
 357         accum[1] += MULT16_16(sinc[j+1], iptr[j+1]);
 358         accum[2] += MULT16_16(sinc[j+2], iptr[j+2]);
 359         accum[3] += MULT16_16(sinc[j+3], iptr[j+3]);
 360       }
 361       sum = accum[0] + accum[1] + accum[2] + accum[3];
 362 */
 363 #else
 364       sum = inner_product_single(sinc, 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    int j;
 398
 399    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 400    {
 401       const spx_word16_t *sinc = & sinc_table[samp_frac_num*N];
 402       const spx_word16_t *iptr = & in[last_sample];
 403
 404 #ifndef OVERRIDE_INNER_PRODUCT_DOUBLE
 405       double accum[4] = {0,0,0,0};
 406
 407       for(j=0;j<N;j+=4) {
 408         accum[0] += sinc[j]*iptr[j];
 409         accum[1] += sinc[j+1]*iptr[j+1];
 410         accum[2] += sinc[j+2]*iptr[j+2];
 411         accum[3] += sinc[j+3]*iptr[j+3];
 412       }
 413       sum = accum[0] + accum[1] + accum[2] + accum[3];
 414 #else
 415       sum = inner_product_double(sinc, 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    int j;
 445    spx_word32_t sum;
 446
 447    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 448    {
 449       const spx_word16_t *iptr = & in[last_sample];
 450
 451       const int offset = samp_frac_num*st->oversample/st->den_rate;
 452 #ifdef FIXED_POINT
 453       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
 454 #else
 455       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
 456 #endif
 457       spx_word16_t interp[4];
 458
 459
 460 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_SINGLE
 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    int j;
 507    spx_word32_t sum;
 508
 509    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 510    {
 511       const spx_word16_t *iptr = & in[last_sample];
 512
 513       const int offset = samp_frac_num*st->oversample/st->den_rate;
 514 #ifdef FIXED_POINT
 515       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
 516 #else
 517       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
 518 #endif
 519       spx_word16_t interp[4];
 520
 521
 522 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE
 523       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 down to make sure we have a multiple of 4 */
 571       st->filt_len &= (~0x3);
 572       if (2*st->den_rate < st->num_rate)
 573          st->oversample >>= 1;
 574       if (4*st->den_rate < st->num_rate)
 575          st->oversample >>= 1;
 576       if (8*st->den_rate < st->num_rate)
 577          st->oversample >>= 1;
 578       if (16*st->den_rate < st->num_rate)
 579          st->oversample >>= 1;
 580       if (st->oversample < 1)
 581          st->oversample = 1;
 582    } else {
 583       /* up-sampling */
 584       st->cutoff = quality_map[st->quality].upsample_bandwidth;
 585    }
 586
 587    /* Choose the resampling type that requires the least amount of memory */
 588    if (st->den_rate <= st->oversample)
 589    {
 590       spx_uint32_t i;
 591       if (!st->sinc_table)
 592          st->sinc_table = (spx_word16_t *)speex_alloc(st->filt_len*st->den_rate*sizeof(spx_word16_t));
 593       else if (st->sinc_table_length < st->filt_len*st->den_rate)
 594       {
 595          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,st->filt_len*st->den_rate*sizeof(spx_word16_t));
 596          st->sinc_table_length = st->filt_len*st->den_rate;
 597       }
 598       for (i=0;i<st->den_rate;i++)
 599       {
 600          spx_int32_t j;
 601          for (j=0;j<st->filt_len;j++)
 602          {
 603             st->sinc_table[i*st->filt_len+j] = sinc(st->cutoff,((j-(spx_int32_t)st->filt_len/2+1)-((float)i)/st->den_rate), st->filt_len, quality_map[st->quality].window_func);
 604          }
 605       }
 606 #ifdef FIXED_POINT
 607       st->resampler_ptr = resampler_basic_direct_single;
 608 #else
 609       if (st->quality>8)
 610          st->resampler_ptr = resampler_basic_direct_double;
 611       else
 612          st->resampler_ptr = resampler_basic_direct_single;
 613 #endif
 614       /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff);*/
 615    } else {
 616       spx_int32_t i;
 617       if (!st->sinc_table)
 618          st->sinc_table = (spx_word16_t *)speex_alloc((st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
 619       else if (st->sinc_table_length < st->filt_len*st->oversample+8)
 620       {
 621          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,(st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
 622          st->sinc_table_length = st->filt_len*st->oversample+8;
 623       }
 624       for (i=-4;i<(spx_int32_t)(st->oversample*st->filt_len+4);i++)
 625          st->sinc_table[i+4] = sinc(st->cutoff,(i/(float)st->oversample - st->filt_len/2), st->filt_len, quality_map[st->quality].window_func);
 626 #ifdef FIXED_POINT
 627       st->resampler_ptr = resampler_basic_interpolate_single;
 628 #else
 629       if (st->quality>8)
 630          st->resampler_ptr = resampler_basic_interpolate_double;
 631       else
 632          st->resampler_ptr = resampler_basic_interpolate_single;
 633 #endif
 634       /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff);*/
 635    }
 636    st->int_advance = st->num_rate/st->den_rate;
 637    st->frac_advance = st->num_rate%st->den_rate;
 638
 639
 640    /* Here's the place where we update the filter memory to take into account
 641       the change in filter length. It's probably the messiest part of the code
 642       due to handling of lots of corner cases. */
 643    if (!st->mem)
 644    {
 645       spx_uint32_t i;
 646       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 647       st->mem = (spx_word16_t*)speex_alloc(st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 648       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
 649          st->mem[i] = 0;
 650       /*speex_warning("init filter");*/
 651    } else if (!st->started)
 652    {
 653       spx_uint32_t i;
 654       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 655       st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 656       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
 657          st->mem[i] = 0;
 658       /*speex_warning("reinit filter");*/
 659    } else if (st->filt_len > old_length)
 660    {
 661       spx_int32_t i;
 662       /* Increase the filter length */
 663       /*speex_warning("increase filter size");*/
 664       int old_alloc_size = st->mem_alloc_size;
 665       if ((st->filt_len-1 + st->buffer_size) > st->mem_alloc_size)
 666       {
 667          st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 668          st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 669       }
 670       for (i=st->nb_channels-1;i>=0;i--)
 671       {
 672          spx_int32_t j;
 673          spx_uint32_t olen = old_length;
 674          /*if (st->magic_samples[i])*/
 675          {
 676             /* Try and remove the magic samples as if nothing had happened */
 677
 678             /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
 679             olen = old_length + 2*st->magic_samples[i];
 680             for (j=old_length-2+st->magic_samples[i];j>=0;j--)
 681                st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]] = st->mem[i*old_alloc_size+j];
 682             for (j=0;j<st->magic_samples[i];j++)
 683                st->mem[i*st->mem_alloc_size+j] = 0;
 684             st->magic_samples[i] = 0;
 685          }
 686          if (st->filt_len > olen)
 687          {
 688             /* If the new filter length is still bigger than the "augmented" length */
 689             /* Copy data going backward */
 690             for (j=0;j<olen-1;j++)
 691                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = st->mem[i*st->mem_alloc_size+(olen-2-j)];
 692             /* Then put zeros for lack of anything better */
 693             for (;j<st->filt_len-1;j++)
 694                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = 0;
 695             /* Adjust last_sample */
 696             st->last_sample[i] += (st->filt_len - olen)/2;
 697          } else {
 698             /* Put back some of the magic! */
 699             st->magic_samples[i] = (olen - st->filt_len)/2;
 700             for (j=0;j<st->filt_len-1+st->magic_samples[i];j++)
 701                st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
 702          }
 703       }
 704    } else if (st->filt_len < old_length)
 705    {
 706       spx_uint32_t i;
 707       /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"
 708          samples so they can be used directly as input the next time(s) */
 709       for (i=0;i<st->nb_channels;i++)
 710       {
 711          spx_uint32_t j;
 712          spx_uint32_t old_magic = st->magic_samples[i];
 713          st->magic_samples[i] = (old_length - st->filt_len)/2;
 714          /* We must copy some of the memory that's no longer used */
 715          /* Copy data going backward */
 716          for (j=0;j<st->filt_len-1+st->magic_samples[i]+old_magic;j++)
 717             st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
 718          st->magic_samples[i] += old_magic;
 719       }
 720    }
 721
 722 }
 723
 724 EXPORT SpeexResamplerState *speex_resampler_init(spx_uint32_t nb_channels, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
 725 {
 726    return speex_resampler_init_frac(nb_channels, in_rate, out_rate, in_rate, out_rate, quality, err);
 727 }
 728
 729 EXPORT SpeexResamplerState *speex_resampler_init_frac(spx_uint32_t nb_channels, spx_uint32_t ratio_num, spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
 730 {
 731    spx_uint32_t i;
 732    SpeexResamplerState *st;
 733    if (quality > 10 || quality < 0)
 734    {
 735       if (err)
 736          *err = RESAMPLER_ERR_INVALID_ARG;
 737       return NULL;
 738    }
 739    st = (SpeexResamplerState *)speex_alloc(sizeof(SpeexResamplerState));
 740    st->initialised = 0;
 741    st->started = 0;
 742    st->in_rate = 0;
 743    st->out_rate = 0;
 744    st->num_rate = 0;
 745    st->den_rate = 0;
 746    st->quality = -1;
 747    st->sinc_table_length = 0;
 748    st->mem_alloc_size = 0;
 749    st->filt_len = 0;
 750    st->mem = 0;
 751    st->resampler_ptr = 0;
 752
 753    st->cutoff = 1.f;
 754    st->nb_channels = nb_channels;
 755    st->in_stride = 1;
 756    st->out_stride = 1;
 757
 758 #ifdef FIXED_POINT
 759    st->buffer_size = 160;
 760 #else
 761    st->buffer_size = 160;
 762 #endif
 763
 764    /* Per channel data */
 765    st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(int));
 766    st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
 767    st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
 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_len = *out_len;
 964    istride_save = st->in_stride;
 965    ostride_save = st->out_stride;
 966    st->in_stride = st->out_stride = st->nb_channels;
 967    for (i=0;i<st->nb_channels;i++)
 968    {
 969       *out_len = bak_len;
 970       if (in != NULL)
 971          speex_resampler_process_float(st, i, in+i, in_len, out+i, out_len);
 972       else
 973          speex_resampler_process_float(st, i, NULL, in_len, out+i, out_len);
 974    }
 975    st->in_stride = istride_save;
 976    st->out_stride = ostride_save;
 977    return RESAMPLER_ERR_SUCCESS;
 978 }
 979
 980 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)
 981 {
 982    spx_uint32_t i;
 983    int istride_save, ostride_save;
 984    spx_uint32_t bak_len = *out_len;
 985    istride_save = st->in_stride;
 986    ostride_save = st->out_stride;
 987    st->in_stride = st->out_stride = st->nb_channels;
 988    for (i=0;i<st->nb_channels;i++)
 989    {
 990       *out_len = bak_len;
 991       if (in != NULL)
 992          speex_resampler_process_int(st, i, in+i, in_len, out+i, out_len);
 993       else
 994          speex_resampler_process_int(st, i, NULL, in_len, out+i, out_len);
 995    }
 996    st->in_stride = istride_save;
 997    st->out_stride = ostride_save;
 998    return RESAMPLER_ERR_SUCCESS;
 999 }
1000
1001 EXPORT int speex_resampler_set_rate(SpeexResamplerState *st, spx_uint32_t in_rate, spx_uint32_t out_rate)
1002 {
1003    return speex_resampler_set_rate_frac(st, in_rate, out_rate, in_rate, out_rate);
1004 }
1005
1006 EXPORT void speex_resampler_get_rate(SpeexResamplerState *st, spx_uint32_t *in_rate, spx_uint32_t *out_rate)
1007 {
1008    *in_rate = st->in_rate;
1009    *out_rate = st->out_rate;
1010 }
1011
1012 EXPORT int speex_resampler_set_rate_frac(SpeexResamplerState *st, spx_uint32_t ratio_num, spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate)
1013 {
1014    spx_uint32_t fact;
1015    spx_uint32_t old_den;
1016    spx_uint32_t i;
1017    if (st->in_rate == in_rate && st->out_rate == out_rate && st->num_rate == ratio_num && st->den_rate == ratio_den)
1018       return RESAMPLER_ERR_SUCCESS;
1019
1020    old_den = st->den_rate;
1021    st->in_rate = in_rate;
1022    st->out_rate = out_rate;
1023    st->num_rate = ratio_num;
1024    st->den_rate = ratio_den;
1025    /* FIXME: This is terribly inefficient, but who cares (at least for now)? */
1026    for (fact=2;fact<=IMIN(st->num_rate, st->den_rate);fact++)
1027    {
1028       while ((st->num_rate % fact == 0) && (st->den_rate % fact == 0))
1029       {
1030          st->num_rate /= fact;
1031          st->den_rate /= fact;
1032       }
1033    }
1034
1035    if (old_den > 0)
1036    {
1037       for (i=0;i<st->nb_channels;i++)
1038       {
1039          st->samp_frac_num[i]=st->samp_frac_num[i]*st->den_rate/old_den;
1040          /* Safety net */
1041          if (st->samp_frac_num[i] >= st->den_rate)
1042             st->samp_frac_num[i] = st->den_rate-1;
1043       }
1044    }
1045
1046    if (st->initialised)
1047       update_filter(st);
1048    return RESAMPLER_ERR_SUCCESS;
1049 }
1050
1051 EXPORT void speex_resampler_get_ratio(SpeexResamplerState *st, spx_uint32_t *ratio_num, spx_uint32_t *ratio_den)
1052 {
1053    *ratio_num = st->num_rate;
1054    *ratio_den = st->den_rate;
1055 }
1056
1057 EXPORT int speex_resampler_set_quality(SpeexResamplerState *st, int quality)
1058 {
1059    if (quality > 10 || quality < 0)
1060       return RESAMPLER_ERR_INVALID_ARG;
1061    if (st->quality == quality)
1062       return RESAMPLER_ERR_SUCCESS;
1063    st->quality = quality;
1064    if (st->initialised)
1065       update_filter(st);
1066    return RESAMPLER_ERR_SUCCESS;
1067 }
1068
1069 EXPORT void speex_resampler_get_quality(SpeexResamplerState *st, int *quality)
1070 {
1071    *quality = st->quality;
1072 }
1073
1074 EXPORT void speex_resampler_set_input_stride(SpeexResamplerState *st, spx_uint32_t stride)
1075 {
1076    st->in_stride = stride;
1077 }
1078
1079 EXPORT void speex_resampler_get_input_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1080 {
1081    *stride = st->in_stride;
1082 }
1083
1084 EXPORT void speex_resampler_set_output_stride(SpeexResamplerState *st, spx_uint32_t stride)
1085 {
1086    st->out_stride = stride;
1087 }
1088
1089 EXPORT void speex_resampler_get_output_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1090 {
1091    *stride = st->out_stride;
1092 }
1093
1094 EXPORT int speex_resampler_get_input_latency(SpeexResamplerState *st)
1095 {
1096   return st->filt_len / 2;
1097 }
1098
1099 EXPORT int speex_resampler_get_output_latency(SpeexResamplerState *st)
1100 {
1101   return ((st->filt_len / 2) * st->den_rate + (st->num_rate >> 1)) / st->num_rate;
1102 }
1103
1104 EXPORT int speex_resampler_skip_zeros(SpeexResamplerState *st)
1105 {
1106    spx_uint32_t i;
1107    for (i=0;i<st->nb_channels;i++)
1108       st->last_sample[i] = st->filt_len/2;
1109    return RESAMPLER_ERR_SUCCESS;
1110 }
1111
1112 EXPORT int speex_resampler_reset_mem(SpeexResamplerState *st)
1113 {
1114    spx_uint32_t i;
1115    for (i=0;i<st->nb_channels*(st->filt_len-1);i++)
1116       st->mem[i] = 0;
1117    return RESAMPLER_ERR_SUCCESS;
1118 }
1119
1120 EXPORT const char *speex_resampler_strerror(int err)
1121 {
1122    switch (err)
1123    {
1124       case RESAMPLER_ERR_SUCCESS:
1125          return "Success.";
1126       case RESAMPLER_ERR_ALLOC_FAILED:
1127          return "Memory allocation failed.";
1128       case RESAMPLER_ERR_BAD_STATE:
1129          return "Bad resampler state.";
1130       case RESAMPLER_ERR_INVALID_ARG:
1131          return "Invalid argument.";
1132       case RESAMPLER_ERR_PTR_OVERLAP:
1133          return "Input and output buffers overlap.";
1134       default:
1135          return "Unknown error. Bad error code or strange version mismatch.";
1136    }
1137 }