src/spectral/dct_ooura.c

   1 /*
   2   Copyright (C) 2017 Paul Brossier <piem@aubio.org>
   3
   4   This file is part of aubio.
   5
   6   aubio is free software: you can redistribute it and/or modify
   7   it under the terms of the GNU General Public License as published by
   8   the Free Software Foundation, either version 3 of the License, or
   9   (at your option) any later version.
  10
  11   aubio is distributed in the hope that it will be useful,
  12   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14   GNU General Public License for more details.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with aubio.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 */
  20
  21 #include "aubio_priv.h"
  22 #include "fvec.h"
  23 #include "cvec.h"
  24 #include "spectral/dct.h"
  25
  26 #if !defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3) && !defined(HAVE_INTEL_IPP)
  27
  28 extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *);
  29
  30 struct _aubio_dct_t {
  31   uint_t size;
  32   fvec_t *input;
  33   smpl_t *w;
  34   int *ip;
  35   smpl_t scalers[5];
  36 };
  37
  38 aubio_dct_t * new_aubio_dct (uint_t size) {
  39   aubio_dct_t * s = AUBIO_NEW(aubio_dct_t);
  40   if (aubio_is_power_of_two(size) != 1) {
  41     AUBIO_ERR("dct: can only create with sizes power of two, requested %d\n",
  42         size);
  43     goto beach;
  44   }
  45   s->size = size;
  46   s->input = new_fvec(s->size);
  47   s->w = AUBIO_ARRAY(smpl_t, s->size * 5 / 4);
  48   s->ip = AUBIO_ARRAY(int, 3 + (1 << (int)FLOOR(LOG(s->size/2) / LOG(2))) / 2);
  49   s->ip[0] = 0;
  50   s->scalers[0] = 2. * SQRT(1./(4.*s->size));
  51   s->scalers[1] = 2. * SQRT(1./(2.*s->size));
  52   s->scalers[2] = 1. / s->scalers[0];
  53   s->scalers[3] = 1. / s->scalers[1];
  54   s->scalers[4] = 2. / s->size;
  55   return s;
  56 beach:
  57   AUBIO_FREE(s);
  58   return NULL;
  59 }
  60
  61 void del_aubio_dct(aubio_dct_t *s) {
  62   del_fvec(s->input);
  63   AUBIO_FREE(s->ip);
  64   AUBIO_FREE(s->w);
  65   AUBIO_FREE(s);
  66 }
  67
  68 void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
  69   uint_t i = 0;
  70   fvec_copy(input, s->input);
  71   aubio_ooura_ddct(s->size, -1, s->input->data, s->ip, s->w);
  72   // apply orthonormal scaling
  73   s->input->data[0] *= s->scalers[0];
  74   for (i = 1; i < s->input->length; i++) {
  75     s->input->data[i] *= s->scalers[1];
  76   }
  77   fvec_copy(s->input, output);
  78 }
  79
  80 void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
  81   uint_t i = 0;
  82   fvec_copy(input, s->input);
  83   s->input->data[0] *= s->scalers[2];
  84   for (i = 1; i < s->input->length; i++) {
  85     s->input->data[i] *= s->scalers[3];
  86   }
  87   s->input->data[0] *= .5;
  88   aubio_ooura_ddct(s->size, 1, s->input->data, s->ip, s->w);
  89   for (i = 0; i < s->input->length; i++) {
  90     s->input->data[i] *= s->scalers[4];
  91   }
  92   fvec_copy(s->input, output);
  93 }
  94
  95 #endif //!defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3)