src/lib/dl/ext/speex/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       float accum[4] = {0,0,0,0};
 349
 350       for(j=0;j<N;j+=4) {
 351         accum[0] += sinc[j]*iptr[j];
 352         accum[1] += sinc[j+1]*iptr[j+1];
 353         accum[2] += sinc[j+2]*iptr[j+2];
 354         accum[3] += sinc[j+3]*iptr[j+3];
 355       }
 356       sum = accum[0] + accum[1] + accum[2] + accum[3];
 357 #else
 358       sum = inner_product_single(sinc, iptr, N);
 359 #endif
 360
 361       out[out_stride * out_sample++] = PSHR32(sum, 15);
 362       last_sample += int_advance;
 363       samp_frac_num += frac_advance;
 364       if (samp_frac_num >= den_rate)
 365       {
 366          samp_frac_num -= den_rate;
 367          last_sample++;
 368       }
 369    }
 370
 371    st->last_sample[channel_index] = last_sample;
 372    st->samp_frac_num[channel_index] = samp_frac_num;
 373    return out_sample;
 374 }
 375
 376 #ifdef FIXED_POINT
 377 #else
 378 /* This is the same as the previous function, except with a double-precision accumulator */
 379 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)
 380 {
 381    const int N = st->filt_len;
 382    int out_sample = 0;
 383    int last_sample = st->last_sample[channel_index];
 384    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 385    const spx_word16_t *sinc_table = st->sinc_table;
 386    const int out_stride = st->out_stride;
 387    const int int_advance = st->int_advance;
 388    const int frac_advance = st->frac_advance;
 389    const spx_uint32_t den_rate = st->den_rate;
 390    double sum;
 391    int j;
 392
 393    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 394    {
 395       const spx_word16_t *sinc = & sinc_table[samp_frac_num*N];
 396       const spx_word16_t *iptr = & in[last_sample];
 397
 398 #ifndef OVERRIDE_INNER_PRODUCT_DOUBLE
 399       double accum[4] = {0,0,0,0};
 400
 401       for(j=0;j<N;j+=4) {
 402         accum[0] += sinc[j]*iptr[j];
 403         accum[1] += sinc[j+1]*iptr[j+1];
 404         accum[2] += sinc[j+2]*iptr[j+2];
 405         accum[3] += sinc[j+3]*iptr[j+3];
 406       }
 407       sum = accum[0] + accum[1] + accum[2] + accum[3];
 408 #else
 409       sum = inner_product_double(sinc, iptr, N);
 410 #endif
 411
 412       out[out_stride * out_sample++] = PSHR32(sum, 15);
 413       last_sample += int_advance;
 414       samp_frac_num += frac_advance;
 415       if (samp_frac_num >= den_rate)
 416       {
 417          samp_frac_num -= den_rate;
 418          last_sample++;
 419       }
 420    }
 421
 422    st->last_sample[channel_index] = last_sample;
 423    st->samp_frac_num[channel_index] = samp_frac_num;
 424    return out_sample;
 425 }
 426 #endif
 427
 428 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)
 429 {
 430    const int N = st->filt_len;
 431    int out_sample = 0;
 432    int last_sample = st->last_sample[channel_index];
 433    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 434    const int out_stride = st->out_stride;
 435    const int int_advance = st->int_advance;
 436    const int frac_advance = st->frac_advance;
 437    const spx_uint32_t den_rate = st->den_rate;
 438    int j;
 439    spx_word32_t sum;
 440
 441    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 442    {
 443       const spx_word16_t *iptr = & in[last_sample];
 444
 445       const int offset = samp_frac_num*st->oversample/st->den_rate;
 446 #ifdef FIXED_POINT
 447       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
 448 #else
 449       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
 450 #endif
 451       spx_word16_t interp[4];
 452
 453
 454 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_SINGLE
 455       spx_word32_t accum[4] = {0,0,0,0};
 456
 457       for(j=0;j<N;j++) {
 458         const spx_word16_t curr_in=iptr[j];
 459         accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
 460         accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
 461         accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
 462         accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
 463       }
 464
 465       cubic_coef(frac, interp);
 466       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]);
 467 #else
 468       cubic_coef(frac, interp);
 469       sum = interpolate_product_single(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
 470 #endif
 471
 472       out[out_stride * out_sample++] = PSHR32(sum,15);
 473       last_sample += int_advance;
 474       samp_frac_num += frac_advance;
 475       if (samp_frac_num >= den_rate)
 476       {
 477          samp_frac_num -= den_rate;
 478          last_sample++;
 479       }
 480    }
 481
 482    st->last_sample[channel_index] = last_sample;
 483    st->samp_frac_num[channel_index] = samp_frac_num;
 484    return out_sample;
 485 }
 486
 487 #ifdef FIXED_POINT
 488 #else
 489 /* This is the same as the previous function, except with a double-precision accumulator */
 490 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)
 491 {
 492    const int N = st->filt_len;
 493    int out_sample = 0;
 494    int last_sample = st->last_sample[channel_index];
 495    spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 496    const int out_stride = st->out_stride;
 497    const int int_advance = st->int_advance;
 498    const int frac_advance = st->frac_advance;
 499    const spx_uint32_t den_rate = st->den_rate;
 500    int j;
 501    spx_word32_t sum;
 502
 503    while (!(last_sample >= (spx_int32_t)*in_len || out_sample >= (spx_int32_t)*out_len))
 504    {
 505       const spx_word16_t *iptr = & in[last_sample];
 506
 507       const int offset = samp_frac_num*st->oversample/st->den_rate;
 508 #ifdef FIXED_POINT
 509       const spx_word16_t frac = PDIV32(SHL32((samp_frac_num*st->oversample) % st->den_rate,15),st->den_rate);
 510 #else
 511       const spx_word16_t frac = ((float)((samp_frac_num*st->oversample) % st->den_rate))/st->den_rate;
 512 #endif
 513       spx_word16_t interp[4];
 514
 515
 516 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE
 517       double accum[4] = {0,0,0,0};
 518
 519       for(j=0;j<N;j++) {
 520         const double curr_in=iptr[j];
 521         accum[0] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-2]);
 522         accum[1] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset-1]);
 523         accum[2] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset]);
 524         accum[3] += MULT16_16(curr_in,st->sinc_table[4+(j+1)*st->oversample-offset+1]);
 525       }
 526
 527       cubic_coef(frac, interp);
 528       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]);
 529 #else
 530       cubic_coef(frac, interp);
 531       sum = interpolate_product_double(iptr, st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample, interp);
 532 #endif
 533
 534       out[out_stride * out_sample++] = PSHR32(sum,15);
 535       last_sample += int_advance;
 536       samp_frac_num += frac_advance;
 537       if (samp_frac_num >= den_rate)
 538       {
 539          samp_frac_num -= den_rate;
 540          last_sample++;
 541       }
 542    }
 543
 544    st->last_sample[channel_index] = last_sample;
 545    st->samp_frac_num[channel_index] = samp_frac_num;
 546    return out_sample;
 547 }
 548 #endif
 549
 550 static void update_filter(SpeexResamplerState *st)
 551 {
 552    spx_uint32_t old_length;
 553
 554    old_length = st->filt_len;
 555    st->oversample = quality_map[st->quality].oversample;
 556    st->filt_len = quality_map[st->quality].base_length;
 557
 558    if (st->num_rate > st->den_rate)
 559    {
 560       /* down-sampling */
 561       st->cutoff = quality_map[st->quality].downsample_bandwidth * st->den_rate / st->num_rate;
 562       /* FIXME: divide the numerator and denominator by a certain amount if they're too large */
 563       st->filt_len = st->filt_len*st->num_rate / st->den_rate;
 564       /* Round down to make sure we have a multiple of 4 */
 565       st->filt_len &= (~0x3);
 566       if (2*st->den_rate < st->num_rate)
 567          st->oversample >>= 1;
 568       if (4*st->den_rate < st->num_rate)
 569          st->oversample >>= 1;
 570       if (8*st->den_rate < st->num_rate)
 571          st->oversample >>= 1;
 572       if (16*st->den_rate < st->num_rate)
 573          st->oversample >>= 1;
 574       if (st->oversample < 1)
 575          st->oversample = 1;
 576    } else {
 577       /* up-sampling */
 578       st->cutoff = quality_map[st->quality].upsample_bandwidth;
 579    }
 580
 581    /* Choose the resampling type that requires the least amount of memory */
 582    if (st->den_rate <= st->oversample)
 583    {
 584       spx_uint32_t i;
 585       if (!st->sinc_table)
 586          st->sinc_table = (spx_word16_t *)speex_alloc(st->filt_len*st->den_rate*sizeof(spx_word16_t));
 587       else if (st->sinc_table_length < st->filt_len*st->den_rate)
 588       {
 589          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,st->filt_len*st->den_rate*sizeof(spx_word16_t));
 590          st->sinc_table_length = st->filt_len*st->den_rate;
 591       }
 592       for (i=0;i<st->den_rate;i++)
 593       {
 594          spx_int32_t j;
 595          for (j=0;j<st->filt_len;j++)
 596          {
 597             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);
 598          }
 599       }
 600 #ifdef FIXED_POINT
 601       st->resampler_ptr = resampler_basic_direct_single;
 602 #else
 603       if (st->quality>8)
 604          st->resampler_ptr = resampler_basic_direct_double;
 605       else
 606          st->resampler_ptr = resampler_basic_direct_single;
 607 #endif
 608       /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff);*/
 609    } else {
 610       spx_int32_t i;
 611       if (!st->sinc_table)
 612          st->sinc_table = (spx_word16_t *)speex_alloc((st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
 613       else if (st->sinc_table_length < st->filt_len*st->oversample+8)
 614       {
 615          st->sinc_table = (spx_word16_t *)speex_realloc(st->sinc_table,(st->filt_len*st->oversample+8)*sizeof(spx_word16_t));
 616          st->sinc_table_length = st->filt_len*st->oversample+8;
 617       }
 618       for (i=-4;i<(spx_int32_t)(st->oversample*st->filt_len+4);i++)
 619          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);
 620 #ifdef FIXED_POINT
 621       st->resampler_ptr = resampler_basic_interpolate_single;
 622 #else
 623       if (st->quality>8)
 624          st->resampler_ptr = resampler_basic_interpolate_double;
 625       else
 626          st->resampler_ptr = resampler_basic_interpolate_single;
 627 #endif
 628       /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff);*/
 629    }
 630    st->int_advance = st->num_rate/st->den_rate;
 631    st->frac_advance = st->num_rate%st->den_rate;
 632
 633
 634    /* Here's the place where we update the filter memory to take into account
 635       the change in filter length. It's probably the messiest part of the code
 636       due to handling of lots of corner cases. */
 637    if (!st->mem)
 638    {
 639       spx_uint32_t i;
 640       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 641       st->mem = (spx_word16_t*)speex_alloc(st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 642       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
 643          st->mem[i] = 0;
 644       /*speex_warning("init filter");*/
 645    } else if (!st->started)
 646    {
 647       spx_uint32_t i;
 648       st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 649       st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 650       for (i=0;i<st->nb_channels*st->mem_alloc_size;i++)
 651          st->mem[i] = 0;
 652       /*speex_warning("reinit filter");*/
 653    } else if (st->filt_len > old_length)
 654    {
 655       spx_int32_t i;
 656       /* Increase the filter length */
 657       /*speex_warning("increase filter size");*/
 658       int old_alloc_size = st->mem_alloc_size;
 659       if ((st->filt_len-1 + st->buffer_size) > st->mem_alloc_size)
 660       {
 661          st->mem_alloc_size = st->filt_len-1 + st->buffer_size;
 662          st->mem = (spx_word16_t*)speex_realloc(st->mem, st->nb_channels*st->mem_alloc_size * sizeof(spx_word16_t));
 663       }
 664       for (i=st->nb_channels-1;i>=0;i--)
 665       {
 666          spx_int32_t j;
 667          spx_uint32_t olen = old_length;
 668          /*if (st->magic_samples[i])*/
 669          {
 670             /* Try and remove the magic samples as if nothing had happened */
 671
 672             /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
 673             olen = old_length + 2*st->magic_samples[i];
 674             for (j=old_length-2+st->magic_samples[i];j>=0;j--)
 675                st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]] = st->mem[i*old_alloc_size+j];
 676             for (j=0;j<st->magic_samples[i];j++)
 677                st->mem[i*st->mem_alloc_size+j] = 0;
 678             st->magic_samples[i] = 0;
 679          }
 680          if (st->filt_len > olen)
 681          {
 682             /* If the new filter length is still bigger than the "augmented" length */
 683             /* Copy data going backward */
 684             for (j=0;j<olen-1;j++)
 685                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = st->mem[i*st->mem_alloc_size+(olen-2-j)];
 686             /* Then put zeros for lack of anything better */
 687             for (;j<st->filt_len-1;j++)
 688                st->mem[i*st->mem_alloc_size+(st->filt_len-2-j)] = 0;
 689             /* Adjust last_sample */
 690             st->last_sample[i] += (st->filt_len - olen)/2;
 691          } else {
 692             /* Put back some of the magic! */
 693             st->magic_samples[i] = (olen - st->filt_len)/2;
 694             for (j=0;j<st->filt_len-1+st->magic_samples[i];j++)
 695                st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
 696          }
 697       }
 698    } else if (st->filt_len < old_length)
 699    {
 700       spx_uint32_t i;
 701       /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"
 702          samples so they can be used directly as input the next time(s) */
 703       for (i=0;i<st->nb_channels;i++)
 704       {
 705          spx_uint32_t j;
 706          spx_uint32_t old_magic = st->magic_samples[i];
 707          st->magic_samples[i] = (old_length - st->filt_len)/2;
 708          /* We must copy some of the memory that's no longer used */
 709          /* Copy data going backward */
 710          for (j=0;j<st->filt_len-1+st->magic_samples[i]+old_magic;j++)
 711             st->mem[i*st->mem_alloc_size+j] = st->mem[i*st->mem_alloc_size+j+st->magic_samples[i]];
 712          st->magic_samples[i] += old_magic;
 713       }
 714    }
 715
 716 }
 717
 718 EXPORT SpeexResamplerState *speex_resampler_init(spx_uint32_t nb_channels, spx_uint32_t in_rate, spx_uint32_t out_rate, int quality, int *err)
 719 {
 720    return speex_resampler_init_frac(nb_channels, in_rate, out_rate, in_rate, out_rate, quality, err);
 721 }
 722
 723 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)
 724 {
 725    spx_uint32_t i;
 726    SpeexResamplerState *st;
 727    if (quality > 10 || quality < 0)
 728    {
 729       if (err)
 730          *err = RESAMPLER_ERR_INVALID_ARG;
 731       return NULL;
 732    }
 733    st = (SpeexResamplerState *)speex_alloc(sizeof(SpeexResamplerState));
 734    st->initialised = 0;
 735    st->started = 0;
 736    st->in_rate = 0;
 737    st->out_rate = 0;
 738    st->num_rate = 0;
 739    st->den_rate = 0;
 740    st->quality = -1;
 741    st->sinc_table_length = 0;
 742    st->mem_alloc_size = 0;
 743    st->filt_len = 0;
 744    st->mem = 0;
 745    st->resampler_ptr = 0;
 746
 747    st->cutoff = 1.f;
 748    st->nb_channels = nb_channels;
 749    st->in_stride = 1;
 750    st->out_stride = 1;
 751
 752 #ifdef FIXED_POINT
 753    st->buffer_size = 160;
 754 #else
 755    st->buffer_size = 160;
 756 #endif
 757
 758    /* Per channel data */
 759    st->last_sample = (spx_int32_t*)speex_alloc(nb_channels*sizeof(int));
 760    st->magic_samples = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
 761    st->samp_frac_num = (spx_uint32_t*)speex_alloc(nb_channels*sizeof(int));
 762    for (i=0;i<nb_channels;i++)
 763    {
 764       st->last_sample[i] = 0;
 765       st->magic_samples[i] = 0;
 766       st->samp_frac_num[i] = 0;
 767    }
 768
 769    speex_resampler_set_quality(st, quality);
 770    speex_resampler_set_rate_frac(st, ratio_num, ratio_den, in_rate, out_rate);
 771
 772
 773    update_filter(st);
 774
 775    st->initialised = 1;
 776    if (err)
 777       *err = RESAMPLER_ERR_SUCCESS;
 778
 779    return st;
 780 }
 781
 782 EXPORT void speex_resampler_destroy(SpeexResamplerState *st)
 783 {
 784    speex_free(st->mem);
 785    speex_free(st->sinc_table);
 786    speex_free(st->last_sample);
 787    speex_free(st->magic_samples);
 788    speex_free(st->samp_frac_num);
 789    speex_free(st);
 790 }
 791
 792 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)
 793 {
 794    int j=0;
 795    const int N = st->filt_len;
 796    int out_sample = 0;
 797    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
 798    spx_uint32_t ilen;
 799
 800    st->started = 1;
 801
 802    /* Call the right resampler through the function ptr */
 803    out_sample = st->resampler_ptr(st, channel_index, mem, in_len, out, out_len);
 804
 805    if (st->last_sample[channel_index] < (spx_int32_t)*in_len)
 806       *in_len = st->last_sample[channel_index];
 807    *out_len = out_sample;
 808    st->last_sample[channel_index] -= *in_len;
 809
 810    ilen = *in_len;
 811
 812    for(j=0;j<N-1;++j)
 813      mem[j] = mem[j+ilen];
 814
 815    return RESAMPLER_ERR_SUCCESS;
 816 }
 817
 818 static int speex_resampler_magic(SpeexResamplerState *st, spx_uint32_t channel_index, spx_word16_t **out, spx_uint32_t out_len) {
 819    spx_uint32_t tmp_in_len = st->magic_samples[channel_index];
 820    spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
 821    const int N = st->filt_len;
 822
 823    speex_resampler_process_native(st, channel_index, &tmp_in_len, *out, &out_len);
 824
 825    st->magic_samples[channel_index] -= tmp_in_len;
 826
 827    /* If we couldn't process all "magic" input samples, save the rest for next time */
 828    if (st->magic_samples[channel_index])
 829    {
 830       spx_uint32_t i;
 831       for (i=0;i<st->magic_samples[channel_index];i++)
 832          mem[N-1+i]=mem[N-1+i+tmp_in_len];
 833    }
 834    *out += out_len*st->out_stride;
 835    return out_len;
 836 }
 837
 838 #ifdef FIXED_POINT
 839 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)
 840 #else
 841 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)
 842 #endif
 843 {
 844    int j;
 845    spx_uint32_t ilen = *in_len;
 846    spx_uint32_t olen = *out_len;
 847    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
 848    const int filt_offs = st->filt_len - 1;
 849    const spx_uint32_t xlen = st->mem_alloc_size - filt_offs;
 850    const int istride = st->in_stride;
 851
 852    if (st->magic_samples[channel_index])
 853       olen -= speex_resampler_magic(st, channel_index, &out, olen);
 854    if (! st->magic_samples[channel_index]) {
 855       while (ilen && olen) {
 856         spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
 857         spx_uint32_t ochunk = olen;
 858
 859         if (in) {
 860            for(j=0;j<ichunk;++j)
 861               x[j+filt_offs]=in[j*istride];
 862         } else {
 863           for(j=0;j<ichunk;++j)
 864             x[j+filt_offs]=0;
 865         }
 866         speex_resampler_process_native(st, channel_index, &ichunk, out, &ochunk);
 867         ilen -= ichunk;
 868         olen -= ochunk;
 869         out += ochunk * st->out_stride;
 870         if (in)
 871            in += ichunk * istride;
 872       }
 873    }
 874    *in_len -= ilen;
 875    *out_len -= olen;
 876    return RESAMPLER_ERR_SUCCESS;
 877 }
 878
 879 #ifdef FIXED_POINT
 880 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)
 881 #else
 882 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)
 883 #endif
 884 {
 885    int j;
 886    const int istride_save = st->in_stride;
 887    const int ostride_save = st->out_stride;
 888    spx_uint32_t ilen = *in_len;
 889    spx_uint32_t olen = *out_len;
 890    spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
 891    const spx_uint32_t xlen = st->mem_alloc_size - (st->filt_len - 1);
 892 #ifdef VAR_ARRAYS
 893    const unsigned int ylen = (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;
 894    VARDECL(spx_word16_t *ystack);
 895    ALLOC(ystack, ylen, spx_word16_t);
 896 #else
 897    const unsigned int ylen = FIXED_STACK_ALLOC;
 898    spx_word16_t ystack[FIXED_STACK_ALLOC];
 899 #endif
 900
 901    st->out_stride = 1;
 902
 903    while (ilen && olen) {
 904      spx_word16_t *y = ystack;
 905      spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
 906      spx_uint32_t ochunk = (olen > ylen) ? ylen : olen;
 907      spx_uint32_t omagic = 0;
 908
 909      if (st->magic_samples[channel_index]) {
 910        omagic = speex_resampler_magic(st, channel_index, &y, ochunk);
 911        ochunk -= omagic;
 912        olen -= omagic;
 913      }
 914      if (! st->magic_samples[channel_index]) {
 915        if (in) {
 916          for(j=0;j<ichunk;++j)
 917 #ifdef FIXED_POINT
 918            x[j+st->filt_len-1]=WORD2INT(in[j*istride_save]);
 919 #else
 920            x[j+st->filt_len-1]=in[j*istride_save];
 921 #endif
 922        } else {
 923          for(j=0;j<ichunk;++j)
 924            x[j+st->filt_len-1]=0;
 925        }
 926
 927        speex_resampler_process_native(st, channel_index, &ichunk, y, &ochunk);
 928      } else {
 929        ichunk = 0;
 930        ochunk = 0;
 931      }
 932
 933      for (j=0;j<ochunk+omagic;++j)
 934 #ifdef FIXED_POINT
 935         out[j*ostride_save] = ystack[j];
 936 #else
 937         out[j*ostride_save] = WORD2INT(ystack[j]);
 938 #endif
 939
 940      ilen -= ichunk;
 941      olen -= ochunk;
 942      out += (ochunk+omagic) * ostride_save;
 943      if (in)
 944        in += ichunk * istride_save;
 945    }
 946    st->out_stride = ostride_save;
 947    *in_len -= ilen;
 948    *out_len -= olen;
 949
 950    return RESAMPLER_ERR_SUCCESS;
 951 }
 952
 953 EXPORT int speex_resampler_process_interleaved_float(SpeexResamplerState *st, const float *in, spx_uint32_t *in_len, float *out, spx_uint32_t *out_len)
 954 {
 955    spx_uint32_t i;
 956    int istride_save, ostride_save;
 957    spx_uint32_t bak_len = *out_len;
 958    istride_save = st->in_stride;
 959    ostride_save = st->out_stride;
 960    st->in_stride = st->out_stride = st->nb_channels;
 961    for (i=0;i<st->nb_channels;i++)
 962    {
 963       *out_len = bak_len;
 964       if (in != NULL)
 965          speex_resampler_process_float(st, i, in+i, in_len, out+i, out_len);
 966       else
 967          speex_resampler_process_float(st, i, NULL, in_len, out+i, out_len);
 968    }
 969    st->in_stride = istride_save;
 970    st->out_stride = ostride_save;
 971    return RESAMPLER_ERR_SUCCESS;
 972 }
 973
 974 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)
 975 {
 976    spx_uint32_t i;
 977    int istride_save, ostride_save;
 978    spx_uint32_t bak_len = *out_len;
 979    istride_save = st->in_stride;
 980    ostride_save = st->out_stride;
 981    st->in_stride = st->out_stride = st->nb_channels;
 982    for (i=0;i<st->nb_channels;i++)
 983    {
 984       *out_len = bak_len;
 985       if (in != NULL)
 986          speex_resampler_process_int(st, i, in+i, in_len, out+i, out_len);
 987       else
 988          speex_resampler_process_int(st, i, NULL, in_len, out+i, out_len);
 989    }
 990    st->in_stride = istride_save;
 991    st->out_stride = ostride_save;
 992    return RESAMPLER_ERR_SUCCESS;
 993 }
 994
 995 EXPORT int speex_resampler_set_rate(SpeexResamplerState *st, spx_uint32_t in_rate, spx_uint32_t out_rate)
 996 {
 997    return speex_resampler_set_rate_frac(st, in_rate, out_rate, in_rate, out_rate);
 998 }
 999
