2 Copyright (C) 2017 Paul Brossier <piem@aubio.org>
4 This file is part of aubio.
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.
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.
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/>.
21 #include "aubio_priv.h"
23 #include "spectral/dct.h"
32 #error "Using aubio in double precision with fftw3 in single precision"
33 #endif /* HAVE_AUBIO_DOUBLE */
34 #else /* HAVE_FFTW3F */
35 #if !HAVE_AUBIO_DOUBLE
36 #error "Using aubio in single precision with fftw3 in double precision"
37 #endif /* HAVE_AUBIO_DOUBLE */
38 #endif /* HAVE_FFTW3F */
41 #define fftw_malloc fftwf_malloc
42 #define fftw_free fftwf_free
43 #define fftw_execute fftwf_execute
44 #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
45 #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
46 #define fftw_plan_r2r_1d fftwf_plan_r2r_1d
47 #define fftw_plan fftwf_plan
48 #define fftw_destroy_plan fftwf_destroy_plan
51 // defined in src/spectral/fft.c
52 extern pthread_mutex_t aubio_fftw_mutex;
54 extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *);
64 aubio_dct_t * new_aubio_dct (uint_t size) {
65 aubio_dct_t * s = AUBIO_NEW(aubio_dct_t);
70 s->in = new_fvec(size);
71 s->out = new_fvec(size);
72 pthread_mutex_lock(&aubio_fftw_mutex);
73 s->data = (smpl_t *)fftw_malloc(sizeof(smpl_t) * size);
74 s->pfw = fftw_plan_r2r_1d(size, s->in->data, s->data, FFTW_REDFT10,
76 s->pbw = fftw_plan_r2r_1d(size, s->data, s->out->data, FFTW_REDFT01,
78 pthread_mutex_unlock(&aubio_fftw_mutex);
79 s->scalers[0] = SQRT(1./(4.*s->size));
80 s->scalers[1] = SQRT(1./(2.*s->size));
81 s->scalers[2] = 1. / s->scalers[0];
82 s->scalers[3] = 1. / s->scalers[1];
83 s->scalers[4] = .5 / s->size;
90 void del_aubio_dct(aubio_dct_t *s) {
91 pthread_mutex_lock(&aubio_fftw_mutex);
92 fftw_destroy_plan(s->pfw);
93 fftw_destroy_plan(s->pbw);
95 pthread_mutex_unlock(&aubio_fftw_mutex);
101 void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
103 fvec_copy(input, s->in);
104 fftw_execute(s->pfw);
105 //fvec_copy(s->out, output);
106 s->data[0] *= s->scalers[0];
107 for (i = 1; i < s->size; i++) {
108 s->data[i] *= s->scalers[1];
110 memcpy(output->data, s->data, output->length * sizeof(smpl_t));
113 void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
115 memcpy(s->data, input->data, input->length * sizeof(smpl_t));
117 s->data[0] *= s->scalers[2];
118 for (i = 1; i < s->size; i++) {
119 s->data[i] *= s->scalers[3];
121 fftw_execute(s->pbw);
122 for (i = 0; i < s->size; i++) {
123 s->out->data[i] *= s->scalers[4];
125 fvec_copy(s->out, output);