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