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 #include <glib.h>

    68 #include <glibconfig.h>

    69 #include <e32def.h>

    70 

    71 #ifndef __SYMBIAN32__

    72 #define EXPORT EXPORT_C

    73 #else

    74 #define EXPORT G_GNUC_INTERNAL

    75 #endif

    76 

    77 static inline void *

    78 speex_alloc (int size)

    79 {

    80   return g_malloc0 (size);

    81 }

    82 

    83 static inline void *

    84 speex_realloc (void *ptr, int size)

    85 {

    86   return g_realloc (ptr, size);

    87 }

    88 

    89 static inline void

    90 speex_free (void *ptr)

    91 {

    92   g_free (ptr);

    93 }

    94 

    95 #include "speex_resampler.h"

    96 #include "arch_float.h"

    97 #else /* OUTSIDE_SPEEX */

    98 

    99 #include "speex_resampler.h"

   100 #include "arch.h"

   101 //#include "os_support.h"

   102 #endif /* OUTSIDE_SPEEX */

   103 

   104 #include <math.h>

   105 

   106 #ifndef M_PI

   107 #define M_PI 3.14159263

   108 #endif

   109 

   110 #ifdef FIXED_POINT

   111 #define WORD2INT(x) ((x) < -32767 ? -32768 : ((x) > 32766 ? 32767 : (x)))

   112 #else

   113 #define WORD2INT(x) ((x) < -32767.5f ? -32768 : ((x) > 32766.5f ? 32767 : floor(.5+(x))))

   114 #endif

   115 

   116 #define IMAX(a,b) ((a) > (b) ? (a) : (b))

   117 #define IMIN(a,b) ((a) < (b) ? (a) : (b))

   118 

   119 #ifndef NULL

   120 #define NULL 0

   121 #endif

   122 

   123 #ifdef _USE_SSE

   124 #include "resample_sse.h"

   125 #endif

   126 

   127 /* Numer of elements to allocate on the stack */

   128 #ifdef VAR_ARRAYS

   129 #define FIXED_STACK_ALLOC 8192

   130 #else

   131 #define FIXED_STACK_ALLOC 1024

   132 #endif

   133 

   134 typedef int (*resampler_basic_func) (SpeexResamplerState *, spx_uint32_t,

   135     const spx_word16_t *, spx_uint32_t *, spx_word16_t *, spx_uint32_t *);

   136 

   137 struct SpeexResamplerState_

   138 {

   139   spx_uint32_t in_rate;

   140   spx_uint32_t out_rate;

   141   spx_uint32_t num_rate;

   142   spx_uint32_t den_rate;

   143 

   144   int quality;

   145   spx_uint32_t nb_channels;

   146   spx_uint32_t filt_len;

   147   spx_uint32_t mem_alloc_size;

   148   spx_uint32_t buffer_size;

   149   int int_advance;

   150   int frac_advance;

   151   float cutoff;

   152   spx_uint32_t oversample;

   153   int initialised;

   154   int started;

   155 

   156   /* These are per-channel */

   157   spx_int32_t *last_sample;

   158   spx_uint32_t *samp_frac_num;

   159   spx_uint32_t *magic_samples;

   160 

   161   spx_word16_t *mem;

   162   spx_word16_t *sinc_table;

   163   spx_uint32_t sinc_table_length;

   164   resampler_basic_func resampler_ptr;

   165 

   166   int in_stride;

   167   int out_stride;

   168 };

   169 

   170 static double kaiser12_table[68] = {

   171   0.99859849, 1.00000000, 0.99859849, 0.99440475, 0.98745105, 0.97779076,

   172   0.96549770, 0.95066529, 0.93340547, 0.91384741, 0.89213598, 0.86843014,

   173   0.84290116, 0.81573067, 0.78710866, 0.75723148, 0.72629970, 0.69451601,

   174   0.66208321, 0.62920216, 0.59606986, 0.56287762, 0.52980938, 0.49704014,

   175   0.46473455, 0.43304576, 0.40211431, 0.37206735, 0.34301800, 0.31506490,

   176   0.28829195, 0.26276832, 0.23854851, 0.21567274, 0.19416736, 0.17404546,

   177   0.15530766, 0.13794294, 0.12192957, 0.10723616, 0.09382272, 0.08164178,

   178   0.07063950, 0.06075685, 0.05193064, 0.04409466, 0.03718069, 0.03111947,

   179   0.02584161, 0.02127838, 0.01736250, 0.01402878, 0.01121463, 0.00886058,

   180   0.00691064, 0.00531256, 0.00401805, 0.00298291, 0.00216702, 0.00153438,

   181   0.00105297, 0.00069463, 0.00043489, 0.00025272, 0.00013031, 0.0000527734,

   182   0.00001000, 0.00000000

   183 };

   184 

   185 /*

   186 static double kaiser12_table[36] = {

   187    0.99440475, 1.00000000, 0.99440475, 0.97779076, 0.95066529, 0.91384741,

   188    0.86843014, 0.81573067, 0.75723148, 0.69451601, 0.62920216, 0.56287762,

   189    0.49704014, 0.43304576, 0.37206735, 0.31506490, 0.26276832, 0.21567274,

   190    0.17404546, 0.13794294, 0.10723616, 0.08164178, 0.06075685, 0.04409466,

   191    0.03111947, 0.02127838, 0.01402878, 0.00886058, 0.00531256, 0.00298291,

   192    0.00153438, 0.00069463, 0.00025272, 0.0000527734, 0.00000500, 0.00000000};

   193 */

   194 static double kaiser10_table[36] = {

   195   0.99537781, 1.00000000, 0.99537781, 0.98162644, 0.95908712, 0.92831446,

   196   0.89005583, 0.84522401, 0.79486424, 0.74011713, 0.68217934, 0.62226347,

   197   0.56155915, 0.50119680, 0.44221549, 0.38553619, 0.33194107, 0.28205962,

   198   0.23636152, 0.19515633, 0.15859932, 0.12670280, 0.09935205, 0.07632451,

   199   0.05731132, 0.04193980, 0.02979584, 0.02044510, 0.01345224, 0.00839739,

   200   0.00488951, 0.00257636, 0.00115101, 0.00035515, 0.00000000, 0.00000000

   201 };

   202 

   203 static double kaiser8_table[36] = {

   204   0.99635258, 1.00000000, 0.99635258, 0.98548012, 0.96759014, 0.94302200,

   205   0.91223751, 0.87580811, 0.83439927, 0.78875245, 0.73966538, 0.68797126,

   206   0.63451750, 0.58014482, 0.52566725, 0.47185369, 0.41941150, 0.36897272,

   207   0.32108304, 0.27619388, 0.23465776, 0.19672670, 0.16255380, 0.13219758,

   208   0.10562887, 0.08273982, 0.06335451, 0.04724088, 0.03412321, 0.02369490,

   209   0.01563093, 0.00959968, 0.00527363, 0.00233883, 0.00050000, 0.00000000

   210 };

   211 

   212 static double kaiser6_table[36] = {

   213   0.99733006, 1.00000000, 0.99733006, 0.98935595, 0.97618418, 0.95799003,

   214   0.93501423, 0.90755855, 0.87598009, 0.84068475, 0.80211977, 0.76076565,

   215   0.71712752, 0.67172623, 0.62508937, 0.57774224, 0.53019925, 0.48295561,

   216   0.43647969, 0.39120616, 0.34752997, 0.30580127, 0.26632152, 0.22934058,

   217   0.19505503, 0.16360756, 0.13508755, 0.10953262, 0.08693120, 0.06722600,

   218   0.05031820, 0.03607231, 0.02432151, 0.01487334, 0.00752000, 0.00000000

   219 };

   220 

   221 struct FuncDef

   222 {

   223   double *table;

   224   int oversample;

   225 };

   226 

   227 static struct FuncDef _KAISER12 = { kaiser12_table, 64 };

   228 

   229 #define KAISER12 (&_KAISER12)

   230 /*static struct FuncDef _KAISER12 = {kaiser12_table, 32};

   231 #define KAISER12 (&_KAISER12)*/

   232 static struct FuncDef _KAISER10 = { kaiser10_table, 32 };

   233 

   234 #define KAISER10 (&_KAISER10)

   235 static struct FuncDef _KAISER8 = { kaiser8_table, 32 };

   236 

   237 #define KAISER8 (&_KAISER8)

   238 static struct FuncDef _KAISER6 = { kaiser6_table, 32 };

   239 

   240 #define KAISER6 (&_KAISER6)

   241 

   242 struct QualityMapping

   243 {

   244   int base_length;

   245   int oversample;

   246   float downsample_bandwidth;

   247   float upsample_bandwidth;

   248   struct FuncDef *window_func;

   249 };

   250 

   251 

   252 /* This table maps conversion quality to internal parameters. There are two

   253    reasons that explain why the up-sampling bandwidth is larger than the 

   254    down-sampling bandwidth:

   255    1) When up-sampling, we can assume that the spectrum is already attenuated

   256       close to the Nyquist rate (from an A/D or a previous resampling filter)

   257    2) Any aliasing that occurs very close to the Nyquist rate will be masked

   258       by the sinusoids/noise just below the Nyquist rate (guaranteed only for

   259       up-sampling).

   260 */

   261 static const struct QualityMapping quality_map[11] = {

   262   {8, 4, 0.830f, 0.860f, KAISER6},      /* Q0 */

   263   {16, 4, 0.850f, 0.880f, KAISER6},     /* Q1 */

   264   {32, 4, 0.882f, 0.910f, KAISER6},     /* Q2 *//* 82.3% cutoff ( ~60 dB stop) 6  */

   265   {48, 8, 0.895f, 0.917f, KAISER8},     /* Q3 *//* 84.9% cutoff ( ~80 dB stop) 8  */

   266   {64, 8, 0.921f, 0.940f, KAISER8},     /* Q4 *//* 88.7% cutoff ( ~80 dB stop) 8  */

   267   {80, 16, 0.922f, 0.940f, KAISER10},   /* Q5 *//* 89.1% cutoff (~100 dB stop) 10 */

   268   {96, 16, 0.940f, 0.945f, KAISER10},   /* Q6 *//* 91.5% cutoff (~100 dB stop) 10 */

   269   {128, 16, 0.950f, 0.950f, KAISER10},  /* Q7 *//* 93.1% cutoff (~100 dB stop) 10 */

   270   {160, 16, 0.960f, 0.960f, KAISER10},  /* Q8 *//* 94.5% cutoff (~100 dB stop) 10 */

   271   {192, 32, 0.968f, 0.968f, KAISER12},  /* Q9 *//* 95.5% cutoff (~100 dB stop) 10 */

   272   {256, 32, 0.975f, 0.975f, KAISER12},  /* Q10 *//* 96.6% cutoff (~100 dB stop) 10 */

   273 };

   274 

   275 /*8,24,40,56,80,104,128,160,200,256,320*/

   276 #ifdef DOUBLE_PRECISION

   277 static double

   278 compute_func (double x, struct FuncDef *func)

   279 {

   280   double y, frac;

   281 #else

   282 static double

   283 compute_func (float x, struct FuncDef *func)

   284 {

   285   float y, frac;

   286 #endif

   287   double interp[4];

   288   int ind;

   289   y = x * func->oversample;

   290   ind = (int) floor (y);

   291   frac = (y - ind);

   292   /* CSE with handle the repeated powers */

   293   interp[3] = -0.1666666667 * frac + 0.1666666667 * (frac * frac * frac);

   294   interp[2] = frac + 0.5 * (frac * frac) - 0.5 * (frac * frac * frac);

   295   /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac; */

   296   interp[0] =

   297       -0.3333333333 * frac + 0.5 * (frac * frac) -

   298       0.1666666667 * (frac * frac * frac);

   299   /* Just to make sure we don't have rounding problems */

   300   interp[1] = 1.f - interp[3] - interp[2] - interp[0];

   301 

   302   /*sum = frac*accum[1] + (1-frac)*accum[2]; */

   303   return interp[0] * func->table[ind] + interp[1] * func->table[ind + 1] +

   304       interp[2] * func->table[ind + 2] + interp[3] * func->table[ind + 3];

   305 }

   306 

   307 #if 0

   308 #include <stdio.h>

   309 int

   310 main (int argc, char **argv)

   311 {

   312   int i;

   313   for (i = 0; i < 256; i++) {

   314     printf ("%f\n", compute_func (i / 256., KAISER12));

   315   }

   316   return 0;

   317 }

   318 #endif

   319 

   320 #ifdef FIXED_POINT

   321 /* The slow way of computing a sinc for the table. Should improve that some day */

   322 static spx_word16_t

   323 sinc (float cutoff, float x, int N, struct FuncDef *window_func)

   324 {

   325   /*fprintf (stderr, "%f ", x); */

   326   float xx = x * cutoff;

   327   if (fabs (x) < 1e-6f)

   328     return WORD2INT (32768. * cutoff);

   329   else if (fabs (x) > .5f * N)

   330     return 0;

   331   /*FIXME: Can it really be any slower than this? */

   332   return WORD2INT (32768. * cutoff * sin (M_PI * xx) / (M_PI * xx) *

   333       compute_func (fabs (2. * x / N), window_func));

   334 }

   335 #else

   336 /* The slow way of computing a sinc for the table. Should improve that some day */

   337 #ifdef DOUBLE_PRECISION

   338 static spx_word16_t

   339 sinc (double cutoff, double x, int N, struct FuncDef *window_func)

   340 {

   341   /*fprintf (stderr, "%f ", x); */

   342   double xx = x * cutoff;

   343 #else

   344 static spx_word16_t

   345 sinc (float cutoff, float x, int N, struct FuncDef *window_func)

   346 {

   347   /*fprintf (stderr, "%f ", x); */

   348   float xx = x * cutoff;

   349 #endif

   350   if (fabs (x) < 1e-6)

   351     return cutoff;

   352   else if (fabs (x) > .5 * N)

   353     return 0;

   354   /*FIXME: Can it really be any slower than this? */

   355   return cutoff * sin (M_PI * xx) / (M_PI * xx) * compute_func (fabs (2. * x /

   356           N), window_func);

   357 }

   358 #endif

   359 

   360 #ifdef FIXED_POINT

   361 static void

   362 cubic_coef (spx_word16_t x, spx_word16_t interp[4])

   363 {

   364   /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation

   365      but I know it's MMSE-optimal on a sinc */

   366   spx_word16_t x2, x3;

   367   x2 = MULT16_16_P15 (x, x);

   368   x3 = MULT16_16_P15 (x, x2);

   369   interp[0] =

   370       PSHR32 (MULT16_16 (QCONST16 (-0.16667f, 15),

   371           x) + MULT16_16 (QCONST16 (0.16667f, 15), x3), 15);

   372   interp[1] =

   373       EXTRACT16 (EXTEND32 (x) + SHR32 (SUB32 (EXTEND32 (x2), EXTEND32 (x3)),

   374           1));

   375   interp[3] =

   376       PSHR32 (MULT16_16 (QCONST16 (-0.33333f, 15),

   377           x) + MULT16_16 (QCONST16 (.5f, 15),

   378           x2) - MULT16_16 (QCONST16 (0.16667f, 15), x3), 15);

   379   /* Just to make sure we don't have rounding problems */

   380   interp[2] = Q15_ONE - interp[0] - interp[1] - interp[3];

   381   if (interp[2] < 32767)

   382     interp[2] += 1;

   383 }

   384 #else

   385 static void

   386 cubic_coef (spx_word16_t frac, spx_word16_t interp[4])

   387 {

   388   /* Compute interpolation coefficients. I'm not sure whether this corresponds to cubic interpolation

   389      but I know it's MMSE-optimal on a sinc */

   390   interp[0] = -0.16667f * frac + 0.16667f * frac * frac * frac;

   391   interp[1] = frac + 0.5f * frac * frac - 0.5f * frac * frac * frac;

   392   /*interp[2] = 1.f - 0.5f*frac - frac*frac + 0.5f*frac*frac*frac; */

   393   interp[3] =

   394       -0.33333f * frac + 0.5f * frac * frac - 0.16667f * frac * frac * frac;

   395   /* Just to make sure we don't have rounding problems */

   396   interp[2] = 1. - interp[0] - interp[1] - interp[3];

   397 }

   398 #endif

   399 

   400 #ifndef DOUBLE_PRECISION

   401 static int

   402 resampler_basic_direct_single (SpeexResamplerState * st,

   403     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,

   404     spx_word16_t * out, spx_uint32_t * out_len)

   405 {

   406   const int N = st->filt_len;

   407   int out_sample = 0;

   408   int last_sample = st->last_sample[channel_index];

   409   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];

   410   const spx_word16_t *sinc_table = st->sinc_table;

   411   const int out_stride = st->out_stride;

   412   const int int_advance = st->int_advance;

   413   const int frac_advance = st->frac_advance;

   414   const spx_uint32_t den_rate = st->den_rate;

   415   spx_word32_t sum;

   416   int j;

   417 

   418   while (!(last_sample >= (spx_int32_t) * in_len

   419           || out_sample >= (spx_int32_t) * out_len)) {

   420     const spx_word16_t *sinc = &sinc_table[samp_frac_num * N];

   421     const spx_word16_t *iptr = &in[last_sample];

   422 

   423 #ifndef OVERRIDE_INNER_PRODUCT_SINGLE

   424     float accum[4] = { 0, 0, 0, 0 };

   425 

   426     for (j = 0; j < N; j += 4) {

   427       accum[0] += sinc[j] * iptr[j];

   428       accum[1] += sinc[j + 1] * iptr[j + 1];

   429       accum[2] += sinc[j + 2] * iptr[j + 2];

   430       accum[3] += sinc[j + 3] * iptr[j + 3];

   431     }

   432     sum = accum[0] + accum[1] + accum[2] + accum[3];

   433 #else

   434     sum = inner_product_single (sinc, iptr, N);

   435 #endif

   436 

   437     out[out_stride * out_sample++] = PSHR32 (sum, 15);

   438     last_sample += int_advance;

   439     samp_frac_num += frac_advance;

   440     if (samp_frac_num >= den_rate) {

   441       samp_frac_num -= den_rate;

   442       last_sample++;

   443     }

   444   }

   445 

   446   st->last_sample[channel_index] = last_sample;

   447   st->samp_frac_num[channel_index] = samp_frac_num;

   448   return out_sample;

   449 }

   450 #endif

   451 

   452 #ifdef FIXED_POINT

   453 #else

   454 /* This is the same as the previous function, except with a double-precision accumulator */

   455 static int

   456 resampler_basic_direct_double (SpeexResamplerState * st,

   457     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,

   458     spx_word16_t * out, spx_uint32_t * out_len)

   459 {

   460   const int N = st->filt_len;

   461   int out_sample = 0;

   462   int last_sample = st->last_sample[channel_index];

   463   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];

   464   const spx_word16_t *sinc_table = st->sinc_table;

   465   const int out_stride = st->out_stride;

   466   const int int_advance = st->int_advance;

   467   const int frac_advance = st->frac_advance;

   468   const spx_uint32_t den_rate = st->den_rate;

   469   double sum;

   470   int j;

   471 

   472   while (!(last_sample >= (spx_int32_t) * in_len

   473           || out_sample >= (spx_int32_t) * out_len)) {

   474     const spx_word16_t *sinc = &sinc_table[samp_frac_num * N];

   475     const spx_word16_t *iptr = &in[last_sample];

   476 

   477 #ifndef OVERRIDE_INNER_PRODUCT_DOUBLE

   478     double accum[4] = { 0, 0, 0, 0 };

   479 

   480     for (j = 0; j < N; j += 4) {

   481       accum[0] += sinc[j] * iptr[j];

   482       accum[1] += sinc[j + 1] * iptr[j + 1];

   483       accum[2] += sinc[j + 2] * iptr[j + 2];

   484       accum[3] += sinc[j + 3] * iptr[j + 3];

   485     }

   486     sum = accum[0] + accum[1] + accum[2] + accum[3];

   487 #else

   488     sum = inner_product_double (sinc, iptr, N);

   489 #endif

   490 

   491     out[out_stride * out_sample++] = PSHR32 (sum, 15);

   492     last_sample += int_advance;

   493     samp_frac_num += frac_advance;

   494     if (samp_frac_num >= den_rate) {

   495       samp_frac_num -= den_rate;

   496       last_sample++;

   497     }

   498   }

   499 

   500   st->last_sample[channel_index] = last_sample;

   501   st->samp_frac_num[channel_index] = samp_frac_num;

   502   return out_sample;

   503 }

   504 #endif

   505 

   506 #ifndef DOUBLE_PRECISION

   507 static int

   508 resampler_basic_interpolate_single (SpeexResamplerState * st,

   509     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,

   510     spx_word16_t * out, spx_uint32_t * out_len)

   511 {

   512   const int N = st->filt_len;

   513   int out_sample = 0;

   514   int last_sample = st->last_sample[channel_index];

   515   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];

   516   const int out_stride = st->out_stride;

   517   const int int_advance = st->int_advance;

   518   const int frac_advance = st->frac_advance;

   519   const spx_uint32_t den_rate = st->den_rate;

   520   int j;

   521   spx_word32_t sum;

   522 

   523   while (!(last_sample >= (spx_int32_t) * in_len

   524           || out_sample >= (spx_int32_t) * out_len)) {

   525     const spx_word16_t *iptr = &in[last_sample];

   526 

   527     const int offset = samp_frac_num * st->oversample / st->den_rate;

   528 #ifdef FIXED_POINT

   529     const spx_word16_t frac =

   530         PDIV32 (SHL32 ((samp_frac_num * st->oversample) % st->den_rate, 15),

   531         st->den_rate);

   532 #else

   533     const spx_word16_t frac =

   534         ((float) ((samp_frac_num * st->oversample) % st->den_rate)) /

   535         st->den_rate;

   536 #endif

   537     spx_word16_t interp[4];

   538 

   539 

   540 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_SINGLE

   541     spx_word32_t accum[4] = { 0, 0, 0, 0 };

   542 

   543     for (j = 0; j < N; j++) {

   544       const spx_word16_t curr_in = iptr[j];

   545       accum[0] +=

   546           MULT16_16 (curr_in,

   547           st->sinc_table[4 + (j + 1) * st->oversample - offset - 2]);

   548       accum[1] +=

   549           MULT16_16 (curr_in,

   550           st->sinc_table[4 + (j + 1) * st->oversample - offset - 1]);

   551       accum[2] +=

   552           MULT16_16 (curr_in,

   553           st->sinc_table[4 + (j + 1) * st->oversample - offset]);

   554       accum[3] +=

   555           MULT16_16 (curr_in,

   556           st->sinc_table[4 + (j + 1) * st->oversample - offset + 1]);

   557     }

   558 

   559     cubic_coef (frac, interp);

   560     sum =

   561         MULT16_32_Q15 (interp[0], accum[0]) + MULT16_32_Q15 (interp[1],

   562         accum[1]) + MULT16_32_Q15 (interp[2],

   563         accum[2]) + MULT16_32_Q15 (interp[3], accum[3]);

   564 #else

   565     cubic_coef (frac, interp);

   566     sum =

   567         interpolate_product_single (iptr,

   568         st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample,

   569         interp);

   570 #endif

   571 

   572     out[out_stride * out_sample++] = PSHR32 (sum, 15);

   573     last_sample += int_advance;

   574     samp_frac_num += frac_advance;

   575     if (samp_frac_num >= den_rate) {

   576       samp_frac_num -= den_rate;

   577       last_sample++;

   578     }

   579   }

   580 

   581   st->last_sample[channel_index] = last_sample;

   582   st->samp_frac_num[channel_index] = samp_frac_num;

   583   return out_sample;

   584 }

   585 #endif

   586 

   587 #ifdef FIXED_POINT

   588 #else

   589 /* This is the same as the previous function, except with a double-precision accumulator */

   590 static int

   591 resampler_basic_interpolate_double (SpeexResamplerState * st,

   592     spx_uint32_t channel_index, const spx_word16_t * in, spx_uint32_t * in_len,

   593     spx_word16_t * out, spx_uint32_t * out_len)

   594 {

   595   const int N = st->filt_len;

   596   int out_sample = 0;

   597   int last_sample = st->last_sample[channel_index];

   598   spx_uint32_t samp_frac_num = st->samp_frac_num[channel_index];

   599   const int out_stride = st->out_stride;

   600   const int int_advance = st->int_advance;

   601   const int frac_advance = st->frac_advance;

   602   const spx_uint32_t den_rate = st->den_rate;

   603   int j;

   604   spx_word32_t sum;

   605 

   606   while (!(last_sample >= (spx_int32_t) * in_len

   607           || out_sample >= (spx_int32_t) * out_len)) {

   608     const spx_word16_t *iptr = &in[last_sample];

   609 

   610     const int offset = samp_frac_num * st->oversample / st->den_rate;

   611 #ifdef FIXED_POINT

   612     const spx_word16_t frac =

   613         PDIV32 (SHL32 ((samp_frac_num * st->oversample) % st->den_rate, 15),

   614         st->den_rate);

   615 #else

   616 #ifdef DOUBLE_PRECISION

   617     const spx_word16_t frac =

   618         ((double) ((samp_frac_num * st->oversample) % st->den_rate)) /

   619         st->den_rate;

   620 #else

   621     const spx_word16_t frac =

   622         ((float) ((samp_frac_num * st->oversample) % st->den_rate)) /

   623         st->den_rate;

   624 #endif

   625 #endif

   626     spx_word16_t interp[4];

   627 

   628 

   629 #ifndef OVERRIDE_INTERPOLATE_PRODUCT_DOUBLE

   630     double accum[4] = { 0, 0, 0, 0 };

   631 

   632     for (j = 0; j < N; j++) {

   633       const double curr_in = iptr[j];

   634       accum[0] +=

   635           MULT16_16 (curr_in,

   636           st->sinc_table[4 + (j + 1) * st->oversample - offset - 2]);

   637       accum[1] +=

   638           MULT16_16 (curr_in,

   639           st->sinc_table[4 + (j + 1) * st->oversample - offset - 1]);

   640       accum[2] +=

   641           MULT16_16 (curr_in,

   642           st->sinc_table[4 + (j + 1) * st->oversample - offset]);

   643       accum[3] +=

   644           MULT16_16 (curr_in,

   645           st->sinc_table[4 + (j + 1) * st->oversample - offset + 1]);

   646     }

   647 

   648     cubic_coef (frac, interp);

   649     sum =

   650         MULT16_32_Q15 (interp[0], accum[0]) + MULT16_32_Q15 (interp[1],

   651         accum[1]) + MULT16_32_Q15 (interp[2],

   652         accum[2]) + MULT16_32_Q15 (interp[3], accum[3]);

   653 #else

   654     cubic_coef (frac, interp);

   655     sum =

   656         interpolate_product_double (iptr,

   657         st->sinc_table + st->oversample + 4 - offset - 2, N, st->oversample,

   658         interp);

   659 #endif

   660 

   661     out[out_stride * out_sample++] = PSHR32 (sum, 15);

   662     last_sample += int_advance;

   663     samp_frac_num += frac_advance;

   664     if (samp_frac_num >= den_rate) {

   665       samp_frac_num -= den_rate;

   666       last_sample++;

   667     }

   668   }

   669 

   670   st->last_sample[channel_index] = last_sample;

   671   st->samp_frac_num[channel_index] = samp_frac_num;

   672   return out_sample;

   673 }

   674 #endif

   675 

   676 static void

   677 update_filter (SpeexResamplerState * st)

   678 {

   679   spx_uint32_t old_length;

   680 

   681   old_length = st->filt_len;

   682   st->oversample = quality_map[st->quality].oversample;

   683   st->filt_len = quality_map[st->quality].base_length;

   684 

   685   if (st->num_rate > st->den_rate) {

   686     /* down-sampling */

   687     st->cutoff =

   688         quality_map[st->quality].downsample_bandwidth * st->den_rate /

   689         st->num_rate;

   690     /* FIXME: divide the numerator and denominator by a certain amount if they're too large */

   691     st->filt_len = st->filt_len * st->num_rate / st->den_rate;

   692     /* Round down to make sure we have a multiple of 4 */

   693     st->filt_len &= (~0x3);

   694     if (2 * st->den_rate < st->num_rate)

   695       st->oversample >>= 1;

   696     if (4 * st->den_rate < st->num_rate)

   697       st->oversample >>= 1;

   698     if (8 * st->den_rate < st->num_rate)

   699       st->oversample >>= 1;

   700     if (16 * st->den_rate < st->num_rate)

   701       st->oversample >>= 1;

   702     if (st->oversample < 1)

   703       st->oversample = 1;

   704   } else {

   705     /* up-sampling */

   706     st->cutoff = quality_map[st->quality].upsample_bandwidth;

   707   }

   708 

   709   /* Choose the resampling type that requires the least amount of memory */

   710   if (st->den_rate <= st->oversample) {

   711     spx_uint32_t i;

   712     if (!st->sinc_table)

   713       st->sinc_table =

   714           (spx_word16_t *) speex_alloc (st->filt_len * st->den_rate *

   715           sizeof (spx_word16_t));

   716     else if (st->sinc_table_length < st->filt_len * st->den_rate) {

   717       st->sinc_table =

   718           (spx_word16_t *) speex_realloc (st->sinc_table,

   719           st->filt_len * st->den_rate * sizeof (spx_word16_t));

   720       st->sinc_table_length = st->filt_len * st->den_rate;

   721     }

   722     for (i = 0; i < st->den_rate; i++) {

   723       spx_int32_t j;

   724       for (j = 0; j < st->filt_len; j++) {

   725         st->sinc_table[i * st->filt_len + j] =

   726             sinc (st->cutoff, ((j - (spx_int32_t) st->filt_len / 2 + 1) -

   727 #ifdef DOUBLE_PRECISION

   728                 ((double) i) / st->den_rate), st->filt_len,

   729 #else

   730                 ((float) i) / st->den_rate), st->filt_len,

   731 #endif

   732             quality_map[st->quality].window_func);

   733       }

   734     }

   735 #ifdef FIXED_POINT

   736     st->resampler_ptr = resampler_basic_direct_single;

   737 #else

   738 #ifdef DOUBLE_PRECISION

   739     st->resampler_ptr = resampler_basic_direct_double;

   740 #else

   741     if (st->quality > 8)

   742       st->resampler_ptr = resampler_basic_direct_double;

   743     else

   744       st->resampler_ptr = resampler_basic_direct_single;

   745 #endif

   746 #endif

   747     /*fprintf (stderr, "resampler uses direct sinc table and normalised cutoff %f\n", cutoff); */

   748   } else {

   749     spx_int32_t i;

   750     if (!st->sinc_table)

   751       st->sinc_table =

   752           (spx_word16_t *) speex_alloc ((st->filt_len * st->oversample +

   753               8) * sizeof (spx_word16_t));

   754     else if (st->sinc_table_length < st->filt_len * st->oversample + 8) {

   755       st->sinc_table =

   756           (spx_word16_t *) speex_realloc (st->sinc_table,

   757           (st->filt_len * st->oversample + 8) * sizeof (spx_word16_t));

   758       st->sinc_table_length = st->filt_len * st->oversample + 8;

   759     }

   760     for (i = -4; i < (spx_int32_t) (st->oversample * st->filt_len + 4); i++)

   761       st->sinc_table[i + 4] =

   762 #ifdef DOUBLE_PRECISION

   763           sinc (st->cutoff, (i / (double) st->oversample - st->filt_len / 2),

   764 #else

   765           sinc (st->cutoff, (i / (float) st->oversample - st->filt_len / 2),

   766 #endif

   767           st->filt_len, quality_map[st->quality].window_func);

   768 #ifdef FIXED_POINT

   769     st->resampler_ptr = resampler_basic_interpolate_single;

   770 #else

   771 #ifdef DOUBLE_PRECISION

   772     st->resampler_ptr = resampler_basic_interpolate_double;

   773 #else

   774     if (st->quality > 8)

   775       st->resampler_ptr = resampler_basic_interpolate_double;

   776     else

   777       st->resampler_ptr = resampler_basic_interpolate_single;

   778 #endif

   779 #endif

   780     /*fprintf (stderr, "resampler uses interpolated sinc table and normalised cutoff %f\n", cutoff); */

   781   }

   782   st->int_advance = st->num_rate / st->den_rate;

   783   st->frac_advance = st->num_rate % st->den_rate;

   784 

   785 

   786   /* Here's the place where we update the filter memory to take into account

   787      the change in filter length. It's probably the messiest part of the code

   788      due to handling of lots of corner cases. */

   789   if (!st->mem) {

   790     spx_uint32_t i;

   791     st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;

   792     st->mem =

   793         (spx_word16_t *) speex_alloc (st->nb_channels * st->mem_alloc_size *

   794         sizeof (spx_word16_t));

   795     for (i = 0; i < st->nb_channels * st->mem_alloc_size; i++)

   796       st->mem[i] = 0;

   797     /*speex_warning("init filter"); */

   798   } else if (!st->started) {

   799     spx_uint32_t i;

   800     st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;

   801     st->mem =

   802         (spx_word16_t *) speex_realloc (st->mem,

   803         st->nb_channels * st->mem_alloc_size * sizeof (spx_word16_t));

   804     for (i = 0; i < st->nb_channels * st->mem_alloc_size; i++)

   805       st->mem[i] = 0;

   806     /*speex_warning("reinit filter"); */

   807   } else if (st->filt_len > old_length) {

   808     spx_int32_t i;

   809     /* Increase the filter length */

   810     /*speex_warning("increase filter size"); */

   811     int old_alloc_size = st->mem_alloc_size;

   812     if ((st->filt_len - 1 + st->buffer_size) > st->mem_alloc_size) {

   813       st->mem_alloc_size = st->filt_len - 1 + st->buffer_size;

   814       st->mem =

   815           (spx_word16_t *) speex_realloc (st->mem,

   816           st->nb_channels * st->mem_alloc_size * sizeof (spx_word16_t));

   817     }

   818     for (i = st->nb_channels - 1; i >= 0; i--) {

   819       spx_int32_t j;

   820       spx_uint32_t olen = old_length;

   821       /*if (st->magic_samples[i]) */

   822       {

   823         /* Try and remove the magic samples as if nothing had happened */

   824 

   825         /* FIXME: This is wrong but for now we need it to avoid going over the array bounds */

   826         olen = old_length + 2 * st->magic_samples[i];

   827         for (j = old_length - 2 + st->magic_samples[i]; j >= 0; j--)

   828           st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]] =

   829               st->mem[i * old_alloc_size + j];

   830         for (j = 0; j < st->magic_samples[i]; j++)

   831           st->mem[i * st->mem_alloc_size + j] = 0;

   832         st->magic_samples[i] = 0;

   833       }

   834       if (st->filt_len > olen) {

   835         /* If the new filter length is still bigger than the "augmented" length */

   836         /* Copy data going backward */

   837         for (j = 0; j < olen - 1; j++)

   838           st->mem[i * st->mem_alloc_size + (st->filt_len - 2 - j)] =

   839               st->mem[i * st->mem_alloc_size + (olen - 2 - j)];

   840         /* Then put zeros for lack of anything better */

   841         for (; j < st->filt_len - 1; j++)

   842           st->mem[i * st->mem_alloc_size + (st->filt_len - 2 - j)] = 0;

   843         /* Adjust last_sample */

   844         st->last_sample[i] += (st->filt_len - olen) / 2;

   845       } else {

   846         /* Put back some of the magic! */

   847         st->magic_samples[i] = (olen - st->filt_len) / 2;

   848         for (j = 0; j < st->filt_len - 1 + st->magic_samples[i]; j++)

   849           st->mem[i * st->mem_alloc_size + j] =

   850               st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]];

   851       }

   852     }

   853   } else if (st->filt_len < old_length) {

   854     spx_uint32_t i;

   855     /* Reduce filter length, this a bit tricky. We need to store some of the memory as "magic"

   856        samples so they can be used directly as input the next time(s) */

   857     for (i = 0; i < st->nb_channels; i++) {

   858       spx_uint32_t j;

   859       spx_uint32_t old_magic = st->magic_samples[i];

   860       st->magic_samples[i] = (old_length - st->filt_len) / 2;

   861       /* We must copy some of the memory that's no longer used */

   862       /* Copy data going backward */

   863       for (j = 0; j < st->filt_len - 1 + st->magic_samples[i] + old_magic; j++)

   864         st->mem[i * st->mem_alloc_size + j] =

   865             st->mem[i * st->mem_alloc_size + j + st->magic_samples[i]];

   866       st->magic_samples[i] += old_magic;

   867     }

   868   }

   869 

   870 }

   871 

   872 EXPORT SpeexResamplerState *

   873 speex_resampler_init (spx_uint32_t nb_channels, spx_uint32_t in_rate,

   874     spx_uint32_t out_rate, int quality, int *err)

   875 {

   876   return speex_resampler_init_frac (nb_channels, in_rate, out_rate, in_rate,

   877       out_rate, quality, err);

   878 }

   879 

   880 EXPORT SpeexResamplerState *

   881 speex_resampler_init_frac (spx_uint32_t nb_channels, spx_uint32_t ratio_num,

   882     spx_uint32_t ratio_den, spx_uint32_t in_rate, spx_uint32_t out_rate,

   883     int quality, int *err)

   884 {

   885   spx_uint32_t i;

   886   SpeexResamplerState *st;

   887   if (quality > 10 || quality < 0) {

   888     if (err)

   889       *err = RESAMPLER_ERR_INVALID_ARG;

   890     return NULL;

   891   }

   892   st = (SpeexResamplerState *) speex_alloc (sizeof (SpeexResamplerState));

   893   st->initialised = 0;

   894   st->started = 0;

   895   st->in_rate = 0;

   896   st->out_rate = 0;

   897   st->num_rate = 0;

   898   st->den_rate = 0;

   899   st->quality = -1;

   900   st->sinc_table_length = 0;

   901   st->mem_alloc_size = 0;

   902   st->filt_len = 0;

   903   st->mem = 0;

   904   st->resampler_ptr = 0;

   905 

   906   st->cutoff = 1.f;

   907   st->nb_channels = nb_channels;

   908   st->in_stride = 1;

   909   st->out_stride = 1;

   910 

   911 #ifdef FIXED_POINT

   912   st->buffer_size = 160;

   913 #else

   914   st->buffer_size = 160;

   915 #endif

   916 

   917   /* Per channel data */

   918   st->last_sample = (spx_int32_t *) speex_alloc (nb_channels * sizeof (int));

   919   st->magic_samples = (spx_uint32_t *) speex_alloc (nb_channels * sizeof (int));

   920   st->samp_frac_num = (spx_uint32_t *) speex_alloc (nb_channels * sizeof (int));

   921   for (i = 0; i < nb_channels; i++) {

   922     st->last_sample[i] = 0;

   923     st->magic_samples[i] = 0;

   924     st->samp_frac_num[i] = 0;

   925   }

   926 

   927   speex_resampler_set_quality (st, quality);

   928   speex_resampler_set_rate_frac (st, ratio_num, ratio_den, in_rate, out_rate);

   929 

   930 

   931   update_filter (st);

   932 

   933   st->initialised = 1;

   934   if (err)

   935     *err = RESAMPLER_ERR_SUCCESS;

   936 

   937   return st;

   938 }

   939 

   940 EXPORT void

   941 speex_resampler_destroy (SpeexResamplerState * st)

   942 {

   943   speex_free (st->mem);

   944   speex_free (st->sinc_table);

   945   speex_free (st->last_sample);

   946   speex_free (st->magic_samples);

   947   speex_free (st->samp_frac_num);

   948   speex_free (st);

   949 }

   950 

   951 static int

   952 speex_resampler_process_native (SpeexResamplerState * st,

   953     spx_uint32_t channel_index, spx_uint32_t * in_len, spx_word16_t * out,

   954     spx_uint32_t * out_len)

   955 {

   956   int j = 0;

   957   const int N = st->filt_len;

   958   int out_sample = 0;

   959   spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;

   960   spx_uint32_t ilen;

   961 

   962   st->started = 1;

   963 

   964   /* Call the right resampler through the function ptr */

   965   out_sample = st->resampler_ptr (st, channel_index, mem, in_len, out, out_len);

   966 

   967   if (st->last_sample[channel_index] < (spx_int32_t) * in_len)

   968     *in_len = st->last_sample[channel_index];

   969   *out_len = out_sample;

   970   st->last_sample[channel_index] -= *in_len;

   971 

   972   ilen = *in_len;

   973 

   974   for (j = 0; j < N - 1; ++j)

   975     mem[j] = mem[j + ilen];

   976 

   977   return RESAMPLER_ERR_SUCCESS;

   978 }

   979 

   980 static int

   981 speex_resampler_magic (SpeexResamplerState * st, spx_uint32_t channel_index,

   982     spx_word16_t ** out, spx_uint32_t out_len)

   983 {

   984   spx_uint32_t tmp_in_len = st->magic_samples[channel_index];

   985   spx_word16_t *mem = st->mem + channel_index * st->mem_alloc_size;

   986   const int N = st->filt_len;

   987 

   988   speex_resampler_process_native (st, channel_index, &tmp_in_len, *out,

   989       &out_len);

   990 

   991   st->magic_samples[channel_index] -= tmp_in_len;

   992 

   993   /* If we couldn't process all "magic" input samples, save the rest for next time */

   994   if (st->magic_samples[channel_index]) {

   995     spx_uint32_t i;

   996     for (i = 0; i < st->magic_samples[channel_index]; i++)

   997       mem[N - 1 + i] = mem[N - 1 + i + tmp_in_len];

   998   }

   999   *out += out_len * st->out_stride;

  1000   return out_len;

  1001 }

  1002 

  1003 #ifdef FIXED_POINT

  1004 EXPORT int

  1005 speex_resampler_process_int (SpeexResamplerState * st,

  1006     spx_uint32_t channel_index, const spx_int16_t * in, spx_uint32_t * in_len,

  1007     spx_int16_t * out, spx_uint32_t * out_len)

  1008 #else

  1009 #ifdef DOUBLE_PRECISION

  1010 EXPORT int

  1011 speex_resampler_process_float (SpeexResamplerState * st,

  1012     spx_uint32_t channel_index, const double *in, spx_uint32_t * in_len,

  1013     double *out, spx_uint32_t * out_len)

  1014 #else

  1015 EXPORT int

  1016 speex_resampler_process_float (SpeexResamplerState * st,

  1017     spx_uint32_t channel_index, const float *in, spx_uint32_t * in_len,

  1018     float *out, spx_uint32_t * out_len)

  1019 #endif

  1020 #endif

  1021 {

  1022   int j;

  1023   spx_uint32_t ilen = *in_len;

  1024   spx_uint32_t olen = *out_len;

  1025   spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;

  1026   const int filt_offs = st->filt_len - 1;

  1027   const spx_uint32_t xlen = st->mem_alloc_size - filt_offs;

  1028   const int istride = st->in_stride;

  1029 

  1030   if (st->magic_samples[channel_index])

  1031     olen -= speex_resampler_magic (st, channel_index, &out, olen);

  1032   if (!st->magic_samples[channel_index]) {

  1033     while (ilen && olen) {

  1034       spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;

  1035       spx_uint32_t ochunk = olen;

  1036 

  1037       if (in) {

  1038         for (j = 0; j < ichunk; ++j)

  1039           x[j + filt_offs] = in[j * istride];

  1040       } else {

  1041         for (j = 0; j < ichunk; ++j)

  1042           x[j + filt_offs] = 0;

  1043       }

  1044       speex_resampler_process_native (st, channel_index, &ichunk, out, &ochunk);

  1045       ilen -= ichunk;

  1046       olen -= ochunk;

  1047       out += ochunk * st->out_stride;

  1048       if (in)

  1049         in += ichunk * istride;

  1050     }

  1051   }

  1052   *in_len -= ilen;

  1053   *out_len -= olen;

  1054   return RESAMPLER_ERR_SUCCESS;

  1055 }

  1056 

  1057 #ifdef FIXED_POINT

  1058 EXPORT int

  1059 speex_resampler_process_float (SpeexResamplerState * st,

  1060     spx_uint32_t channel_index, const float *in, spx_uint32_t * in_len,

  1061     float *out, spx_uint32_t * out_len)

  1062 #else

  1063 EXPORT int

  1064 speex_resampler_process_int (SpeexResamplerState * st,

  1065     spx_uint32_t channel_index, const spx_int16_t * in, spx_uint32_t * in_len,

  1066     spx_int16_t * out, spx_uint32_t * out_len)

  1067 #endif

  1068 {

  1069   int j;

  1070   const int istride_save = st->in_stride;

  1071   const int ostride_save = st->out_stride;

  1072   spx_uint32_t ilen = *in_len;

  1073   spx_uint32_t olen = *out_len;

  1074   spx_word16_t *x = st->mem + channel_index * st->mem_alloc_size;

  1075   const spx_uint32_t xlen = st->mem_alloc_size - (st->filt_len - 1);

  1076 #ifdef VAR_ARRAYS

  1077   const unsigned int ylen =

  1078       (olen < FIXED_STACK_ALLOC) ? olen : FIXED_STACK_ALLOC;

  1079   VARDECL (spx_word16_t * ystack);

  1080   ALLOC (ystack, ylen, spx_word16_t);

  1081 #else

  1082   const unsigned int ylen = FIXED_STACK_ALLOC;

  1083   spx_word16_t ystack[FIXED_STACK_ALLOC];

  1084 #endif

  1085 

  1086   st->out_stride = 1;

  1087 

  1088   while (ilen && olen) {

  1089     spx_word16_t *y = ystack;

  1090     spx_uint32_t ichunk = (ilen > xlen) ? xlen : ilen;

  1091     spx_uint32_t ochunk = (olen > ylen) ? ylen : olen;

  1092     spx_uint32_t omagic = 0;

  1093 

  1094     if (st->magic_samples[channel_index]) {

  1095       omagic = speex_resampler_magic (st, channel_index, &y, ochunk);

  1096       ochunk -= omagic;

  1097       olen -= omagic;

  1098     }

  1099     if (!st->magic_samples[channel_index]) {

  1100       if (in) {

  1101         for (j = 0; j < ichunk; ++j)

  1102 #ifdef FIXED_POINT

  1103           x[j + st->filt_len - 1] = WORD2INT (in[j * istride_save]);

  1104 #else

  1105           x[j + st->filt_len - 1] = in[j * istride_save];

  1106 #endif

  1107       } else {

  1108         for (j = 0; j < ichunk; ++j)

  1109           x[j + st->filt_len - 1] = 0;

  1110       }

  1111 

  1112       speex_resampler_process_native (st, channel_index, &ichunk, y, &ochunk);

  1113     } else {

  1114       ichunk = 0;

  1115       ochunk = 0;

  1116     }

  1117 

  1118     for (j = 0; j < ochunk + omagic; ++j)

  1119 #ifdef FIXED_POINT

  1120       out[j * ostride_save] = ystack[j];

  1121 #else

  1122       out[j * ostride_save] = WORD2INT (ystack[j]);

  1123 #endif

  1124 

  1125     ilen -= ichunk;

  1126     olen -= ochunk;

  1127     out += (ochunk + omagic) * ostride_save;

  1128     if (in)

  1129       in += ichunk * istride_save;

  1130   }

  1131   st->out_stride = ostride_save;

  1132   *in_len -= ilen;

  1133   *out_len -= olen;

  1134