/* GStreamer

 * Copyright (C) <2015> Wim Taymans <wim.taymans@gmail.com>

 *

 * This library is free software; you can redistribute it and/or

 * modify it under the terms of the GNU Library General Public

 * License as published by the Free Software Foundation; either

 * version 2 of the License, or (at your option) any later version.

 *

 * This library is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * Library General Public License for more details.

 *

 * You should have received a copy of the GNU Library General Public

 * License along with this library; if not, write to the

 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,

 * Boston, MA 02110-1301, USA.

 */

#ifdef HAVE_CONFIG_H

#  include "config.h"

#endif

#include <string.h>

#include <stdio.h>

#include <math.h>

#ifdef HAVE_ORC

#include <orc/orc.h>

#endif

#include "audio-resampler.h"

#include "audio-resampler-private.h"

#include "audio-resampler-macros.h"

#define MEM_ALIGN(m,a) ((gint8 *)((guintptr)((gint8 *)(m) + ((a)-1)) & ~((a)-1)))

#define ALIGN 16

#define TAPS_OVERREAD 16

GST_DEBUG_CATEGORY_STATIC (audio_resampler_debug);

#define GST_CAT_DEFAULT audio_resampler_debug

/**

 * SECTION:gstaudioresampler

 * @title: GstAudioResampler

 * @short_description: Utility structure for resampler information

 *

 * #GstAudioResampler is a structure which holds the information

 * required to perform various kinds of resampling filtering.

 *

 */

static const gint oversample_qualities[] = {

  4, 4, 4, 8, 8, 16, 16, 16, 16, 32, 32

};

typedef struct

{

  gdouble cutoff;

  gdouble downsample_cutoff_factor;

  gdouble stopband_attenuation;

  gdouble transition_bandwidth;

} KaiserQualityMap;

static const KaiserQualityMap kaiser_qualities[] = {

  {0.860, 0.96511, 60, 0.7},    /* 8 taps */

  {0.880, 0.96591, 65, 0.29},   /* 16 taps */

  {0.910, 0.96923, 70, 0.145},  /* 32 taps */

  {0.920, 0.97600, 80, 0.105},  /* 48 taps */

  {0.940, 0.97979, 85, 0.087},  /* 64 taps default quality */

  {0.940, 0.98085, 95, 0.077},  /* 80 taps */

  {0.945, 0.99471, 100, 0.068}, /* 96 taps */

  {0.950, 1.0, 105, 0.055},     /* 128 taps */

  {0.960, 1.0, 110, 0.045},     /* 160 taps */

  {0.968, 1.0, 115, 0.039},     /* 192 taps */

  {0.975, 1.0, 120, 0.0305}     /* 256 taps */

};

typedef struct

{

  gint n_taps;

  gdouble cutoff;

} BlackmanQualityMap;

static const BlackmanQualityMap blackman_qualities[] = {

  {8, 0.5,},

  {16, 0.6,},

  {24, 0.72,},

  {32, 0.8,},

  {48, 0.85,},                  /* default */

  {64, 0.90,},

  {80, 0.92,},

  {96, 0.933,},

  {128, 0.950,},

  {148, 0.955,},

  {160, 0.960,}

};

#define DEFAULT_RESAMPLER_METHOD GST_AUDIO_RESAMPLER_METHOD_KAISER

#define DEFAULT_QUALITY GST_AUDIO_RESAMPLER_QUALITY_DEFAULT

#define DEFAULT_OPT_CUBIC_B 1.0

#define DEFAULT_OPT_CUBIC_C 0.0

#define DEFAULT_OPT_FILTER_MODE GST_AUDIO_RESAMPLER_FILTER_MODE_AUTO

#define DEFAULT_OPT_FILTER_MODE_THRESHOLD 1048576

#define DEFAULT_OPT_FILTER_INTERPOLATION GST_AUDIO_RESAMPLER_FILTER_INTERPOLATION_CUBIC

#define DEFAULT_OPT_FILTER_OVERSAMPLE 8

#define DEFAULT_OPT_MAX_PHASE_ERROR 0.1

static gdouble

get_opt_double (GstStructure * options, const gchar * name, gdouble def)

{

  gdouble res;

  if (!options || !gst_structure_get_double (options, name, &res))

    res = def;

  return res;

}

static gint

get_opt_int (GstStructure * options, const gchar * name, gint def)

{

  gint res;

  if (!options || !gst_structure_get_int (options, name, &res))

    res = def;

  return res;

}

static gint

get_opt_enum (GstStructure * options, const gchar * name, GType type, gint def)

{

  gint res;

  if (!options || !gst_structure_get_enum (options, name, type, &res))

    res = def;

  return res;

}

#define GET_OPT_CUTOFF(options,def) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_CUTOFF,def)

#define GET_OPT_DOWN_CUTOFF_FACTOR(options,def) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_DOWN_CUTOFF_FACTOR, def)

