[4daa7ec] | 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) |
---|
[0fbbb067] | 43 | typedef struct _aubio_dct_accelerate_t aubio_dct_accelerate_t; |
---|
[4daa7ec] | 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) |
---|
[0fbbb067] | 49 | typedef struct _aubio_dct_fftw_t aubio_dct_fftw_t; |
---|
[4daa7ec] | 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 defined(HAVE_ACCELERATE) |
---|
| 85 | // vDSP supports sizes = f * 2 ** n, where n >= 4 and f in [1, 3, 5, 15] |
---|
| 86 | // see https://developer.apple.com/documentation/accelerate/1449930-vdsp_dct_createsetup |
---|
[8b69453] | 87 | { |
---|
| 88 | uint_t radix = size; |
---|
| 89 | uint_t order = 0; |
---|
[e680926] | 90 | while ((radix >= 1) && ((radix / 2) * 2 == radix)) { |
---|
[8b69453] | 91 | radix /= 2; |
---|
| 92 | order++; |
---|
| 93 | } |
---|
| 94 | if (order < 4 || (radix != 1 && radix != 3 && radix != 5 && radix != 15)) { |
---|
| 95 | goto plain; |
---|
| 96 | } |
---|
[4daa7ec] | 97 | } |
---|
| 98 | s->dct = (void *)new_aubio_dct_accelerate (size); |
---|
| 99 | if (s->dct) { |
---|
| 100 | s->dct_do = (aubio_dct_do_t)aubio_dct_accelerate_do; |
---|
| 101 | s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_accelerate_rdo; |
---|
| 102 | s->del_dct = (del_aubio_dct_t)del_aubio_dct_accelerate; |
---|
| 103 | return s; |
---|
| 104 | } |
---|
| 105 | #elif defined(HAVE_FFTW3) |
---|
| 106 | // fftw supports any positive integer size |
---|
| 107 | s->dct = (void *)new_aubio_dct_fftw (size); |
---|
| 108 | if (s->dct) { |
---|
| 109 | s->dct_do = (aubio_dct_do_t)aubio_dct_fftw_do; |
---|
| 110 | s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_fftw_rdo; |
---|
| 111 | s->del_dct = (del_aubio_dct_t)del_aubio_dct_fftw; |
---|
| 112 | return s; |
---|
| 113 | } else { |
---|
[5c3061e] | 114 | AUBIO_WRN("dct: unexpected error while creating dct_fftw with size %d\n", |
---|
[4daa7ec] | 115 | size); |
---|
| 116 | goto plain; |
---|
| 117 | } |
---|
| 118 | #elif defined(HAVE_INTEL_IPP) |
---|
| 119 | // unclear from the docs, but intel ipp seems to support any size |
---|
| 120 | s->dct = (void *)new_aubio_dct_ipp (size); |
---|
| 121 | if (s->dct) { |
---|
| 122 | s->dct_do = (aubio_dct_do_t)aubio_dct_ipp_do; |
---|
| 123 | s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_ipp_rdo; |
---|
| 124 | s->del_dct = (del_aubio_dct_t)del_aubio_dct_ipp; |
---|
| 125 | return s; |
---|
| 126 | } else { |
---|
[5c3061e] | 127 | AUBIO_WRN("dct: unexpected error while creating dct_ipp with size %d\n", |
---|
[4daa7ec] | 128 | size); |
---|
| 129 | goto plain; |
---|
| 130 | } |
---|
| 131 | #else |
---|
[e458af4] | 132 | // ooura support sizes that are power of 2 |
---|
[aad7702] | 133 | if (aubio_is_power_of_two(size) != 1 || size == 1) { |
---|
[4daa7ec] | 134 | goto plain; |
---|
| 135 | } |
---|
| 136 | s->dct = (void *)new_aubio_dct_ooura (size); |
---|
| 137 | if (s->dct) { |
---|
| 138 | s->dct_do = (aubio_dct_do_t)aubio_dct_ooura_do; |
---|
| 139 | s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_ooura_rdo; |
---|
| 140 | s->del_dct = (del_aubio_dct_t)del_aubio_dct_ooura; |
---|
| 141 | return s; |
---|
| 142 | } |
---|
| 143 | #endif |
---|
| 144 | // falling back to plain mode |
---|
[d3440e7] | 145 | AUBIO_WRN("dct: no optimised implementation could be created for size %d\n", |
---|
[4daa7ec] | 146 | size); |
---|
| 147 | plain: |
---|
| 148 | s->dct = (void *)new_aubio_dct_plain (size); |
---|
| 149 | if (s->dct) { |
---|
| 150 | s->dct_do = (aubio_dct_do_t)aubio_dct_plain_do; |
---|
| 151 | s->dct_rdo = (aubio_dct_rdo_t)aubio_dct_plain_rdo; |
---|
| 152 | s->del_dct = (del_aubio_dct_t)del_aubio_dct_plain; |
---|
| 153 | return s; |
---|
| 154 | } else { |
---|
| 155 | goto beach; |
---|
| 156 | } |
---|
| 157 | beach: |
---|
[7c85c15] | 158 | AUBIO_ERROR("dct: failed creating with size %d, should be > 0\n", size); |
---|
| 159 | del_aubio_dct(s); |
---|
[4daa7ec] | 160 | return NULL; |
---|
| 161 | } |
---|
| 162 | |
---|
| 163 | void del_aubio_dct(aubio_dct_t *s) { |
---|
| 164 | if (s->dct && s->del_dct) s->del_dct (s->dct); |
---|
| 165 | AUBIO_FREE (s); |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { |
---|
| 169 | s->dct_do ((void *)s->dct, input, output); |
---|
| 170 | } |
---|
| 171 | |
---|
| 172 | void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { |
---|
| 173 | s->dct_rdo ((void *)s->dct, input, output); |
---|
| 174 | } |
---|