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