#define GET_OPT_STOP_ATTENUATION(options,def) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_STOP_ATTENUATION, def)

#define GET_OPT_TRANSITION_BANDWIDTH(options,def) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_TRANSITION_BANDWIDTH, def)

#define GET_OPT_CUBIC_B(options) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_CUBIC_B, DEFAULT_OPT_CUBIC_B)

#define GET_OPT_CUBIC_C(options) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_CUBIC_C, DEFAULT_OPT_CUBIC_C)

#define GET_OPT_N_TAPS(options,def) get_opt_int(options, \

    GST_AUDIO_RESAMPLER_OPT_N_TAPS, def)

#define GET_OPT_FILTER_MODE(options) get_opt_enum(options, \

    GST_AUDIO_RESAMPLER_OPT_FILTER_MODE, GST_TYPE_AUDIO_RESAMPLER_FILTER_MODE, \

    DEFAULT_OPT_FILTER_MODE)

#define GET_OPT_FILTER_MODE_THRESHOLD(options) get_opt_int(options, \

    GST_AUDIO_RESAMPLER_OPT_FILTER_MODE_THRESHOLD, DEFAULT_OPT_FILTER_MODE_THRESHOLD)

#define GET_OPT_FILTER_INTERPOLATION(options) get_opt_enum(options, \

    GST_AUDIO_RESAMPLER_OPT_FILTER_INTERPOLATION, GST_TYPE_AUDIO_RESAMPLER_FILTER_INTERPOLATION, \

    DEFAULT_OPT_FILTER_INTERPOLATION)

#define GET_OPT_FILTER_OVERSAMPLE(options) get_opt_int(options, \

    GST_AUDIO_RESAMPLER_OPT_FILTER_OVERSAMPLE, DEFAULT_OPT_FILTER_OVERSAMPLE)

#define GET_OPT_MAX_PHASE_ERROR(options) get_opt_double(options, \

    GST_AUDIO_RESAMPLER_OPT_MAX_PHASE_ERROR, DEFAULT_OPT_MAX_PHASE_ERROR)

#include "dbesi0.c"

#define bessel dbesi0

static inline gdouble

get_linear_tap (gdouble x, gint n_taps)

{

  gdouble res = GST_ROUND_UP_2 (n_taps) / 2 - fabs (x);

  return res;

}

static inline gdouble

get_cubic_tap (gdouble x, gint n_taps, gdouble b, gdouble c)

{

  gdouble res, a, a2, a3;

  a = fabs (x * 4.0) / n_taps;

  a2 = a * a;

  a3 = a2 * a;

  if (a <= 1.0)

    res = ((12.0 - 9.0 * b - 6.0 * c) * a3 +

        (-18.0 + 12.0 * b + 6.0 * c) * a2 + (6.0 - 2.0 * b)) / 6.0;

  else if (a <= 2.0)

    res = ((-b - 6.0 * c) * a3 +

        (6.0 * b + 30.0 * c) * a2 +

        (-12.0 * b - 48.0 * c) * a + (8.0 * b + 24.0 * c)) / 6.0;

  else

    res = 0.0;

  return res;

}

static inline gdouble

get_blackman_nuttall_tap (gdouble x, gint n_taps, gdouble Fc)

{

  gdouble s, y, w;

  y = G_PI * x;

  s = (y == 0.0 ? Fc : sin (y * Fc) / y);

  w = 2.0 * y / n_taps + G_PI;

  return s * (0.3635819 - 0.4891775 * cos (w) + 0.1365995 * cos (2 * w) -

      0.0106411 * cos (3 * w));

}

static inline gdouble

get_kaiser_tap (gdouble x, gint n_taps, gdouble Fc, gdouble beta)

{

  gdouble s, y, w;

  y = G_PI * x;

  s = (y == 0.0 ? Fc : sin (y * Fc) / y);

  w = 2.0 * x / n_taps;

  return s * bessel (beta * sqrt (MAX (1 - w * w, 0)));

}

#define MAKE_CONVERT_TAPS_INT_FUNC(type, precision)                     \

static void                                                             \

convert_taps_##type##_c (gdouble *tmp_taps, gpointer taps,              \

    gdouble weight, gint n_taps)                                        \

{                                                                       \

  gint64 one = (1LL << precision) - 1;                                  \

  type *t = taps;                                                       \

  gdouble multiplier = one;                                             \

  gint i, j;                                                            \

  gdouble offset, l_offset, h_offset;                                   \

  gboolean exact = FALSE;                                               \

  /* Round to integer, but with an adjustable bias that we use to */    \

  /* eliminate the DC error. */                                         \

  l_offset = 0.0;                                                       \

  h_offset = 1.0;                                                       \

  offset = 0.5;                                                         \

  for (i = 0; i < 32; i++) {                                            \

    gint64 sum = 0;                                                     \

    for (j = 0; j < n_taps; j++)                                        \

      sum += floor (offset + tmp_taps[j] * multiplier / weight);        \

    if (sum == one) {                                                   \

      exact = TRUE;                                                     \

      break;                                                            \

    }                                                                   \

    if (l_offset == h_offset)                                           \

      break;                                                            \

    if (sum < one) {                                                    \

      if (offset > l_offset)                                            \

        l_offset = offset;                                              \

      offset += (h_offset - l_offset) / 2;                              \

    } else {                                                            \

      if (offset < h_offset)                                            \

        h_offset = offset;                                              \

      offset -= (h_offset - l_offset) / 2;                              \

    }                                                                   \

  }                                                                     \

  for (j = 0; j < n_taps; j++)                                          \

    t[j] = floor (offset + tmp_taps[j] * multiplier / weight);          \

  if (!exact)                                                           \

    GST_WARNING ("can't find exact taps");                              \

}

