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