source: src/spectral/dct.c @ e2010b3

feature/crepe_org
Last change on this file since e2010b3 was e680926, checked in by Paul Brossier <piem@piem.org>, 6 years ago

[dct] fix radix check for accelerate

  • Property mode set to 100644
File size: 6.1 KB
RevLine 
[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
38typedef void (*aubio_dct_do_t)(aubio_dct_t * s, const fvec_t * input, fvec_t * output);
39typedef void (*aubio_dct_rdo_t)(aubio_dct_t * s, const fvec_t * input, fvec_t * output);
40typedef void (*del_aubio_dct_t)(aubio_dct_t * s);
41
42#if defined(HAVE_ACCELERATE)
[0fbbb067]43typedef struct _aubio_dct_accelerate_t aubio_dct_accelerate_t;
[4daa7ec]44extern aubio_dct_accelerate_t * new_aubio_dct_accelerate (uint_t size);
45extern void aubio_dct_accelerate_do(aubio_dct_accelerate_t *s, const fvec_t *input, fvec_t *output);
46extern void aubio_dct_accelerate_rdo(aubio_dct_accelerate_t *s, const fvec_t *input, fvec_t *output);
47extern void del_aubio_dct_accelerate (aubio_dct_accelerate_t *s);
48#elif defined(HAVE_FFTW3)
[0fbbb067]49typedef struct _aubio_dct_fftw_t aubio_dct_fftw_t;
[4daa7ec]50extern aubio_dct_fftw_t * new_aubio_dct_fftw (uint_t size);
51extern void aubio_dct_fftw_do(aubio_dct_fftw_t *s, const fvec_t *input, fvec_t *output);
52extern void aubio_dct_fftw_rdo(aubio_dct_fftw_t *s, const fvec_t *input, fvec_t *output);
53extern void del_aubio_dct_fftw (aubio_dct_fftw_t *s);
54#elif defined(HAVE_INTEL_IPP)
55typedef struct _aubio_dct_ipp_t aubio_dct_ipp_t;
56extern aubio_dct_ipp_t * new_aubio_dct_ipp (uint_t size);
57extern void aubio_dct_ipp_do(aubio_dct_ipp_t *s, const fvec_t *input, fvec_t *output);
58extern void aubio_dct_ipp_rdo(aubio_dct_ipp_t *s, const fvec_t *input, fvec_t *output);
59extern void del_aubio_dct_ipp (aubio_dct_ipp_t *s);
60#else
61typedef struct _aubio_dct_ooura_t aubio_dct_ooura_t;
62extern aubio_dct_ooura_t * new_aubio_dct_ooura (uint_t size);
63extern void aubio_dct_ooura_do(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output);
64extern void aubio_dct_ooura_rdo(aubio_dct_ooura_t *s, const fvec_t *input, fvec_t *output);
65extern void del_aubio_dct_ooura (aubio_dct_ooura_t *s);
66#endif
67
68// plain mode
69typedef struct _aubio_dct_plain_t aubio_dct_plain_t;
70extern aubio_dct_plain_t * new_aubio_dct_plain (uint_t size);
71extern void aubio_dct_plain_do(aubio_dct_plain_t *s, const fvec_t *input, fvec_t *output);
72extern void aubio_dct_plain_rdo(aubio_dct_plain_t *s, const fvec_t *input, fvec_t *output);
73extern void del_aubio_dct_plain (aubio_dct_plain_t *s);
74
75struct _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
82aubio_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);
147plain:
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  }
157beach:
[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
163void 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
168void 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
172void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
173  s->dct_rdo ((void *)s->dct, input, output);
174}
Note: See TracBrowser for help on using the repository browser.