#define MAKE_CONVERT_TAPS_FLOAT_FUNC(type)                              \

static void                                                             \

convert_taps_##type##_c (gdouble *tmp_taps, gpointer taps,              \

    gdouble weight, gint n_taps)                                        \

{                                                                       \

  gint i;                                                               \

  type *t = taps;                                                       \

  for (i = 0; i < n_taps; i++)                                          \

    t[i] = tmp_taps[i] / weight;                                        \

}

MAKE_CONVERT_TAPS_INT_FUNC (gint16, PRECISION_S16);

MAKE_CONVERT_TAPS_INT_FUNC (gint32, PRECISION_S32);

MAKE_CONVERT_TAPS_FLOAT_FUNC (gfloat);

MAKE_CONVERT_TAPS_FLOAT_FUNC (gdouble);

static ConvertTapsFunc convert_taps_funcs[] = {

  convert_taps_gint16_c,

  convert_taps_gint32_c,

  convert_taps_gfloat_c,

  convert_taps_gdouble_c

};

#define convert_taps_gint16   convert_taps_funcs[0]

#define convert_taps_gint32   convert_taps_funcs[1]

#define convert_taps_gfloat   convert_taps_funcs[2]

#define convert_taps_gdouble  convert_taps_funcs[3]

static void

make_taps (GstAudioResampler * resampler, gdouble * res, gdouble x, gint n_taps)

{

  gdouble weight = 0.0, *tmp_taps = resampler->tmp_taps;

  gint i;

  switch (resampler->method) {

    case GST_AUDIO_RESAMPLER_METHOD_NEAREST:

      break;

    case GST_AUDIO_RESAMPLER_METHOD_LINEAR:

      for (i = 0; i < n_taps; i++)

        weight += tmp_taps[i] = get_linear_tap (x + i, resampler->n_taps);

      break;

    case GST_AUDIO_RESAMPLER_METHOD_CUBIC:

      for (i = 0; i < n_taps; i++)

        weight += tmp_taps[i] = get_cubic_tap (x + i, resampler->n_taps,

            resampler->b, resampler->c);

      break;

    case GST_AUDIO_RESAMPLER_METHOD_BLACKMAN_NUTTALL:

      for (i = 0; i < n_taps; i++)

        weight += tmp_taps[i] =

            get_blackman_nuttall_tap (x + i,

            resampler->n_taps, resampler->cutoff);

      break;

    case GST_AUDIO_RESAMPLER_METHOD_KAISER:

      for (i = 0; i < n_taps; i++)

        weight += tmp_taps[i] =

            get_kaiser_tap (x + i, resampler->n_taps,

            resampler->cutoff, resampler->kaiser_beta);

      break;

  }

  resampler->convert_taps (tmp_taps, res, weight, n_taps);

}

#define MAKE_COEFF_LINEAR_INT_FUNC(type,type2,prec)                     \

static inline void                                                      \

make_coeff_##type##_linear (gint num, gint denom, type *icoeff)         \

{                                                                       \

  type x = ((gint64)num << prec) / denom;                               \

  icoeff[0] = icoeff[2] = x;                                            \

  icoeff[1] = icoeff[3] = (type)(((type2)1 << prec)-1)  - x;            \

}

#define MAKE_COEFF_LINEAR_FLOAT_FUNC(type)                              \

static inline void                                                      \

make_coeff_##type##_linear (gint num, gint denom, type *icoeff)         \

{                                                                       \

  type x = (type)num / denom;                                           \

  icoeff[0] = icoeff[2] = x;                                            \

  icoeff[1] = icoeff[3] = (type)1.0 - x;                                \

}

MAKE_COEFF_LINEAR_INT_FUNC (gint16, gint32, PRECISION_S16);

MAKE_COEFF_LINEAR_INT_FUNC (gint32, gint64, PRECISION_S32);

MAKE_COEFF_LINEAR_FLOAT_FUNC (gfloat);

MAKE_COEFF_LINEAR_FLOAT_FUNC (gdouble);

#define MAKE_COEFF_CUBIC_INT_FUNC(type,type2,prec)                      \

static inline void                                                      \

make_coeff_##type##_cubic (gint num, gint denom, type *icoeff)          \

