gst/audioresample/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
  67 #ifdef HAVE_STRING_H
  68 #include <string.h>
  69 #endif
  70
  71 #include <glib.h>
  72
  73 #ifdef HAVE_ORC
  74 #include <orc/orc.h>
  75 #endif
  76
  77 #define EXPORT G_GNUC_INTERNAL
  78
  79 #ifdef _USE_SSE
  80 #ifndef HAVE_XMMINTRIN_H
  81 #undef _USE_SSE
  82 #endif
  83 #endif
  84
  85 #ifdef _USE_SSE2
  86 #ifndef HAVE_EMMINTRIN_H
  87 #undef _USE_SSE2
  88 #endif
  89 #endif
  90
  91 #ifdef _USE_NEON
  92 #ifndef HAVE_ARM_NEON
  93 #undef _USE_NEON
  94 #endif
  95 #endif
  96
  97 static inline void *
  98 speex_alloc (int size)
  99 {
 100   return g_malloc0 (size);
 101 }
 102
 103 static inline void *
 104 speex_realloc (void *ptr, int size)
 105 {
 106   return g_realloc (ptr, size);
 107 }
 108
 109 static inline void
 110 speex_free (void *ptr)
 111 {
 112   g_free (ptr);
 113 }
 114
 115 #include "speex_resampler.h"
 116 #include "arch.h"
 117 #else /* OUTSIDE_SPEEX */
 118
 119 #include "../include/speex/speex_resampler.h"
 120 #include "arch.h"
 121 #include "os_support.h"
 122 #endif /* OUTSIDE_SPEEX */
 123
 124 #include <math.h>
 125
 126 #ifdef FIXED_POINT
 127 #define WORD2INT(x) ((x) < -32767 ? -32768 : ((x) > 32766 ? 32767 : (x)))
 128 #else
 129 #define WORD2INT(x) ((x) < -32767.5f ? -32768 : ((x) > 32766.5f ? 32767 : floor(.5+(x))))
 130 #endif
 131
 132 #define IMAX(a,b) ((a) > (b) ? (a) : (b))
 133 #define IMIN(a,b) ((a) < (b) ? (a) : (b))
 134
 135 #ifndef NULL
 136 #define NULL 0
 137 #endif
 138
 139 #if defined _USE_SSE || defined _USE_SSE2
 140 #include "resample_sse.h"
 141 #endif
 142
 143 #ifdef _USE_NEON
 144 #include "resample_neon.h"
 145 #endif
 146
 147 /* Numer of elements to allocate on the stack */
 148 #ifdef VAR_ARRAYS
 149 #define FIXED_STACK_ALLOC 8192
 150 #else
 151 #define FIXED_STACK_ALLOC 1024
 152 #endif
 153
 154 /* Allow selecting SSE or not when compiled with SSE support */
 155 #ifdef _USE_SSE
 156 #define SSE_FALLBACK(macro) \
 157   if (st->use_sse) goto sse_##macro##_sse; {
 158 #define SSE_IMPLEMENTATION(macro) \
 159   goto sse_##macro##_end; } sse_##macro##_sse: {
 160 #define SSE_END(macro) sse_##macro##_end:; }
 161 #else
 162 #define SSE_FALLBACK(macro)
 163 #endif
 164
 165 #ifdef _USE_SSE2
 166 #define SSE2_FALLBACK(macro) \
 167   if (st->use_sse2) goto sse2_##macro##_sse2; {
 168 #define SSE2_IMPLEMENTATION(macro) \
 169   goto sse2_##macro##_end; } sse2_##macro##_sse2: {
 170 #define SSE2_END(macro) sse2_##macro##_end:; }
 171 #else
 172 #define SSE2_FALLBACK(macro)
 173 #endif
 174
 175 #ifdef _USE_NEON
 176 #define NEON_FALLBACK(macro) \
 177   if (st->use_neon) goto neon_##macro##_neon; {
 178 #define NEON_IMPLEMENTATION(macro) \
 179   goto neon_##macro##_end; } neon_##macro##_neon: {
 180 #define NEON_END(macro) neon_##macro##_end:; }
 181 #else
 182 #define NEON_FALLBACK(macro)
 183 #endif
 184
 185
 186 typedef int (*resampler_basic_func) (SpeexResamplerState *, spx_uint32_t,
 187     const spx_word16_t *, spx_uint32_t *, spx_word16_t *, spx_uint32_t *);
 188
 189 struct SpeexResamplerState_
 190 {
 191   spx_uint32_t in_rate;
 192   spx_uint32_t out_rate;
 193   spx_uint32_t num_rate;
 194   spx_uint32_t den_rate;
 195
 196   int quality;
 197   spx_uint32_t nb_channels;
 198   spx_uint32_t filt_len;
 199   spx_uint32_t mem_alloc_size;
 200   spx_uint32_t buffer_size;
 201   int int_advance;
 202   int frac_advance;
 203   float cutoff;
 204   spx_uint32_t oversample;
 205   int initialised;
 206   int started;
 207   int use_full_sinc_table;
 208
 209   /* These are per-channel */
 210   spx_int32_t *last_sample;
 211   spx_uint32_t *samp_frac_num;
 212   spx_uint32_t *magic_samples;
 213
 214   spx_word16_t *mem;
 215   spx_word16_t *sinc_table;
 216   spx_uint32_t sinc_table_length;
 217   resampler_basic_func resampler_ptr;
 218
 219   int in_stride;
 220   int out_stride;
 221
 222   int use_sse:1;
 223   int use_sse2:1;
 224   int use_neon:1;
 225 };
 226
 227 static double kaiser12_table[68] = {
 228   0.99859849, 1.00000000, 0.99859849, 0.99440475, 0.98745105, 0.97779076,
 229   0.96549770, 0.95066529, 0.93340547, 0.91384741, 0.89213598, 0.86843014,
 230   0.84290116, 0.81573067, 0.78710866, 0.75723148, 0.72629970, 0.69451601,
 231   0.66208321, 0.62920216, 0.59606986, 0.56287762, 0.52980938, 0.49704014,
 232   0.46473455, 0.43304576, 0.40211431, 0.37206735, 0.34301800, 0.31506490,
 233   0.28829195, 0.26276832, 0.23854851, 0.21567274, 0.19416736, 0.17404546,
 234   0.15530766, 0.13794294, 0.12192957, 0.10723616, 0.09382272, 0.08164178,
 235   0.07063950, 0.06075685, 0.05193064, 0.04409466, 0.03718069, 0.03111947,
 236   0.02584161, 0.02127838, 0.01736250, 0.01402878, 0.01121463, 0.00886058,
 237   0.00691064, 0.00531256, 0.00401805, 0.00298291, 0.00216702, 0.00153438,
 238   0.00105297, 0.00069463, 0.00043489, 0.00025272, 0.00013031, 0.0000527734,
 239   0.00001000, 0.00000000
 240 };
 241
 242 /*
 243 static double kaiser12_table[36] = {
 244    0.99440475, 1.00000000, 0.99440475, 0.97779076, 0.95066529, 0.91384741,
 245    0.86843014, 0.81573067, 0.75723148, 0.69451601, 0.62920216, 0.56287762,
 246    0.49704014, 0.43304576, 0.37206735, 0.31506490, 0.26276832, 0.21567274,
 247    0.17404546, 0.13794294, 0.10723616, 0.08164178, 0.06075685, 0.04409466,
 248    0.03111947, 0.02127838, 0.01402878, 0.00886058, 0.00531256, 0.00298291,
 249    0.00153438, 0.00069463, 0.00025272, 0.0000527734, 0.00000500, 0.00000000};
 250 */
 251 static double kaiser10_table[36] = {
 252   0.99537781, 1.00000000, 0.99537781, 0.98162644, 0.95908712, 0.92831446,
 253   0.89005583, 0.84522401, 0.79486424, 0.74011713, 0.68217934, 0.62226347,
 254   0.56155915, 0.50119680, 0.44221549, 0.38553619, 0.33194107, 0.28205962,
 255   0.23636152, 0.19515633, 0.15859932, 0.12670280, 0.09935205, 0.07632451,
 256   0.05731132, 0.04193980, 0.02979584, 0.02044510, 0.01345224, 0.00839739,
 257   0.00488951, 0.00257636, 0.00115101, 0.00035515, 0.00000000, 0.00000000
 258 };
 259
 260 static double kaiser8_table[36] = {
 261   0.99635258, 1.00000000, 0.99635258, 0.98548012, 0.96759014, 0.94302200,
 262   0.91223751, 0.87580811, 0.83439927, 0.78875245, 0.73966538, 0.68797126,
 263   0.63451750, 0.58014482, 0.52566725, 0.47185369, 0.41941150, 0.36897272,
 264   0.32108304, 0.27619388, 0.23465776, 0.19672670, 0.16255380, 0.13219758,
 265   0.10562887, 0.08273982, 0.06335451, 0.04724088, 0.03412321, 0.02369490,
 266   0.01563093, 0.00959968, 0.00527363, 0.00233883, 0.00050000, 0.00000000
 267 };
 268
 269 static double kaiser6_table[36] = {
 270   0.99733006, 1.00000000, 0.99733006, 0.98935595, 0.97618418, 0.95799003,
 271   0.93501423, 0.90755855, 0.87598009, 0.84068475, 0.80211977, 0.76076565,
 272   0.71712752, 0.67172623, 0.62508937, 0.57774224, 0.53019925, 0.48295561,
 273   0.43647969, 0.39120616, 0.34752997, 0.30580127, 0.26632152, 0.22934058,
 274   0.19505503, 0.16360756, 0.13508755, 0.10953262, 0.08693120, 0.06722600,
 275   0.05031820, 0.03607231, 0.02432151, 0.01487334, 0.00752000, 0.00000000
 276 };
 277
 278 struct FuncDef
 279 {
 280   double *table;
 281   int oversample;
 282 };
 283
 284 static struct FuncDef _KAISER12 = { kaiser12_table, 64 };
 285
 286 #define KAISER12 (&_KAISER12)
 287 /*static struct FuncDef _KAISER12 = {kaiser12_table, 32};
 288 #define KAISER12 (&_KAISER12)*/
 289 static struct FuncDef _KAISER10 = { kaiser10_table, 32 };
 290
 291 #define KAISER10 (&_KAISER10)
 292 static struct FuncDef _KAISER8 = { kaiser8_table, 32 };
 293
 294 #define KAISER8 (&_KAISER8)
 295 static struct FuncDef _KAISER6 = { kaiser6_table, 32 };
 296
 297 #define KAISER6 (&_KAISER6)
 298
 299 struct QualityMapping
 300 {
 301   int base_length;
 302   int oversample;
 303   float downsample_bandwidth;
 304   float upsample_bandwidth;
 305   struct FuncDef *window_func;
 306 };
 307
 308
 309 /* This table maps conversion quality to internal parameters. There are two
 310    reasons that explain why the up-sampling bandwidth is larger than the
 311    down-sampling bandwidth:
 312    1) When up-sampling, we can assume that the spectrum is already attenuated
 313       close to the Nyquist rate (from an A/D or a previous resampling filter)
 314    2) Any aliasing that occurs very close to the Nyquist rate will be masked
 315       by the sinusoids/noise just below the Nyquist rate (guaranteed only for
 316       up-sampling).
 317 */
 318 static const struct QualityMapping quality_map[11] = {
 319   {8, 4, 0.830f, 0.860f, KAISER6},      /* Q0 */
 320   {16, 4, 0.850f, 0.880f, KAISER6},     /* Q1 */
 321   {32, 4, 0.882f, 0.910f, KAISER6},     /* Q2 *//* 82.3% cutoff ( ~60 dB stop) 6  */
 322   {48, 8, 0.895f, 0.917f, KAISER8},     /* Q3 *//* 84.9% cutoff ( ~80 dB stop) 8  */
 323   {64, 8, 0.921f, 0.940f, KAISER8},     /* Q4 *//* 88.7% cutoff ( ~80 dB stop) 8  */
 324   {80, 16, 0.922f, 0.940f, KAISER10},   /* Q5 *//* 89.1% cutoff (~100 dB stop) 10 */
 325   {96, 16, 0.940f, 0.945f, KAISER10},   /* Q6 *//* 91.5% cutoff (~100 dB stop) 10 */
 326   {128, 16, 0.950f, 0.950f, KAISER10},  /* Q7 *//* 93.1% cutoff (~100 dB stop) 10 */
 327   {160, 16, 0.960f, 0.960f, KAISER10},  /* Q8 *//* 94.5% cutoff (~100 dB stop) 10 */
 328   {192, 32, 0.968f, 0.968f, KAISER12},  /* Q9 *//* 95.5% cutoff (~100 dB stop) 10 */
 329   {256, 32, 0.975f, 0.975f, KAISER12},  /* Q10 *//* 96.6% cutoff (~100 dB stop) 10 */
 330 };
 331
 332 /*8,24,40,56,80,104,128,160,200,256,320*/
 333 #ifdef DOUBLE_PRECISION
 334 static double
 335 compute_func (double x, struct FuncDef *func)
 336 {
 337   double y, frac;
 338 #else
 339 static double
 340 compute_func (float x, struct FuncDef *func)
 341 {
 342   float y, frac;
 343 #endif
 344   double interp[4];
 345   int ind;
 346   y = x * func->oversample;
 347   ind = (int) floor (y);
 348   frac = (y - ind);
 349   /* CSE with handle the repeated powers */
 350   interp[3] = -0.1666666667 * frac + 0.1666666667 * (frac * frac * frac);
 351   interp[2] = frac + 0.5 * (frac * frac) - 0.5 * (frac * frac * frac);
 352   /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac; */
 353   interp[0] =
 354       -0.3333333333 * frac + 0.5 * (frac * frac) -
 355       0.1666666667 * (frac * frac * frac);
 356   /* Just to make sure we don't have rounding problems */
 357   interp[1] = 1.f - interp[3] - interp[2] - interp[0];
 358
 359   /*sum = frac*accum[1] + (1-frac)*accum[2]; */
 360   return interp[0] * func->table[ind] + interp[1] * func->table[ind + 1] +
 361       interp[2] * func->table[ind + 2] + interp[3] * func->table[ind + 3];
 362 }
 363
 364 #if 0
 365 #include <stdio.h>
 366 int
 367 main (int argc, char **argv)
 368 {
 369   int i;
 370   for (i = 0; i < 256; i++) {
 371     printf ("%f\n", compute_func (i / 256., KAISER12));
 372   }
 373   return 0;
 374 }
 375 #endif
 376
 377 #ifdef FIXED_POINT
 378 /* The slow way of computing a sinc for the table. Should improve that some day */
 379 static spx_word16_t
 380 sinc (float cutoff, float x, int N, struct FuncDef *window_func)
 381 {
 382   /*fprintf (stderr, "%f ", x); */
 383   float xx = x * cutoff;
 384   if (fabs (x) < 1e-6f)
 385     return WORD2INT (32768. * cutoff);
 386   else if (fabs (x) > .5f * N)
 387     return 0;
 388   /*FIXME: Can it really be any slower than this? */
 389   return WORD2INT (32768. * cutoff * sin (G_PI * xx) / (G_PI * xx) *
 390       compute_func (fabs (2. * x / N), window_func));
 391 }
 392 #else
 393 /* The slow way of computing a sinc for the table. Should improve that some day */
 394 #ifdef DOUBLE_PRECISION
 395 static spx_word16_t
 396 sinc (double cutoff, double x, int N, struct FuncDef *window_func)
 397 {
 398   /*fprintf (stderr, "%f ", x); */
 399   double xx = x * cutoff;
 400 #else
 401 static spx_word16_t
 402 sinc (float cutoff, float x, int N, struct FuncDef *window_func)
 403 {
 404   /*fprintf (stderr, "%f ", x); */
 405   float xx = x * cutoff;
 406 #endif
 407   if (fabs (x) < 1e-6)
 408     return cutoff;
 409   else if (fabs (x) > .5 * N)
 410     return 0;
 411   /*FIXME: Can it really be any slower than this? */
 412   return cutoff * sin (G_PI * xx) / (G_PI * xx) * compute_func (fabs (2. * x /
 413           N), window_func);
 414 }
 415 #endif
 416
 417 #ifdef FIXED_POINT
 418 static void
 419 cubic_coef (spx_word16_t x, spx_word16_t interp[4])
 420 {
 421   /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
 422      but I know it's MMSE-optimal on a sinc */
 423   spx_word16_t x2, x3;
 424   x2 = MULT16_16_P15 (x, x);
 425   x3 = MULT16_16_P15 (x, x2);
 426   interp[0] =
 427       PSHR32 (MULT16_16 (QCONST16 (-0.16667f, 15),
 428           x) + MULT16_16 (QCONST16 (0.16667f, 15), x3), 15);
 429   interp[1] =
 430       EXTRACT16 (EXTEND32 (x) + SHR32 (SUB32 (EXTEND32 (x2), EXTEND32 (x3)),
 431           1));
 432   interp[3] =
 433       PSHR32 (MULT16_16 (QCONST16 (-0.33333f, 15),
 434           x) + MULT16_16 (QCONST16 (.5f, 15),
 435           x2) - MULT16_16 (QCONST16 (0.16667f, 15), x3), 15);
 436   /* Just to make sure we don't have rounding problems */
 437   interp[2] = Q15_ONE - interp[0] - interp[1] - interp[3];
 438   if (interp[2] < 32767)
 439     interp[2] += 1;
 440 }
 441 #else
 442 static void
 443 cubic_coef (spx_word16_t frac, spx_word16_t interp[4])
 444 {
 445   /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation
 446      but I know it's MMSE-optimal on a sinc */
 447   interp[0] = -0.16667f * frac + 0.16667f * frac * frac * frac;
 448   interp[1] = frac + 0.5f * frac * frac - 0.5f * frac * frac * frac;
 449   /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac; */
 450   interp[3] =
 451       -0.33333f * frac + 0.5f * frac * frac - 0.16667f * frac * frac * frac;
 452   /* Just to make sure we don't have rounding problems */
 453   interp[2] = 1. - interp[0] - interp[1] - interp[3];
 454 }
 455 #endif
 456
 457 #ifndef DOUBLE_PRECISION
 458 static int
 459 resampler_basic_direct_single (SpeexResamplerState * st,
 460     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,
 461     spx_word16_t * out, spx_uint32_t * out_len)
 462 {
 463   const int N = st->filt_len;
 464   int out_sample = 0;
 465   int last_sample = st->last_sample[channel_index];
 466   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 467   const spx_word16_t *sinc_table = st->sinc_table;
 468   const int out_stride = st->out_stride;
 469   const int int_advance = st->int_advance;
 470   const int frac_advance = st->frac_advance;
 471   const spx_uint32_t den_rate = st->den_rate;
 472   spx_word32_t sum;
 473   int j;
 474
 475   while (!(last_sample >= (spx_int32_t) * in_len
 476           || out_sample >= (spx_int32_t) * out_len)) {
 477     const spx_word16_t *sinc = &sinc_table[samp_frac_num * N];
 478     const spx_word16_t *iptr = &in[last_sample];
 479
 480     SSE_FALLBACK (INNER_PRODUCT_SINGLE)
 481     NEON_FALLBACK (INNER_PRODUCT_SINGLE)
 482         sum = 0;
 483     for (j = 0; j < N; j++)
 484       sum += MULT16_16 (sinc[j], iptr[j]);
 485
 486 /*    This code is slower on most DSPs which have only 2 accumulators.
 487       Plus this forces truncation to 32 bits and you lose the HW guard bits.
 488       I think we can trust the compiler and let it vectorize and/or unroll itself.
 489       spx_word32_t accum[4] = {0,0,0,0};
 490       for(j=0;j<N;j+=4) {
 491         accum[0] += MULT16_16(sinc[j], iptr[j]);
 492         accum[1] += MULT16_16(sinc[j+1], iptr[j+1]);
 493         accum[2] += MULT16_16(sinc[j+2], iptr[j+2]);
 494         accum[3] += MULT16_16(sinc[j+3], iptr[j+3]);
 495       }
 496       sum = accum[0] + accum[1] + accum[2] + accum[3];
 497 */
 498 #if defined(OVERRIDE_INNER_PRODUCT_SINGLE) && defined(_USE_NEON)
 499     NEON_IMPLEMENTATION (INNER_PRODUCT_SINGLE)
 500     sum = inner_product_single (sinc, iptr, N);
 501     NEON_END(INNER_PRODUCT_SINGLE)
 502 #elif defined(OVERRIDE_INNER_PRODUCT_SINGLE) && defined(_USE_SSE)
 503     SSE_IMPLEMENTATION (INNER_PRODUCT_SINGLE)
 504         sum = inner_product_single (sinc, iptr, N);
 505     SSE_END (INNER_PRODUCT_SINGLE)
 506 #endif
 507     out[out_stride * out_sample++] = SATURATE32PSHR(sum, 15, 32767);
 508     last_sample += int_advance;
 509     samp_frac_num += frac_advance;
 510     if (samp_frac_num >= den_rate) {
 511       samp_frac_num -= den_rate;
 512       last_sample++;
 513     }
 514   }
 515
 516   st->last_sample[channel_index] = last_sample;
 517   st->samp_frac_num[channel_index] = samp_frac_num;
 518   return out_sample;
 519 }
 520 #endif
 521
 522 #ifdef FIXED_POINT
 523 #else
 524 /* This is the same as the previous function, except with a double-precision accumulator */
 525 static int
 526 resampler_basic_direct_double (SpeexResamplerState * st,
 527     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,
 528     spx_word16_t * out, spx_uint32_t * out_len)
 529 {
 530   const int N = st->filt_len;
 531   int out_sample = 0;
 532   int last_sample = st->last_sample[channel_index];
 533   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 534   const spx_word16_t *sinc_table = st->sinc_table;
 535   const int out_stride = st->out_stride;
 536   const int int_advance = st->int_advance;
 537   const int frac_advance = st->frac_advance;
 538   const spx_uint32_t den_rate = st->den_rate;
 539   double sum;
 540   int j;
 541
 542   while (!(last_sample >= (spx_int32_t) * in_len
 543           || out_sample >= (spx_int32_t) * out_len)) {
 544     const spx_word16_t *sinc = &sinc_table[samp_frac_num * N];
 545     const spx_word16_t *iptr = &in[last_sample];
 546
 547     SSE2_FALLBACK (INNER_PRODUCT_DOUBLE)
 548     double accum[4] = { 0, 0, 0, 0 };
 549
 550     for (j = 0; j < N; j += 4) {
 551       accum[0] += sinc[j] * iptr[j];
 552       accum[1] += sinc[j + 1] * iptr[j + 1];
 553       accum[2] += sinc[j + 2] * iptr[j + 2];
 554       accum[3] += sinc[j + 3] * iptr[j + 3];
 555     }
 556     sum = accum[0] + accum[1] + accum[2] + accum[3];
 557 #if defined(OVERRIDE_INNER_PRODUCT_DOUBLE) && defined(_USE_SSE2)
 558     SSE2_IMPLEMENTATION (INNER_PRODUCT_DOUBLE)
 559         sum = inner_product_double (sinc, iptr, N);
 560     SSE2_END (INNER_PRODUCT_DOUBLE)
 561 #endif
 562         out[out_stride * out_sample++] = PSHR32 (sum, 15);
 563     last_sample += int_advance;
 564     samp_frac_num += frac_advance;
 565     if (samp_frac_num >= den_rate) {
 566       samp_frac_num -= den_rate;
 567       last_sample++;
 568     }
 569   }
 570
 571   st->last_sample[channel_index] = last_sample;
 572   st->samp_frac_num[channel_index] = samp_frac_num;
 573   return out_sample;
 574 }
 575 #endif
 576
 577 #ifndef DOUBLE_PRECISION
 578 static int
 579 resampler_basic_interpolate_single (SpeexResamplerState * st,
 580     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,
 581     spx_word16_t * out, spx_uint32_t * out_len)
 582 {
 583   const int N = st->filt_len;
 584   int out_sample = 0;
 585   int last_sample = st->last_sample[channel_index];
 586   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 587   const int out_stride = st->out_stride;
 588   const int int_advance = st->int_advance;
 589   const int frac_advance = st->frac_advance;
 590   const spx_uint32_t den_rate = st->den_rate;
 591   int j;
 592   spx_word32_t sum;
 593
 594   while (!(last_sample >= (spx_int32_t) * in_len
 595           || out_sample >= (spx_int32_t) * out_len)) {
 596     const spx_word16_t *iptr = &in[last_sample];
 597
 598     const int offset = samp_frac_num * st->oversample / st->den_rate;
 599 #ifdef FIXED_POINT
 600     const spx_word16_t frac =
 601         PDIV32 (SHL32 ((samp_frac_num * st->oversample) % st->den_rate, 15),
 602         st->den_rate);
 603 #else
 604     const spx_word16_t frac =
 605         ((float) ((samp_frac_num * st->oversample) % st->den_rate)) /
 606         st->den_rate;
 607 #endif
 608     spx_word16_t interp[4];
 609
 610
 611     SSE_FALLBACK (INTERPOLATE_PRODUCT_SINGLE)
 612     spx_word32_t accum[4] = { 0, 0, 0, 0 };
 613
 614     for (j = 0; j < N; j++) {
 615       const spx_word16_t curr_in = iptr[j];
 616       accum[0] +=
 617           MULT16_16 (curr_in,
 618           st->sinc_table[4 + (j + 1) * st->oversample - offset - 2]);
 619       accum[1] +=
 620           MULT16_16 (curr_in,
 621           st->sinc_table[4 + (j + 1) * st->oversample - offset - 1]);
 622       accum[2] +=
 623           MULT16_16 (curr_in,
 624           st->sinc_table[4 + (j + 1) * st->oversample - offset]);
 625       accum[3] +=
 626           MULT16_16 (curr_in,
 627           st->sinc_table[4 + (j + 1) * st->oversample - offset + 1]);
 628     }
 629
 630     cubic_coef (frac, interp);
 631     sum =
 632         MULT16_32_Q15 (interp[0], SHR32 (accum[0],
 633             1)) + MULT16_32_Q15 (interp[1], SHR32 (accum[1],
 634             1)) + MULT16_32_Q15 (interp[2], SHR32 (accum[2],
 635             1)) + MULT16_32_Q15 (interp[3], SHR32 (accum[3], 1));
 636 #if defined(OVERRIDE_INTERPOLATE_PRODUCT_SINGLE) && defined(_USE_SSE)
 637     SSE_IMPLEMENTATION (INTERPOLATE_PRODUCT_SINGLE)
 638         cubic_coef (frac, interp);
 639     sum =
 640         interpolate_product_single (iptr,
 641         st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample,
 642         interp);
 643     SSE_END (INTERPOLATE_PRODUCT_SINGLE)
 644 #endif
 645     out[out_stride * out_sample++] = SATURATE32PSHR(sum, 14, 32767);
 646     last_sample += int_advance;
 647     samp_frac_num += frac_advance;
 648     if (samp_frac_num >= den_rate) {
 649       samp_frac_num -= den_rate;
 650       last_sample++;
 651     }
 652   }
 653
 654   st->last_sample[channel_index] = last_sample;
 655   st->samp_frac_num[channel_index] = samp_frac_num;
 656   return out_sample;
 657 }
 658 #endif
 659
 660 #ifdef FIXED_POINT
 661 #else
 662 /* This is the same as the previous function, except with a double-precision accumulator */
 663 static int
 664 resampler_basic_interpolate_double (SpeexResamplerState * st,
 665     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,
 666     spx_word16_t * out, spx_uint32_t * out_len)
 667 {
 668   const int N = st->filt_len;
 669   int out_sample = 0;
 670   int last_sample = st->last_sample[channel_index];
 671   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];
 672   const int out_stride = st->out_stride;
 673   const int int_advance = st->int_advance;
 674   const int frac_advance = st->frac_advance;
 675   const spx_uint32_t den_rate = st->den_rate;
 676   int j;
 677   spx_word32_t sum;
 678
 679   while (!(last_sample >= (spx_int32_t) * in_len
 680           || out_sample >= (spx_int32_t) * out_len)) {
 681     const spx_word16_t *iptr = &in[last_sample];
 682
 683     const int offset = samp_frac_num * st->oversample / st->den_rate;
 684 #ifdef FIXED_POINT
 685     const spx_word16_t frac =
 686         PDIV32 (SHL32 ((samp_frac_num * st->oversample) % st->den_rate, 15),
 687         st->den_rate);
 688 #else
 689 #ifdef DOUBLE_PRECISION
 690     const spx_word16_t frac =
 691         ((double) ((samp_frac_num * st->oversample) % st->den_rate)) /
 692         st->den_rate;
 693 #else
 694     const spx_word16_t frac =
 695         ((float) ((samp_frac_num * st->oversample) % st->den_rate)) /
 696         st->den_rate;
 697 #endif
 698 #endif
 699     spx_word16_t interp[4];
 700
 701
 702     SSE2_FALLBACK (INTERPOLATE_PRODUCT_DOUBLE)
 703     double accum[4] = { 0, 0, 0, 0 };
 704
 705     for (j = 0; j < N; j++) {
 706       const double curr_in = iptr[j];
 707       accum[0] +=
 708           MULT16_16 (curr_in,
 709           st->sinc_table[4 + (j + 1) * st->oversample - offset - 2]);
 710       accum[1] +=
 711           MULT16_16 (curr_in,
 712           st->sinc_table[4 + (j + 1) * st->oversample - offset - 1]);
 713       accum[2] +=
 714           MULT16_16 (curr_in,
 715           st->sinc_table[4 + (j + 1) * st->oversample - offset]);
 716       accum[3] +=
 717           MULT16_16 (curr_in,
 718           st->sinc_table[4 + (j + 1) * st->oversample - offset + 1]);
 719     }
 720
 721     cubic_coef (frac, interp);
 722     sum =
 723         MULT16_32_Q15 (interp[0], accum[0]) + MULT16_32_Q15 (interp[1],
 724         accum[1]) + MULT16_32_Q15 (interp[2],
 725         accum[2]) + MULT16_32_Q15 (interp[3], accum[3]);
 726 #if defined(OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE) && defined(_USE_SSE2)
 727     SSE2_IMPLEMENTATION (INTERPOLATE_PRODUCT_DOUBLE)
 728         cubic_coef (frac, interp);
 729     sum =
 730         interpolate_product_double (iptr,
 731         st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample,
 732         interp);
 733     SSE2_END (INTERPOLATE_PRODUCT_DOUBLE)
 734 #endif
 735         out[out_stride * out_sample++] = PSHR32 (sum, 15);
 736     last_sample += int_advance;
 737     samp_frac_num += frac_advance;
 738     if (samp_frac_num >= den_rate) {
 739       samp_frac_num -= den_rate;
 740       last_sample++;
 741     }
 742   }
 743
 744   st->last_sample[channel_index] = last_sample;
 745   st->samp_frac_num[channel_index] = samp_frac_num;
 746   return out_sample;
 747 }
 748 #endif
 749
 750 static void
 751 update_filter (SpeexResamplerState * st)
 752 {
 753   spx_uint32_t old_length;
 754
 755   old_length = st->filt_len;
 756   st->oversample = quality_map[st->quality].oversample;
 757   st->filt_len = quality_map[st->quality].base_length;
 758
 759   if (st->num_rate > st->den_rate) {
 760     /* down-sampling */
 761     st->cutoff =
 762         quality_map[st->quality].downsample_bandwidth * st->den_rate /
 763         st->num_rate;
 764     /* FIXME: divide the numerator and denominator by a certain amount if they're too large */
 765     st->filt_len = st->filt_len * st->num_rate / st->den_rate;
 766     /* Round down to make sure we have a multiple of 4 */
 767     st->filt_len &= (~0x3);
 768     if (2 * st->den_rate < st->num_rate)
 769       st->oversample >>= 1;
 770     if (4 * st->den_rate < st->num_rate)
 771       st->oversample >>= 1;
 772     if (8 * st->den_rate < st->num_rate)
 773       st->oversample >>= 1;
 774     if (16 * st->den_rate < st->num_rate)
 775       st->oversample >>= 1;
 776     if (st->oversample < 1)
 777       st->oversample = 1;
 778   } else {
 779     /* up-sampling */
 780     st->cutoff = quality_map[st->quality].upsample_bandwidth;
 781   }
 782
 783   /* Choose the resampling type that requires the least amount of memory */
 784   /* Or if the full sinc table is explicitely requested, use that */
 785   if (st->use_full_sinc_table || (st->den_rate <= st->oversample)) {
 786     spx_uint32_t i;
 787     if (!st->sinc_table)
 788       st->sinc_table =
 789           (spx_word16_t *) speex_alloc (st->filt_len * st->den_rate *
 790           sizeof (spx_word16_t));
 791     else if (st->sinc_table_length < st->filt_len * st->den_rate) {
 792       st->sinc_table =
 793           (spx_word16_t *) speex_realloc (st->sinc_table,
 794           st->filt_len * st->den_rate * sizeof (spx_word16_t));
 795       st->sinc_table_length = st->filt_len * st->den_rate;
 796     }
 797     for (i = 0; i < st->den_rate; i++) {
 798       spx_int32_t j;
 799       for (j = 0; j < st->filt_len; j++) {
 800         st->sinc_table[i * st->filt_len + j] =
 801             sinc (st->cutoff, ((j - (spx_int32_t) st->filt_len / 2 + 1) -
 802 #ifdef DOUBLE_PRECISION
 803                 ((double) i) / st->den_rate), st->filt_len,
 804 #else
 805                 ((float) i) / st->den_rate), st->filt_len,
 806 #endif
 807             quality_map[st->quality].window_func);
 808       }
 809     }
 810 #ifdef FIXED_POINT
 811     st->resampler_ptr = resampler_basic_direct_single;
 812 #else
 813 #ifdef DOUBLE_PRECISION
 814     st->resampler_ptr = resampler_basic_direct_double;
 815 #else
 816     if (st->quality > 8)
 817       st->resampler_ptr = resampler_basic_direct_double;
 818     else
 819       st->resampler_ptr = resampler_basic_direct_single;
 820 #endif
 821 #endif
 822     /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff); */
 823   } else {
 824     spx_int32_t i;
 825     if (!st->sinc_table)
 826       st->sinc_table =
 827           (spx_word16_t *) speex_alloc ((st->filt_len * st->oversample +
 828               8) * sizeof (spx_word16_t));
 829     else if (st->sinc_table_length < st->filt_len * st->oversample + 8) {
 830       st->sinc_table =
 831           (spx_word16_t *) speex_realloc (st->sinc_table,
 832           (st->filt_len * st->oversample + 8) * sizeof (spx_word16_t));
 833       st->sinc_table_length = st->filt_len * st->oversample + 8;
 834     }
 835     for (i = -4; i < (spx_int32_t) (st->oversample * st->filt_len + 4); i++)
 836       st->sinc_table[i + 4] =
 837 #ifdef DOUBLE_PRECISION
 838           sinc (st->cutoff, (i / (double) st->oversample - st->filt_len / 2),
 839 #else
 840           sinc (st->cutoff, (i / (float) st->oversample - st->filt_len / 2),
 841 #endif
 842           st->filt_len, quality_map[st->quality].window_func);
 843 #ifdef FIXED_POINT
 844     st->resampler_ptr = resampler_basic_interpolate_single;
 845 #else
 846 #ifdef DOUBLE_PRECISION
 847     st->resampler_ptr = resampler_basic_interpolate_double;
 848 #else
 849     if (st->quality > 8)
 850       st->resampler_ptr = resampler_basic_interpolate_double;
 851     else
 852       st->resampler_ptr = resampler_basic_interpolate_single;
 853 #endif
 854 #endif
 855     /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff); */
 856   }
 857   st->int_advance = st->num_rate / st->den_rate;
 858   st->frac_advance = st->num_rate % st->den_rate;
 859
 860
 861   /* Here's the place where we update the filter memory to take into account
 862      the change in filter length. It's probably the messiest part of the code
 863      due to handling of lots of corner cases. */
 864   if (!st->mem) {
 865     spx_uint32_t i;
 866     st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;
 867     st->mem =
 868         (spx_word16_t *) speex_alloc (st->nb_channels * st->mem_alloc_size *
 869         sizeof (spx_word16_t));
 870     for (i = 0; i < st->nb_channels * st->mem_alloc_size; i++)
 871       st->mem[i] = 0;
 872     /*speex_warning("init filter"); */
 873   } else if (!st->started) {
 874     spx_uint32_t i;
 875     st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;
 876     st->mem =
 877         (spx_word16_t *) speex_realloc (st->mem,
 878         st->nb_channels * st->mem_alloc_size * sizeof (spx_word16_t));
 879     for (i = 0; i < st->nb_channels * st->mem_alloc_size; i++)
 880       st->mem[i] = 0;
 881     /*speex_warning("reinit filter"); */
 882   } else if (st->filt_len > old_length) {
 883     spx_int32_t i;
 884     /* Increase the filter length */
 885     /*speex_warning("increase filter size"); */
 886     int old_alloc_size = st->mem_alloc_size;
 887     if ((st->filt_len - 1 + st->buffer_size) > st->mem_alloc_size) {
 888       st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;
 889       st->mem =
 890           (spx_word16_t *) speex_realloc (st->mem,
 891           st->nb_channels * st->mem_alloc_size * sizeof (spx_word16_t));
 892     }
 893     for (i = st->nb_channels - 1; i >= 0; i--) {
 894       spx_int32_t j;
 895       spx_uint32_t olen = old_length;
 896       /*if (st->magic_samples[i]) */
 897       {
 898         /* Try and remove the magic samples as if nothing had happened */
 899
 900         /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */
 901         olen = old_length + 2 * st->magic_samples[i];
 902         for (j = old_length - 2 + st->magic_samples[i]; j >= 0; j--)
 903           st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]] =
 904               st->mem[i * old_alloc_size + j];
 905         for (j = 0; j < st->magic_samples[i]; j++)
 906           st->mem[i * st->mem_alloc_size + j] = 0;
 907         st->magic_samples[i] = 0;
 908       }
 909       if (st->filt_len > olen) {
 910         /* If the new filter length is still bigger than the "augmented" length */
 911         /* Copy data going backward */
 912         for (j = 0; j < olen - 1; j++)
 913           st->mem[i * st->mem_alloc_size + (st->filt_len - 2 - j)] =
 914               st->mem[i * st->mem_alloc_size + (olen - 2 - j)];
 915         /* Then put zeros for lack of anything better */
 916         for (; j < st->filt_len - 1; j++)
 917           st->mem[i * st->mem_alloc_size + (st->filt_len - 2 - j)] = 0;
 918         /* Adjust last_sample */
 919         st->last_sample[i] += (st->filt_len - olen) / 2;
 920       } else {
 921         /* Put back some of the magic! */
 922         st->magic_samples[i] = (olen - st->filt_len) / 2;
 923         for (j = 0; j < st->filt_len - 1 + st->magic_samples[i]; j++)
 924           st->mem[i * st->mem_alloc_size + j] =
 925               st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]];
 926       }
 927     }
 928   } else if (st->filt_len < old_length) {
 929     spx_uint32_t i;
 930     /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"
 931        samples so they can be used directly as input the next time(s) */
 932     for (i = 0; i < st->nb_channels; i++) {
 933       spx_uint32_t j;
 934       spx_uint32_t old_magic = st->magic_samples[i];
 935       st->magic_samples[i] = (old_length - st->filt_len) / 2;
 936       /* We must copy some of the memory that's no longer used */
 937       /* Copy data going backward */
 938       for (j = 0; j < st->filt_len - 1 + st->magic_samples[i] + old_magic; j++)
 939         st->mem[i * st->mem_alloc_size + j] =
 940             st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]];
 941       st->magic_samples[i] += old_magic;
 942     }
 943   }
 944
 945 }
 946
 947 EXPORT SpeexResamplerState *
 948 speex_resampler_init (spx_uint32_t nb_channels, spx_uint32_t in_rate,
 949     spx_uint32_t out_rate, int quality, SpeexResamplerSincFilterMode sinc_filter_mode,
 950     spx_uint32_t sinc_filter_auto_threshold, int *err)
 951 {
 952   return speex_resampler_init_frac (nb_channels, in_rate, out_rate, in_rate,
 953       out_rate, quality, sinc_filter_mode, sinc_filter_auto_threshold, err);
 954 }
 955
 956 #if defined HAVE_ORC && !defined DISABLE_ORC
 957 static void
 958 check_insn_set (SpeexResamplerState * st, const char *name)
 959 {
 960   if (!name)
 961     return;
 962 #ifdef _USE_SSE
 963   if (!strcmp (name, "sse"))
 964     st->use_sse = 1;
 965 #endif
 966 #ifdef _USE_SSE2
 967   if (!strcmp (name, "sse2"))
 968     st->use_sse = st->use_sse2 = 1;
 969 #endif
 970 #ifdef _USE_NEON
 971   if (!strcmp (name, "neon"))
 972     st->use_neon = 1;
 973 #endif
 974 }
 975 #endif
 976
 977 EXPORT SpeexResamplerState *
 978 speex_resampler_init_frac (spx_uint32_t nb_channels, spx_uint32_t ratio_num,
 979     spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate,
 980     int quality, SpeexResamplerSincFilterMode sinc_filter_mode,
 981     spx_uint32_t sinc_filter_auto_threshold, int *err)
 982 {
 983   spx_uint32_t i;
 984   SpeexResamplerState *st;
 985   int use_full_sinc_table = 0;
 986   if (quality > 10 || quality < 0) {
 987     if (err)
 988       *err = RESAMPLER_ERR_INVALID_ARG;
 989     return NULL;
 990   }
 991   switch (sinc_filter_mode) {
 992     case RESAMPLER_SINC_FILTER_INTERPOLATED:
 993       use_full_sinc_table = 0;
 994       break;
 995     case RESAMPLER_SINC_FILTER_FULL:
 996       use_full_sinc_table = 1;
 997       break;
 998     case RESAMPLER_SINC_FILTER_AUTO:
 999       /* Handled below */
1000       break;
1001     default:
1002       if (err)
1003         *err = RESAMPLER_ERR_INVALID_ARG;
1004       return NULL;
1005   }
1006
1007   st = (SpeexResamplerState *) speex_alloc (sizeof (SpeexResamplerState));
1008   st->initialised = 0;
1009   st->started = 0;
1010   st->in_rate = 0;
1011   st->out_rate = 0;
1012   st->num_rate = 0;
1013   st->den_rate = 0;
1014   st->quality = -1;
1015   st->sinc_table_length = 0;
1016   st->mem_alloc_size = 0;
1017   st->filt_len = 0;
1018   st->mem = 0;
1019   st->resampler_ptr = 0;
1020   st->use_full_sinc_table = use_full_sinc_table;
1021
1022   st->cutoff = 1.f;
1023   st->nb_channels = nb_channels;
1024   st->in_stride = 1;
1025   st->out_stride = 1;
1026
1027 #ifdef FIXED_POINT
1028   st->buffer_size = 160;
1029 #else
1030   st->buffer_size = 160;
1031 #endif
1032
1033   st->use_sse = st->use_sse2 = 0;
1034   st->use_neon = 0;
1035 #if defined HAVE_ORC && !defined DISABLE_ORC
1036   orc_init ();
1037   {
1038     OrcTarget *target = orc_target_get_default ();
1039     if (target) {
1040       unsigned int flags = orc_target_get_default_flags (target);
1041       check_insn_set (st, orc_target_get_name (target));
1042       for (i = 0; i < 32; ++i) {
1043         if (flags & (1 << i)) {
1044           check_insn_set (st, orc_target_get_flag_name (target, i));
1045         }
1046       }
1047     }
1048   }
1049 #endif
1050
1051   /* Per channel data */
1052   st->last_sample = (spx_int32_t *) speex_alloc (nb_channels * sizeof (int));
1053   st->magic_samples = (spx_uint32_t *) speex_alloc (nb_channels * sizeof (int));
1054   st->samp_frac_num = (spx_uint32_t *) speex_alloc (nb_channels * sizeof (int));
1055   for (i = 0; i < nb_channels; i++) {
1056     st->last_sample[i] = 0;
1057     st->magic_samples[i] = 0;
1058     st->samp_frac_num[i] = 0;
1059   }
1060
1061   speex_resampler_set_quality (st, quality);
1062   speex_resampler_set_rate_frac (st, ratio_num, ratio_den, in_rate, out_rate);
1063
1064   if (sinc_filter_mode == RESAMPLER_SINC_FILTER_AUTO) {
1065     /*
1066     Estimate how big the filter table would become if the full mode were to be used
1067     calculations used correspond to the ones in update_filter()
1068     if the size is bigger than the threshold, use interpolated sinc instead
1069     */
1070     spx_uint32_t base_filter_length = st->filt_len = quality_map[st->quality].base_length;
1071     spx_uint32_t filter_table_size = base_filter_length * st->den_rate * sizeof(spx_uint16_t);
1072     st->use_full_sinc_table = (filter_table_size > sinc_filter_auto_threshold) ? 0 : 1;
1073   }
1074
1075   update_filter (st);
1076
1077   st->initialised = 1;
1078   if (err)
1079     *err = RESAMPLER_ERR_SUCCESS;
1080
1081   return st;
1082 }
1083
1084 EXPORT void
1085 speex_resampler_destroy (SpeexResamplerState * st)
1086 {
1087   speex_free (st->mem);
1088   speex_free (st->sinc_table);
1089   speex_free (st->last_sample);
1090   speex_free (st->magic_samples);
1091   speex_free (st->samp_frac_num);
1092   speex_free (st);
1093 }
1094
1095 static int
1096 speex_resampler_process_native (SpeexResamplerState * st,
1097     spx_uint32_t channel_index, spx_uint32_t * in_len, spx_word16_t * out,
1098     spx_uint32_t * out_len)
1099 {
1100   int j = 0;
1101   const int N = st->filt_len;
1102   int out_sample = 0;
1103   spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
1104   spx_uint32_t ilen;
1105
1106   st->started = 1;
1107
1108   /* Call the right resampler through the function ptr */
1109   out_sample = st->resampler_ptr (st, channel_index, mem, in_len, out, out_len);
1110
1111   if (st->last_sample[channel_index] < (spx_int32_t) * in_len)
1112     *in_len = st->last_sample[channel_index];
1113   *out_len = out_sample;
1114   st->last_sample[channel_index] -= *in_len;
1115
1116   ilen = *in_len;
1117
1118   for (j = 0; j < N - 1; ++j)
1119     mem[j] = mem[j + ilen];
1120
1121   return RESAMPLER_ERR_SUCCESS;
1122 }
1123
1124 static int
1125 speex_resampler_magic (SpeexResamplerState * st, spx_uint32_t channel_index,
1126     spx_word16_t ** out, spx_uint32_t out_len)
1127 {
1128   spx_uint32_t tmp_in_len = st->magic_samples[channel_index];
1129   spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;
1130   const int N = st->filt_len;
1131
1132   speex_resampler_process_native (st, channel_index, &tmp_in_len, *out,
1133       &out_len);
1134
1135   st->magic_samples[channel_index] -= tmp_in_len;
1136
1137   /* If we couldn't process all "magic" input samples, save the rest for next time */
1138   if (st->magic_samples[channel_index]) {
1139     spx_uint32_t i;
1140     for (i = 0; i < st->magic_samples[channel_index]; i++)
1141       mem[N - 1 + i] = mem[N - 1 + i + tmp_in_len];
1142   }
1143   *out += out_len * st->out_stride;
1144   return out_len;
1145 }
1146
1147 #ifdef FIXED_POINT
1148 EXPORT int
1149 speex_resampler_process_int (SpeexResamplerState * st,
1150     spx_uint32_t channel_index, const spx_int16_t * in, spx_uint32_t * in_len,
1151     spx_int16_t * out, spx_uint32_t * out_len)
1152 #else
1153 #ifdef DOUBLE_PRECISION
1154 EXPORT int
1155 speex_resampler_process_float (SpeexResamplerState * st,
1156     spx_uint32_t channel_index, const double *in, spx_uint32_t * in_len,
1157     double *out, spx_uint32_t * out_len)
1158 #else
1159 EXPORT int
1160 speex_resampler_process_float (SpeexResamplerState * st,
1161     spx_uint32_t channel_index, const float *in, spx_uint32_t * in_len,
1162     float *out, spx_uint32_t * out_len)
1163 #endif
1164 #endif
1165 {
1166   int j;
1167   spx_uint32_t ilen = *in_len;
1168   spx_uint32_t olen = *out_len;
1169   spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
1170   const int filt_offs = st->filt_len - 1;
1171   const spx_uint32_t xlen = st->mem_alloc_size - filt_offs;
1172   const int istride = st->in_stride;
1173
1174   if (st->magic_samples[channel_index])
1175     olen -= speex_resampler_magic (st, channel_index, &out, olen);
1176   if (!st->magic_samples[channel_index]) {
1177     while (ilen && olen) {
1178       spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
1179       spx_uint32_t ochunk = olen;
1180
1181       if (in) {
1182         for (j = 0; j < ichunk; ++j)
1183           x[j + filt_offs] = in[j * istride];
1184       } else {
1185         for (j = 0; j < ichunk; ++j)
1186           x[j + filt_offs] = 0;
1187       }
1188       speex_resampler_process_native (st, channel_index, &ichunk, out, &ochunk);
1189       ilen -= ichunk;
1190       olen -= ochunk;
1191       out += ochunk * st->out_stride;
1192       if (in)
1193         in += ichunk * istride;
1194     }
1195   }
1196   *in_len -= ilen;
1197   *out_len -= olen;
1198   return RESAMPLER_ERR_SUCCESS;
1199 }
1200
1201 #ifdef FIXED_POINT
1202 EXPORT int
1203 speex_resampler_process_float (SpeexResamplerState * st,
1204     spx_uint32_t channel_index, const float *in, spx_uint32_t * in_len,
1205     float *out, spx_uint32_t * out_len)
1206 #else
1207 EXPORT int
1208 speex_resampler_process_int (SpeexResamplerState * st,
1209     spx_uint32_t channel_index, const spx_int16_t * in, spx_uint32_t * in_len,
1210     spx_int16_t * out, spx_uint32_t * out_len)
1211 #endif
1212 {
1213   int j;
1214   const int istride_save = st->in_stride;
1215   const int ostride_save = st->out_stride;
1216   spx_uint32_t ilen = *in_len;
1217   spx_uint32_t olen = *out_len;
1218   spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;
1219   const spx_uint32_t xlen = st->mem_alloc_size - (st->filt_len - 1);
1220 #ifdef VAR_ARRAYS
1221   const unsigned int ylen =
1222       (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;
1223   VARDECL (spx_word16_t * ystack);
1224   ALLOC (ystack, ylen, spx_word16_t);
1225 #else
1226   const unsigned int ylen = FIXED_STACK_ALLOC;
1227   spx_word16_t ystack[FIXED_STACK_ALLOC];
1228 #endif
1229
1230   st->out_stride = 1;
1231
1232   while (ilen && olen) {
1233     spx_word16_t *y = ystack;
1234     spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;
1235     spx_uint32_t ochunk = (olen > ylen) ? ylen : olen;
1236     spx_uint32_t omagic = 0;
1237
1238     if (st->magic_samples[channel_index]) {
1239       omagic = speex_resampler_magic (st, channel_index, &y, ochunk);
1240       ochunk -= omagic;
1241       olen -= omagic;
1242     }
1243     if (!st->magic_samples[channel_index]) {
1244       if (in) {
1245         for (j = 0; j < ichunk; ++j)
1246 #ifdef FIXED_POINT
1247           x[j + st->filt_len - 1] = WORD2INT (in[j * istride_save]);
1248 #else
1249           x[j + st->filt_len - 1] = in[j * istride_save];
1250 #endif
1251       } else {
1252         for (j = 0; j < ichunk; ++j)
1253           x[j + st->filt_len - 1] = 0;
1254       }
1255
1256       speex_resampler_process_native (st, channel_index, &ichunk, y, &ochunk);
1257     } else {
1258       ichunk = 0;
1259       ochunk = 0;
1260     }
1261
1262     for (j = 0; j < ochunk + omagic; ++j)
1263 #ifdef FIXED_POINT
1264       out[j * ostride_save] = ystack[j];
1265 #else
1266       out[j * ostride_save] = WORD2INT (ystack[j]);
1267 #endif
1268
1269     ilen -= ichunk;
1270     olen -= ochunk;
1271     out += (ochunk + omagic) * ostride_save;
1272     if (in)
1273       in += ichunk * istride_save;
1274   }
1275   st->out_stride = ostride_save;
1276   *in_len -= ilen;
1277   *out_len -= olen;
1278
1279   return RESAMPLER_ERR_SUCCESS;
1280 }
1281
1282 #ifdef DOUBLE_PRECISION
1283 EXPORT int
1284 speex_resampler_process_interleaved_float (SpeexResamplerState * st,
1285     const double *in, spx_uint32_t * in_len, double *out,
1286     spx_uint32_t * out_len)
1287 #else
1288 EXPORT int
1289 speex_resampler_process_interleaved_float (SpeexResamplerState * st,
1290     const float *in, spx_uint32_t * in_len, float *out, spx_uint32_t * out_len)
1291 #endif
1292 {
1293   spx_uint32_t i;
1294   int istride_save, ostride_save;
1295   spx_uint32_t bak_len = *out_len;
1296   istride_save = st->in_stride;
1297   ostride_save = st->out_stride;
1298   st->in_stride = st->out_stride = st->nb_channels;
1299   for (i = 0; i < st->nb_channels; i++) {
1300     *out_len = bak_len;
1301     if (in != NULL)
1302       speex_resampler_process_float (st, i, in + i, in_len, out + i, out_len);
1303     else
1304       speex_resampler_process_float (st, i, NULL, in_len, out + i, out_len);
1305   }
1306   st->in_stride = istride_save;
1307   st->out_stride = ostride_save;
1308   return RESAMPLER_ERR_SUCCESS;
1309 }
1310
1311 EXPORT int
1312 speex_resampler_process_interleaved_int (SpeexResamplerState * st,
1313     const spx_int16_t * in, spx_uint32_t * in_len, spx_int16_t * out,
1314     spx_uint32_t * out_len)
1315 {
1316   spx_uint32_t i;
1317   int istride_save, ostride_save;
1318   spx_uint32_t bak_len = *out_len;
1319   istride_save = st->in_stride;
1320   ostride_save = st->out_stride;
1321   st->in_stride = st->out_stride = st->nb_channels;
1322   for (i = 0; i < st->nb_channels; i++) {
1323     *out_len = bak_len;
1324     if (in != NULL)
1325       speex_resampler_process_int (st, i, in + i, in_len, out + i, out_len);
1326     else
1327       speex_resampler_process_int (st, i, NULL, in_len, out + i, out_len);
1328   }
1329   st->in_stride = istride_save;
1330   st->out_stride = ostride_save;
1331   return RESAMPLER_ERR_SUCCESS;
1332 }
1333
1334 EXPORT int
1335 speex_resampler_set_rate (SpeexResamplerState * st, spx_uint32_t in_rate,
1336     spx_uint32_t out_rate)
1337 {
1338   return speex_resampler_set_rate_frac (st, in_rate, out_rate, in_rate,
1339       out_rate);
1340 }
1341
1342 EXPORT void
1343 speex_resampler_get_rate (SpeexResamplerState * st, spx_uint32_t * in_rate,
1344     spx_uint32_t * out_rate)
1345 {
1346   *in_rate = st->in_rate;
1347   *out_rate = st->out_rate;
1348 }
1349
1350 EXPORT int
1351 speex_resampler_set_rate_frac (SpeexResamplerState * st, spx_uint32_t ratio_num,
1352     spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate)
1353 {
1354   spx_uint32_t fact;
1355   spx_uint32_t old_den;
1356   spx_uint32_t i;
1357   if (st->in_rate == in_rate && st->out_rate == out_rate
1358       && st->num_rate == ratio_num && st->den_rate == ratio_den)
1359     return RESAMPLER_ERR_SUCCESS;
1360
1361   old_den = st->den_rate;
1362   st->in_rate = in_rate;
1363   st->out_rate = out_rate;
1364   st->num_rate = ratio_num;
1365   st->den_rate = ratio_den;
1366   /* FIXME: This is terribly inefficient, but who cares (at least for now)? */
1367   for (fact = 2; fact <= IMIN (st->num_rate, st->den_rate); fact++) {
1368     while ((st->num_rate % fact == 0) && (st->den_rate % fact == 0)) {
1369       st->num_rate /= fact;
1370       st->den_rate /= fact;
1371     }
1372   }
1373
1374   if (old_den > 0) {
1375     for (i = 0; i < st->nb_channels; i++) {
1376       st->samp_frac_num[i] = st->samp_frac_num[i] * st->den_rate / old_den;
1377       /* Safety net */
1378       if (st->samp_frac_num[i] >= st->den_rate)
1379         st->samp_frac_num[i] = st->den_rate - 1;
1380     }
1381   }
1382
1383   if (st->initialised)
1384     update_filter (st);
1385   return RESAMPLER_ERR_SUCCESS;
1386 }
1387
1388 EXPORT void
1389 speex_resampler_get_ratio (SpeexResamplerState * st, spx_uint32_t * ratio_num,
1390     spx_uint32_t * ratio_den)
1391 {
1392   *ratio_num = st->num_rate;
1393   *ratio_den = st->den_rate;
1394 }
1395
1396 EXPORT int
1397 speex_resampler_set_quality (SpeexResamplerState * st, int quality)
1398 {
1399   if (quality > 10 || quality < 0)
1400     return RESAMPLER_ERR_INVALID_ARG;
1401   if (st->quality == quality)
1402     return RESAMPLER_ERR_SUCCESS;
1403   st->quality = quality;
1404   if (st->initialised)
1405     update_filter (st);
1406   return RESAMPLER_ERR_SUCCESS;
1407 }
1408
1409 EXPORT void
1410 speex_resampler_get_quality (SpeexResamplerState * st, int *quality)
1411 {
1412   *quality = st->quality;
1413 }
1414
1415 EXPORT void
1416 speex_resampler_set_input_stride (SpeexResamplerState * st, spx_uint32_t stride)
1417 {
1418   st->in_stride = stride;
1419 }
1420
1421 EXPORT void
1422 speex_resampler_get_input_stride (SpeexResamplerState * st,
1423     spx_uint32_t * stride)
1424 {
1425   *stride = st->in_stride;
1426 }
1427
1428 EXPORT void
1429 speex_resampler_set_output_stride (SpeexResamplerState * st,
1430     spx_uint32_t stride)
1431 {
1432   st->out_stride = stride;
1433 }
1434
1435 EXPORT void
1436 speex_resampler_get_output_stride (SpeexResamplerState * st,
1437     spx_uint32_t * stride)
1438 {
1439   *stride = st->out_stride;
1440 }
1441
1442 EXPORT int
1443 speex_resampler_get_input_latency (SpeexResamplerState * st)
1444 {
1445   return st->filt_len / 2;
1446 }
1447
1448 EXPORT int
1449 speex_resampler_get_output_latency (SpeexResamplerState * st)
1450 {
1451   return ((st->filt_len / 2) * st->den_rate +
1452       (st->num_rate >> 1)) / st->num_rate;
1453 }
1454
1455 EXPORT int
1456 speex_resampler_get_filt_len (SpeexResamplerState * st)
1457 {
1458   return st->filt_len;
1459 }
1460
1461 EXPORT int
1462 speex_resampler_get_sinc_filter_mode (SpeexResamplerState * st)
1463 {
1464   return st->use_full_sinc_table;
1465 }
1466
1467 EXPORT int
1468 speex_resampler_skip_zeros (SpeexResamplerState * st)
1469 {
1470   spx_uint32_t i;
1471   for (i = 0; i < st->nb_channels; i++)
1472     st->last_sample[i] = st->filt_len / 2;
1473   return RESAMPLER_ERR_SUCCESS;
1474 }
1475
1476 EXPORT int
1477 speex_resampler_reset_mem (SpeexResamplerState * st)
1478 {
1479   spx_uint32_t i;
1480   for (i = 0; i < st->nb_channels * (st->filt_len - 1); i++)
1481     st->mem[i] = 0;
1482   return RESAMPLER_ERR_SUCCESS;
1483 }
1484
1485 EXPORT const char *
1486 speex_resampler_strerror (int err)
1487 {
1488   switch (err) {
1489     case RESAMPLER_ERR_SUCCESS:
1490       return "Success.";
1491     case RESAMPLER_ERR_ALLOC_FAILED:
1492       return "Memory allocation failed.";
1493     case RESAMPLER_ERR_BAD_STATE:
1494       return "Bad resampler state.";
1495     case RESAMPLER_ERR_INVALID_ARG:
1496       return "Invalid argument.";
1497     case RESAMPLER_ERR_PTR_OVERLAP:
1498       return "Input and output buffers overlap.";
1499     default:
1500       return "Unknown error. Bad error code or strange version mismatch.";
1501   }
1502 }