1000 EXPORT void speex_resampler_get_rate(SpeexResamplerState *st, spx_uint32_t *in_rate, spx_uint32_t *out_rate)
1001 {
1002    *in_rate = st->in_rate;
1003    *out_rate = st->out_rate;
1004 }
1005
1006 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)
1007 {
1008    spx_uint32_t fact;
1009    spx_uint32_t old_den;
1010    spx_uint32_t i;
1011    if (st->in_rate == in_rate && st->out_rate == out_rate && st->num_rate == ratio_num && st->den_rate == ratio_den)
1012       return RESAMPLER_ERR_SUCCESS;
1013
1014    old_den = st->den_rate;
1015    st->in_rate = in_rate;
1016    st->out_rate = out_rate;
1017    st->num_rate = ratio_num;
1018    st->den_rate = ratio_den;
1019    /* FIXME: This is terribly inefficient, but who cares (at least for now)? */
1020    for (fact=2;fact<=IMIN(st->num_rate, st->den_rate);fact++)
1021    {
1022       while ((st->num_rate % fact == 0) && (st->den_rate % fact == 0))
1023       {
1024          st->num_rate /= fact;
1025          st->den_rate /= fact;
1026       }
1027    }
1028
1029    if (old_den > 0)
1030    {
1031       for (i=0;i<st->nb_channels;i++)
1032       {
1033          st->samp_frac_num[i]=st->samp_frac_num[i]*st->den_rate/old_den;
1034          /* Safety net */
1035          if (st->samp_frac_num[i] >= st->den_rate)
1036             st->samp_frac_num[i] = st->den_rate-1;
1037       }
1038    }
1039
1040    if (st->initialised)
1041       update_filter(st);
1042    return RESAMPLER_ERR_SUCCESS;
1043 }
1044
1045 EXPORT void speex_resampler_get_ratio(SpeexResamplerState *st, spx_uint32_t *ratio_num, spx_uint32_t *ratio_den)
1046 {
1047    *ratio_num = st->num_rate;
1048    *ratio_den = st->den_rate;
1049 }
1050
1051 EXPORT int speex_resampler_set_quality(SpeexResamplerState *st, int quality)
1052 {
1053    if (quality > 10 || quality < 0)
1054       return RESAMPLER_ERR_INVALID_ARG;
1055    if (st->quality == quality)
1056       return RESAMPLER_ERR_SUCCESS;
1057    st->quality = quality;
1058    if (st->initialised)
1059       update_filter(st);
1060    return RESAMPLER_ERR_SUCCESS;
1061 }
1062
1063 EXPORT void speex_resampler_get_quality(SpeexResamplerState *st, int *quality)
1064 {
1065    *quality = st->quality;
1066 }
1067
1068 EXPORT void speex_resampler_set_input_stride(SpeexResamplerState *st, spx_uint32_t stride)
1069 {
1070    st->in_stride = stride;
1071 }
1072
1073 EXPORT void speex_resampler_get_input_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1074 {
1075    *stride = st->in_stride;
1076 }
1077
1078 EXPORT void speex_resampler_set_output_stride(SpeexResamplerState *st, spx_uint32_t stride)
1079 {
1080    st->out_stride = stride;
1081 }
1082
1083 EXPORT void speex_resampler_get_output_stride(SpeexResamplerState *st, spx_uint32_t *stride)
1084 {
1085    *stride = st->out_stride;
1086 }
1087
1088 EXPORT int speex_resampler_get_input_latency(SpeexResamplerState *st)
1089 {
1090   return st->filt_len / 2;
1091 }
1092
1093 EXPORT int speex_resampler_get_output_latency(SpeexResamplerState *st)
1094 {
1095   return ((st->filt_len / 2) * st->den_rate + (st->num_rate >> 1)) / st->num_rate;
1096 }
1097
1098 EXPORT int speex_resampler_skip_zeros(SpeexResamplerState *st)
1099 {
1100    spx_uint32_t i;
1101    for (i=0;i<st->nb_channels;i++)
1102       st->last_sample[i] = st->filt_len/2;
1103    return RESAMPLER_ERR_SUCCESS;
1104 }
1105
1106 EXPORT int speex_resampler_reset_mem(SpeexResamplerState *st)
1107 {
1108    spx_uint32_t i;
1109    for (i=0;i<st->nb_channels*(st->filt_len-1);i++)
1110       st->mem[i] = 0;
1111    return RESAMPLER_ERR_SUCCESS;
1112 }
1113
1114 EXPORT const char *speex_resampler_strerror(int err)
1115 {
1116    switch (err)
1117    {
1118       case RESAMPLER_ERR_SUCCESS:
1119          return "Success.";
1120       case RESAMPLER_ERR_ALLOC_FAILED:
1121          return "Memory allocation failed.";
1122       case RESAMPLER_ERR_BAD_STATE:
1123          return "Bad resampler state.";
1124       case RESAMPLER_ERR_INVALID_ARG:
1125          return "Invalid argument.";
1126       case RESAMPLER_ERR_PTR_OVERLAP:
1127          return "Input and output buffers overlap.";
1128       default:
1129          return "Unknown error. Bad error code or strange version mismatch.";
1130    }
1131 }