{                                                                       \

  type2 one = ((type2)1 << prec) - 1;                                   \

  type2 x = ((gint64) num << prec) / denom;                             \

  type2 x2 = (x * x) >> prec;                                           \

  type2 x3 = (x2 * x) >> prec;                                          \

  icoeff[0] = (((x3 - x) << prec) / 6) >> prec;                         \

  icoeff[1] = x + ((x2 - x3) >> 1);                                     \

  icoeff[3] = -(((x << prec) / 3) >> prec) +                            \

            (x2 >> 1) - (((x3 << prec) / 6) >> prec);                   \

  icoeff[2] = one - icoeff[0] - icoeff[1] - icoeff[3];                  \

}

#define MAKE_COEFF_CUBIC_FLOAT_FUNC(type)                               \

static inline void                                                      \

make_coeff_##type##_cubic (gint num, gint denom, type *icoeff)          \

{                                                                       \

  type x = (type) num / denom, x2 = x * x, x3 = x2 * x;                 \

  icoeff[0] = 0.16667f * (x3 - x);                                      \

  icoeff[1] = x + 0.5f * (x2 - x3);                                     \

  icoeff[3] = -0.33333f * x + 0.5f * x2 - 0.16667f * x3;                \

  icoeff[2] = (type)1.0 - icoeff[0] - icoeff[1] - icoeff[3];            \

}

MAKE_COEFF_CUBIC_INT_FUNC (gint16, gint32, PRECISION_S16);

MAKE_COEFF_CUBIC_INT_FUNC (gint32, gint64, PRECISION_S32);

MAKE_COEFF_CUBIC_FLOAT_FUNC (gfloat);

MAKE_COEFF_CUBIC_FLOAT_FUNC (gdouble);

#define INTERPOLATE_INT_LINEAR_FUNC(type,type2,prec,limit)      \

static inline void                                              \

interpolate_##type##_linear_c (gpointer op, const gpointer ap,  \

    gint len, const gpointer icp, gint astride)                 \

{                                                               \

  gint i;                                                       \

  type *o = op, *a = ap, *ic = icp;                             \

  type2 tmp, c0 = ic[0];                                        \

  const type *c[2] = {(type*)((gint8*)a + 0*astride),           \

                      (type*)((gint8*)a + 1*astride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    tmp = ((type2)c[0][i] - (type2)c[1][i]) * c0 +              \

         (((type2)c[1][i]) << (prec));                          \

    o[i] = (tmp + ((type2)1 << ((prec) - 1))) >> (prec);        \

  }                                                             \

}

#define INTERPOLATE_FLOAT_LINEAR_FUNC(type)                     \

static inline void                                              \

interpolate_##type##_linear_c (gpointer op, const gpointer ap,  \

    gint len, const gpointer icp, gint astride)                 \

{                                                               \

  gint i;                                                       \

  type *o = op, *a = ap, *ic = icp;                             \

  type c0 = ic[0];                                              \

  const type *c[2] = {(type*)((gint8*)a + 0*astride),           \

                      (type*)((gint8*)a + 1*astride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    o[i] = (c[0][i] - c[1][i]) * c0 + c[1][i];                  \

  }                                                             \

}

INTERPOLATE_INT_LINEAR_FUNC (gint16, gint32, PRECISION_S16, (gint32) 1 << 15);

INTERPOLATE_INT_LINEAR_FUNC (gint32, gint64, PRECISION_S32, (gint64) 1 << 31);

INTERPOLATE_FLOAT_LINEAR_FUNC (gfloat);

INTERPOLATE_FLOAT_LINEAR_FUNC (gdouble);

#define INTERPOLATE_INT_CUBIC_FUNC(type,type2,prec,limit)       \

static inline void                                              \

interpolate_##type##_cubic_c (gpointer op, const gpointer ap,   \

    gint len, const gpointer icp, gint astride)                 \

{                                                               \

  gint i;                                                       \

  type *o = op, *a = ap, *ic = icp;                             \

  type2 tmp, c0 = ic[0], c1 = ic[1], c2 = ic[2], c3 = ic[3];    \

  const type *c[4] = {(type*)((gint8*)a + 0*astride),           \

                      (type*)((gint8*)a + 1*astride),           \

                      (type*)((gint8*)a + 2*astride),           \

                      (type*)((gint8*)a + 3*astride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    tmp = (type2)c[0][i] * c0 + (type2)c[1][i] * c1 +           \

          (type2)c[2][i] * c2 + (type2)c[3][i] * c3;            \

    tmp = (tmp + ((type2)1 << ((prec) - 1))) >> (prec);         \

    o[i] = CLAMP (tmp, -(limit), (limit) - 1);                  \

  }                                                             \

}

#define INTERPOLATE_FLOAT_CUBIC_FUNC(type)                      \

static inline void                                              \

interpolate_##type##_cubic_c (gpointer op, const gpointer ap,   \

    gint len, const gpointer icp, gint astride)                 \

