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