src/spectral/dct.c

   1 /*
   2   Copyright (C) 2018 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 /** \file
  22
  23   Discrete Cosine Transform
  24
  25   Functions aubio_dct_do() and aubio_dct_rdo() are equivalent to MATLAB/Octave
  26   dct() and idct() functions, as well as scipy.fftpack.dct(x, norm='ortho') and
  27   scipy.fftpack.idct(x, norm='ortho')
  28
  29   \example spectral/test-dct.c
  30
  31 */
  32
  33 #include "aubio_priv.h"
  34 #include "fvec.h"
  35 #include "spectral/dct.h"
  36
  37 // function pointers prototypes
  38 typedef void (*aubio_dct_do_t)(aubio_dct_t * s, const fvec_t * input, fvec_t * output);
  39 typedef void (*aubio_dct_rdo_t)(aubio_dct_t * s, const fvec_t * input, fvec_t * output);
  40 typedef void (*del_aubio_dct_t)(aubio_dct_t * s);
  41
  42 #if defined(HAVE_ACCELERATE)
  43 typedef struct _aubio_dct_accelerate_t aubio_dct_accelerate_t;
  44 extern aubio_dct_accelerate_t * new_aubio_dct_accelerate (uint_t size);
  45 extern void aubio_dct_accelerate_do(aubio_dct_accelerate_t *s, const fvec_t *input, fvec_t *output);
  46 extern void aubio_dct_accelerate_rdo(aubio_dct_accelerate_t *s, const fvec_t *input, fvec_t *output);
  47 extern void del_aubio_dct_accelerate (aubio_dct_accelerate_t *s);
  48 #elif defined(HAVE_FFTW3)
  49 typedef struct _aubio_dct_fftw_t aubio_dct_fftw_t;
  50 extern aubio_dct_fftw_t * new_aubio_dct_fftw (uint_t size);
  51 extern void aubio_dct_fftw_do(aubio_dct_fftw_t *s, const fvec_t *input, fvec_t *output);
  52 extern void aubio_dct_fftw_rdo(aubio_dct_fftw_t *s, const fvec_t *input, fvec_t *output);
  53 extern void del_aubio_dct_fftw (aubio_dct_fftw_t *s);
  54 #elif defined(HAVE_INTEL_IPP)
  55 typedef struct _aubio_dct_ipp_t aubio_dct_ipp_t;
  56 extern aubio_dct_ipp_t * new_aubio_dct_ipp (uint_t size);
  57 extern void aubio_dct_ipp_do(aubio_dct_ipp_t *s, const fvec_t *input, fvec_t *output);
  58 extern void aubio_dct_ipp_rdo(aubio_dct_ipp_t *s, const fvec_t *input, fvec_t *output);
  59 extern void del_aubio_dct_ipp (aubio_dct_ipp_t *s);
  60 #else
  61 typedef struct _aubio_dct_ooura_t aubio_dct_ooura_t;
  62 extern aubio_dct_ooura_t * new_aubio_dct_ooura (uint_t size);
  63 extern void aubio_dct_ooura_do(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output);
  64 extern void aubio_dct_ooura_rdo(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output);
  65 extern void del_aubio_dct_ooura (aubio_dct_ooura_t *s);
  66 #endif
  67
  68 // plain mode
  69 typedef struct _aubio_dct_plain_t aubio_dct_plain_t;
  70 extern aubio_dct_plain_t * new_aubio_dct_plain (uint_t size);
  71 extern void aubio_dct_plain_do(aubio_dct_plain_t *s, const fvec_t *input, fvec_t *output);
  72 extern void aubio_dct_plain_rdo(aubio_dct_plain_t *s, const fvec_t *input, fvec_t *output);
  73 extern void del_aubio_dct_plain (aubio_dct_plain_t *s);
  74
  75 struct _aubio_dct_t {
  76   void *dct;
  77   aubio_dct_do_t dct_do;
  78   aubio_dct_rdo_t dct_rdo;
  79   del_aubio_dct_t del_dct;
  80 };
  81
  82 aubio_dct_t* new_aubio_dct (uint_t size) {
  83   aubio_dct_t * s = AUBIO_NEW(aubio_dct_t);
  84   if ((sint_t)size <= 0) goto beach;
  85 #if defined(HAVE_ACCELERATE)
  86   // vDSP supports sizes = f * 2 ** n, where n >= 4 and f in [1, 3, 5, 15]
  87   // see https://developer.apple.com/documentation/accelerate/1449930-vdsp_dct_createsetup
  88   {
  89     uint_t radix = size;
  90     uint_t order = 0;
  91     while ((radix / 2) * 2 == radix) {
  92       radix /= 2;
  93       order++;
  94     }
  95     if (order < 4 || (radix != 1 && radix != 3 && radix != 5 && radix != 15)) {
  96       goto plain;
  97     }
  98   }
  99   s->dct = (void *)new_aubio_dct_accelerate (size);
 100   if (s->dct) {
 101     s->dct_do = (aubio_dct_do_t)aubio_dct_accelerate_do;
 102     s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_accelerate_rdo;
 103     s->del_dct = (del_aubio_dct_t)del_aubio_dct_accelerate;
 104     return s;
 105   }
 106 #elif defined(HAVE_FFTW3)
 107   // fftw supports any positive integer size
 108   s->dct = (void *)new_aubio_dct_fftw (size);
 109   if (s->dct) {
 110     s->dct_do = (aubio_dct_do_t)aubio_dct_fftw_do;
 111     s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_fftw_rdo;
 112     s->del_dct = (del_aubio_dct_t)del_aubio_dct_fftw;
 113     return s;
 114   } else {
 115     AUBIO_WRN("dct: unexcepected error while creating dct_fftw with size %d",
 116         size);
 117     goto plain;
 118   }
 119 #elif defined(HAVE_INTEL_IPP)
 120   // unclear from the docs, but intel ipp seems to support any size
 121   s->dct = (void *)new_aubio_dct_ipp (size);
 122   if (s->dct) {
 123     s->dct_do = (aubio_dct_do_t)aubio_dct_ipp_do;
 124     s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_ipp_rdo;
 125     s->del_dct = (del_aubio_dct_t)del_aubio_dct_ipp;
 126     return s;
 127   } else {
 128     AUBIO_WRN("dct: unexcepected error while creating dct_ipp with size %d",
 129         size);
 130     goto plain;
 131   }
 132 #else
 133   // ooura support sizes that are power of 2
 134   if (aubio_is_power_of_two(size) != 1 || size == 1) {
 135     goto plain;
 136   }
 137   s->dct = (void *)new_aubio_dct_ooura (size);
 138   if (s->dct) {
 139     s->dct_do = (aubio_dct_do_t)aubio_dct_ooura_do;
 140     s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_ooura_rdo;
 141     s->del_dct = (del_aubio_dct_t)del_aubio_dct_ooura;
 142     return s;
 143   }
 144 #endif
 145   // falling back to plain mode
 146   AUBIO_WRN("dct: d no optimised implementation could be created for size %d",
 147       size);
 148 plain:
 149   s->dct = (void *)new_aubio_dct_plain (size);
 150   if (s->dct) {
 151     s->dct_do = (aubio_dct_do_t)aubio_dct_plain_do;
 152     s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_plain_rdo;
 153     s->del_dct = (del_aubio_dct_t)del_aubio_dct_plain;
 154     return s;
 155   } else {
 156     goto beach;
 157   }
 158 beach:
 159   AUBIO_ERROR("dct: failed creating with size %d, should be > 0", size);
 160   AUBIO_FREE (s);
 161   return NULL;
 162 }
 163
 164 void del_aubio_dct(aubio_dct_t *s) {
 165   if (s->dct && s->del_dct) s->del_dct (s->dct);
 166   AUBIO_FREE (s);
 167 }
 168
 169 void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
 170   s->dct_do ((void *)s->dct, input, output);
 171 }
 172
 173 void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
 174   s->dct_rdo ((void *)s->dct, input, output);
 175 }