{                                                               \

  gint i;                                                       \

  type *o = op, *a = ap, *ic = icp;                             \

  type c0 = ic[0], c1 = ic[1], c2 = ic[2], c3 = ic[3];          \

  const type *c[4] = {(type*)((gint8*)a + 0*astride),           \

                      (type*)((gint8*)a + 1*astride),           \

                      (type*)((gint8*)a + 2*astride),           \

                      (type*)((gint8*)a + 3*astride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    o[i] = c[0][i] * c0 + c[1][i] * c1 +                        \

           c[2][i] * c2 + c[3][i] * c3;                         \

  }                                                             \

}

INTERPOLATE_INT_CUBIC_FUNC (gint16, gint32, PRECISION_S16, (gint32) 1 << 15);

INTERPOLATE_INT_CUBIC_FUNC (gint32, gint64, PRECISION_S32, (gint64) 1 << 31);

INTERPOLATE_FLOAT_CUBIC_FUNC (gfloat);

INTERPOLATE_FLOAT_CUBIC_FUNC (gdouble);

static InterpolateFunc interpolate_funcs[] = {

  interpolate_gint16_linear_c,

  interpolate_gint32_linear_c,

  interpolate_gfloat_linear_c,

  interpolate_gdouble_linear_c,

  interpolate_gint16_cubic_c,

  interpolate_gint32_cubic_c,

  interpolate_gfloat_cubic_c,

  interpolate_gdouble_cubic_c,

};

#define interpolate_gint16_linear  interpolate_funcs[0]

#define interpolate_gint32_linear  interpolate_funcs[1]

#define interpolate_gfloat_linear  interpolate_funcs[2]

#define interpolate_gdouble_linear interpolate_funcs[3]

#define interpolate_gint16_cubic   interpolate_funcs[4]

#define interpolate_gint32_cubic   interpolate_funcs[5]

#define interpolate_gfloat_cubic   interpolate_funcs[6]

#define interpolate_gdouble_cubic  interpolate_funcs[7]

#define GET_TAPS_NEAREST_FUNC(type)                                             \

static inline gpointer                                                          \

get_taps_##type##_nearest (GstAudioResampler * resampler,                       \

    gint *samp_index, gint *samp_phase, type icoeff[4])                         \

{                                                                               \

  gint out_rate = resampler->out_rate;                                          \

  *samp_index += resampler->samp_inc;                                           \

  *samp_phase += resampler->samp_frac;                                          \

  if (*samp_phase >= out_rate) {                                                \

    *samp_phase -= out_rate;                                                    \

    *samp_index += 1;                                                           \

  }                                                                             \

  return NULL;                                                                  \

}

GET_TAPS_NEAREST_FUNC (gint16);

GET_TAPS_NEAREST_FUNC (gint32);

GET_TAPS_NEAREST_FUNC (gfloat);

GET_TAPS_NEAREST_FUNC (gdouble);

#define get_taps_gint16_nearest get_taps_gint16_nearest

#define get_taps_gint32_nearest get_taps_gint32_nearest

#define get_taps_gfloat_nearest get_taps_gfloat_nearest

#define get_taps_gdouble_nearest get_taps_gdouble_nearest

#define GET_TAPS_FULL_FUNC(type)                                                \

DECL_GET_TAPS_FULL_FUNC(type)                                                   \

{                                                                               \

  gpointer res;                                                                 \

  gint out_rate = resampler->out_rate;                                          \

  gint n_phases = resampler->n_phases;                                          \

  gint phase = (n_phases == out_rate ? *samp_phase :                            \

      ((gint64)*samp_phase * n_phases) / out_rate);                             \

                                                                                \

  res = resampler->cached_phases[phase];                                        \

  if (G_UNLIKELY (res == NULL)) {                                               \

    res = (gint8 *) resampler->cached_taps +                                    \

                        phase * resampler->cached_taps_stride;                  \

    switch (resampler->filter_interpolation) {                                  \

      case GST_AUDIO_RESAMPLER_FILTER_INTERPOLATION_NONE:                       \

      {                                                                         \

        gdouble x;                                                              \

        gint n_taps = resampler->n_taps;                                        \

                                                                                \

        x = 1.0 - n_taps / 2 - (gdouble) phase / n_phases;                      \

        make_taps (resampler, res, x, n_taps);                                  \

        break;                                                                  \

      }                                                                         \

      default:                                                                  \

      {                                                                         \

        gint offset, pos, frac;                                                 \

        gint oversample = resampler->oversample;                                \

        gint taps_stride = resampler->taps_stride;                              \

        gint n_taps = resampler->n_taps;                                        \

        type ic[4], *taps;                                                      \

                                                                                \

        pos = phase * oversample;                                               \

        offset = (oversample - 1) - pos / n_phases;                             \

        frac = pos % n_phases;                                                  \

                                                                                \

        taps = (type *) ((gint8 *) resampler->taps + offset * taps_stride);     \

                                                                                \

        switch (resampler->filter_interpolation) {                              \

          default:                                                              \

          case GST_AUDIO_RESAMPLER_FILTER_INTERPOLATION_LINEAR:                 \

            make_coeff_##type##_linear (frac, n_phases, ic);                    \

            break;                                                              \

          case GST_AUDIO_RESAMPLER_FILTER_INTERPOLATION_CUBIC:                  \

            make_coeff_##type##_cubic (frac, n_phases, ic);                     \

            break;                                                              \

        }                                                                       \

        resampler->interpolate (res, taps, n_taps, ic, taps_stride);            \

      }                                                                         \

    }                                                                           \

    resampler->cached_phases[phase] = res;                                      \

  }                                                                             \

  *samp_index += resampler->samp_inc;                                           \

  *samp_phase += resampler->samp_frac;                                          \

  if (*samp_phase >= out_rate) {                                                \

    *samp_phase -= out_rate;                                                    \

    *samp_index += 1;                                                           \

  }                                                                             \

  return res;                                                                   \

}

