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 "spectral/dct.h"
  24
  25 #if !defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3) && !defined(HAVE_INTEL_IPP)
  26
  27 typedef struct _aubio_dct_ooura_t aubio_dct_ooura_t;
  28
  29 extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *);
  30
  31 struct _aubio_dct_ooura_t {
  32   uint_t size;
  33   fvec_t *input;
  34   smpl_t *w;
  35   int *ip;
  36   smpl_t scalers[5];
  37 };
  38
  39 aubio_dct_ooura_t * new_aubio_dct_ooura (uint_t size) {
  40   aubio_dct_ooura_t * s = AUBIO_NEW(aubio_dct_ooura_t);
  41   if (aubio_is_power_of_two(size) != 1 || (sint_t)size <= 0) {
  42     AUBIO_ERR("dct_ooura: can only create with sizes power of two, requested %d\n",
  43         size);
  44     goto beach;
  45   }
  46   s->size = size;
  47   s->input = new_fvec(s->size);
  48   s->w = AUBIO_ARRAY(smpl_t, s->size * 5 / 4);
  49   s->ip = AUBIO_ARRAY(int, 3 + (1 << (int)FLOOR(LOG(s->size/2) / LOG(2))) / 2);
  50   s->ip[0] = 0;
  51   s->scalers[0] = 2. * SQRT(1./(4.*s->size));
  52   s->scalers[1] = 2. * SQRT(1./(2.*s->size));
  53   s->scalers[2] = 1. / s->scalers[0];
  54   s->scalers[3] = 1. / s->scalers[1];
  55   s->scalers[4] = 2. / s->size;
  56   return s;
  57 beach:
  58   AUBIO_FREE(s);
  59   return NULL;
  60 }
  61
  62 void del_aubio_dct_ooura(aubio_dct_ooura_t *s) {
  63   del_fvec(s->input);
  64   AUBIO_FREE(s->ip);
  65   AUBIO_FREE(s->w);
  66   AUBIO_FREE(s);
  67 }
  68
  69 void aubio_dct_ooura_do(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output) {
  70   uint_t i = 0;
  71   fvec_copy(input, s->input);
  72   aubio_ooura_ddct(s->size, -1, s->input->data, s->ip, s->w);
  73   // apply orthonormal scaling
  74   s->input->data[0] *= s->scalers[0];
  75   for (i = 1; i < s->input->length; i++) {
  76     s->input->data[i] *= s->scalers[1];
  77   }
  78   fvec_copy(s->input, output);
  79 }
  80
  81 void aubio_dct_ooura_rdo(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output) {
  82   uint_t i = 0;
  83   fvec_copy(input, s->input);
  84   s->input->data[0] *= s->scalers[2];
  85   for (i = 1; i < s->input->length; i++) {
  86     s->input->data[i] *= s->scalers[3];
  87   }
  88   s->input->data[0] *= .5;
  89   aubio_ooura_ddct(s->size, 1, s->input->data, s->ip, s->w);
  90   for (i = 0; i < s->input->length; i++) {
  91     s->input->data[i] *= s->scalers[4];
  92   }
  93   fvec_copy(s->input, output);
  94 }
  95
  96 #endif //!defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3)