GET_TAPS_FULL_FUNC (gint16);

GET_TAPS_FULL_FUNC (gint32);

GET_TAPS_FULL_FUNC (gfloat);

GET_TAPS_FULL_FUNC (gdouble);

#define GET_TAPS_INTERPOLATE_FUNC(type,inter)                   \

DECL_GET_TAPS_INTERPOLATE_FUNC (type, inter)                    \

{                                                               \

  gpointer res;                                                 \

  gint out_rate = resampler->out_rate;                          \

  gint offset, frac, pos;                                       \

  gint oversample = resampler->oversample;                      \

  gint taps_stride = resampler->taps_stride;                    \

                                                                \

  pos = *samp_phase * oversample;                               \

  offset = (oversample - 1) - pos / out_rate;                   \

  frac = pos % out_rate;                                        \

                                                                \

  res = (gint8 *) resampler->taps + offset * taps_stride;       \

  make_coeff_##type##_##inter (frac, out_rate, icoeff);         \

                                                                \

  *samp_index += resampler->samp_inc;                           \

  *samp_phase += resampler->samp_frac;                          \

  if (*samp_phase >= out_rate) {                                \

    *samp_phase -= out_rate;                                    \

    *samp_index += 1;                                           \

  }                                                             \

  return res;                                                   \

}

GET_TAPS_INTERPOLATE_FUNC (gint16, linear);

GET_TAPS_INTERPOLATE_FUNC (gint32, linear);

GET_TAPS_INTERPOLATE_FUNC (gfloat, linear);

GET_TAPS_INTERPOLATE_FUNC (gdouble, linear);

GET_TAPS_INTERPOLATE_FUNC (gint16, cubic);

GET_TAPS_INTERPOLATE_FUNC (gint32, cubic);

GET_TAPS_INTERPOLATE_FUNC (gfloat, cubic);

GET_TAPS_INTERPOLATE_FUNC (gdouble, cubic);

#define INNER_PRODUCT_NEAREST_FUNC(type)                        \

static inline void                                              \

inner_product_##type##_nearest_1_c (type * o, const type * a,   \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  *o = *a;                                                      \

}

INNER_PRODUCT_NEAREST_FUNC (gint16);

INNER_PRODUCT_NEAREST_FUNC (gint32);

INNER_PRODUCT_NEAREST_FUNC (gfloat);

INNER_PRODUCT_NEAREST_FUNC (gdouble);

#define INNER_PRODUCT_INT_FULL_FUNC(type,type2,prec,limit)      \

static inline void                                              \

inner_product_##type##_full_1_c (type * o, const type * a,      \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type2 res[4] = { 0, 0, 0, 0 };                                \

                                                                \

  for (i = 0; i < len; i += 4) {                                \

    res[0] += (type2) a[i + 0] * (type2) b[i + 0];              \

    res[1] += (type2) a[i + 1] * (type2) b[i + 1];              \

    res[2] += (type2) a[i + 2] * (type2) b[i + 2];              \

    res[3] += (type2) a[i + 3] * (type2) b[i + 3];              \

  }                                                             \

  res[0] = res[0] + res[1] + res[2] + res[3];                   \

  res[0] = (res[0] + ((type2)1 << ((prec) - 1))) >> (prec);     \

  *o = CLAMP (res[0], -(limit), (limit) - 1);                   \

}

INNER_PRODUCT_INT_FULL_FUNC (gint16, gint32, PRECISION_S16, (gint32) 1 << 15);

INNER_PRODUCT_INT_FULL_FUNC (gint32, gint64, PRECISION_S32, (gint64) 1 << 31);

#define INNER_PRODUCT_INT_LINEAR_FUNC(type,type2,prec,limit)    \

static inline void                                              \

inner_product_##type##_linear_1_c (type * o, const type * a,    \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type2 res[4] = { 0, 0, 0, 0 }, c0 = ic[0];                    \

  const type *c[2] = {(type*)((gint8*)b + 0*bstride),           \

                      (type*)((gint8*)b + 1*bstride)};          \

                                                                \

  for (i = 0; i < len; i += 2) {                                \

    res[0] += (type2) a[i + 0] * (type2) c[0][i + 0];           \

    res[1] += (type2) a[i + 0] * (type2) c[1][i + 0];           \

    res[2] += (type2) a[i + 1] * (type2) c[0][i + 1];           \

    res[3] += (type2) a[i + 1] * (type2) c[1][i + 1];           \

  }                                                             \

  res[0] = (res[0] + res[2]) >> (prec);                         \

  res[1] = (res[1] + res[3]) >> (prec);                         \

  res[0] = ((type2)(type)res[0] - (type2)(type)res[1]) * c0 +   \

           ((type2)(type)res[1] << (prec));                     \

  res[0] = (res[0] + ((type2)1 << ((prec) - 1))) >> (prec);     \

  *o = CLAMP (res[0], -(limit), (limit) - 1);                   \

}

INNER_PRODUCT_INT_LINEAR_FUNC (gint16, gint32, PRECISION_S16, (gint32) 1 << 15);

INNER_PRODUCT_INT_LINEAR_FUNC (gint32, gint64, PRECISION_S32, (gint64) 1 << 31);

#define INNER_PRODUCT_INT_CUBIC_FUNC(type,type2,prec,limit)     \

static inline void                                              \

inner_product_##type##_cubic_1_c (type * o, const type * a,     \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type2 res[4] = { 0, 0, 0, 0 };                                \

  const type *c[4] = {(type*)((gint8*)b + 0*bstride),           \

                      (type*)((gint8*)b + 1*bstride),           \

                      (type*)((gint8*)b + 2*bstride),           \

                      (type*)((gint8*)b + 3*bstride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    res[0] += (type2) a[i] * (type2) c[0][i];                   \

    res[1] += (type2) a[i] * (type2) c[1][i];                   \

    res[2] += (type2) a[i] * (type2) c[2][i];                   \

    res[3] += (type2) a[i] * (type2) c[3][i];                   \

  }                                                             \

  res[0] = (type2)(type)(res[0] >> (prec)) * (type2) ic[0] +    \

           (type2)(type)(res[1] >> (prec)) * (type2) ic[1] +    \

           (type2)(type)(res[2] >> (prec)) * (type2) ic[2] +    \

           (type2)(type)(res[3] >> (prec)) * (type2) ic[3];     \

  res[0] = (res[0] + ((type2)1 << ((prec) - 1))) >> (prec);     \

  *o = CLAMP (res[0], -(limit), (limit) - 1);                   \

}

INNER_PRODUCT_INT_CUBIC_FUNC (gint16, gint32, PRECISION_S16, (gint32) 1 << 15);

INNER_PRODUCT_INT_CUBIC_FUNC (gint32, gint64, PRECISION_S32, (gint64) 1 << 31);

#define INNER_PRODUCT_FLOAT_FULL_FUNC(type)                     \

static inline void                                              \

inner_product_##type##_full_1_c (type * o, const type * a,      \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type res[4] = { 0.0, 0.0, 0.0, 0.0 };                         \

                                                                \

  for (i = 0; i < len; i += 4) {                                \

    res[0] += a[i + 0] * b[i + 0];                              \

    res[1] += a[i + 1] * b[i + 1];                              \

    res[2] += a[i + 2] * b[i + 2];                              \

    res[3] += a[i + 3] * b[i + 3];                              \

  }                                                             \

  *o = res[0] + res[1] + res[2] + res[3];                       \

}

INNER_PRODUCT_FLOAT_FULL_FUNC (gfloat);

INNER_PRODUCT_FLOAT_FULL_FUNC (gdouble);

#define INNER_PRODUCT_FLOAT_LINEAR_FUNC(type)                   \

static inline void                                              \

inner_product_##type##_linear_1_c (type * o, const type * a,    \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type res[4] = { 0.0, 0.0, 0.0, 0.0 };                         \

  const type *c[2] = {(type*)((gint8*)b + 0*bstride),           \

                      (type*)((gint8*)b + 1*bstride)};          \

                                                                \

  for (i = 0; i < len; i += 2) {                                \

    res[0] += a[i + 0] * c[0][i + 0];                           \

    res[1] += a[i + 0] * c[1][i + 0];                           \

    res[2] += a[i + 1] * c[0][i + 1];                           \

    res[3] += a[i + 1] * c[1][i + 1];                           \

  }                                                             \

  res[0] += res[2];                                             \

  res[1] += res[3];                                             \

  *o = (res[0] - res[1]) * ic[0] + res[1];                      \

}

INNER_PRODUCT_FLOAT_LINEAR_FUNC (gfloat);

INNER_PRODUCT_FLOAT_LINEAR_FUNC (gdouble);

#define INNER_PRODUCT_FLOAT_CUBIC_FUNC(type)                    \

static inline void                                              \

inner_product_##type##_cubic_1_c (type * o, const type * a,     \

    const type * b, gint len, const type *ic, gint bstride)     \

{                                                               \

  gint i;                                                       \

  type res[4] = { 0.0, 0.0, 0.0, 0.0 };                         \

  const type *c[4] = {(type*)((gint8*)b + 0*bstride),           \

                      (type*)((gint8*)b + 1*bstride),           \

                      (type*)((gint8*)b + 2*bstride),           \

                      (type*)((gint8*)b + 3*bstride)};          \

                                                                \

  for (i = 0; i < len; i++) {                                   \

    res[0] += a[i] * c[0][i];                                   \

    res[1] += a[i] * c[1][i];                                   \

    res[2] += a[i] * c[2][i];                                   \

    res[3] += a[i] * c[3][i];                                   \

  }                                                             \

  *o = res[0] * ic[0] + res[1] * ic[1] +                        \

       res[2] * ic[2] + res[3] * ic[3];                         \

}

INNER_PRODUCT_FLOAT_CUBIC_FUNC (gfloat);

INNER_PRODUCT_FLOAT_CUBIC_FUNC (gdouble);

MAKE_RESAMPLE_FUNC_STATIC (gint16, nearest, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint32, nearest, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gfloat, nearest, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gdouble, nearest, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint16, full, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint32, full, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gfloat, full, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gdouble, full, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint16, linear, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint32, linear, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gfloat, linear, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gdouble, linear, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint16, cubic, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gint32, cubic, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gfloat, cubic, 1, c);

MAKE_RESAMPLE_FUNC_STATIC (gdouble, cubic, 1, c);

static ResampleFunc resample_funcs[] = {

  resample_gint16_nearest_1_c,

  resample_gint32_nearest_1_c,

  resample_gfloat_nearest_1_c,

  resample_gdouble_nearest_1_c,

  resample_gint16_full_1_c,

  resample_gint32_full_1_c,

  resample_gfloat_full_1_c,

  resample_gdouble_full_1_c,

  resample_gint16_linear_1_c,

  resample_gint32_linear_1_c,

  resample_gfloat_linear_1_c,

  resample_gdouble_linear_1_c,

  resample_gint16_cubic_1_c,

  resample_gint32_cubic_1_c,

  resample_gfloat_cubic_1_c,

  resample_gdouble_cubic_1_c,

};

#define resample_gint16_nearest_1 resample_funcs[0]

#define resample_gint32_nearest_1 resample_funcs[1]

#define resample_gfloat_nearest_1 resample_funcs[2]

#define resample_gdouble_nearest_1 resample_funcs[3]

#define resample_gint16_full_1 resample_funcs[4]

#define resample_gint32_full_1 resample_funcs[5]

#define resample_gfloat_full_1 resample_funcs[6]

#define resample_gdouble_full_1 resample_funcs[7]

#define resample_gint16_linear_1 resample_funcs[8]

#define resample_gint32_linear_1 resample_funcs[9]

#define resample_gfloat_linear_1 resample_funcs[10]

#define resample_gdouble_linear_1 resample_funcs[11]

#define resample_gint16_cubic_1 resample_funcs[12]

#define resample_gint32_cubic_1 resample_funcs[13]

#define resample_gfloat_cubic_1 resample_funcs[14]

#define resample_gdouble_cubic_1 resample_funcs[15]

#if defined HAVE_ORC && !defined DISABLE_ORC

# if defined (HAVE_ARM_NEON)

#  define CHECK_NEON

#  include "audio-resampler-neon.h"

# endif

# if defined (__i386__) || defined (__x86_64__)

#  define CHECK_X86

#  include "audio-resampler-x86.h"

# endif

#endif

static void

audio_resampler_init (void)

{

  static gsize init_gonce = 0;

  if (g_once_init_enter (&init_gonce)) {

    GST_DEBUG_CATEGORY_INIT (audio_resampler_debug, "audio-resampler", 0,

        "audio-resampler object");

#if defined HAVE_ORC && !defined DISABLE_ORC

    orc_init ();

    {

      OrcTarget *target = orc_target_get_default ();

      gint i;

      if (target) {

        const gchar *name;

        unsigned int flags = orc_target_get_default_flags (target);

        for (i = -1; i < 32; ++i) {

          if (i == -1) {

            name = orc_target_get_name (target);

            GST_DEBUG ("target %s, default flags %08x", name, flags);

          } else if (flags & (1U << i)) {

            name = orc_target_get_flag_name (target, i);

            GST_DEBUG ("target flag %s", name);

          } else

            name = NULL;

          if (name) {

#ifdef CHECK_X86

            audio_resampler_check_x86 (name);

#endif

#ifdef CHECK_NEON

            audio_resampler_check_neon (name);

#endif

          }

        }

      }

    }

#endif

    g_once_init_leave (&init_gonce, 1);

  }