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;
62 aubio_dct_t * new_aubio_dct (uint_t size) {
63 aubio_dct_t * s = AUBIO_NEW(aubio_dct_t);
68 s->in = new_fvec(size);
69 s->out = new_fvec(size);
70 pthread_mutex_lock(&aubio_fftw_mutex);
71 s->data = (smpl_t *)fftw_malloc(sizeof(smpl_t) * size);
72 s->pfw = fftw_plan_r2r_1d(size, s->in->data, s->data, FFTW_REDFT10,
74 s->pbw = fftw_plan_r2r_1d(size, s->data, s->out->data, FFTW_REDFT01,
76 pthread_mutex_unlock(&aubio_fftw_mutex);
77 s->scalers[0] = SQRT(1./(4.*s->size));
78 s->scalers[1] = SQRT(1./(2.*s->size));
79 s->scalers[2] = 1. / s->scalers[0];
80 s->scalers[3] = 1. / s->scalers[1];
81 s->scalers[4] = .5 / s->size;
88 void del_aubio_dct(aubio_dct_t *s) {
89 pthread_mutex_lock(&aubio_fftw_mutex);
90 fftw_destroy_plan(s->pfw);
91 fftw_destroy_plan(s->pbw);
93 pthread_mutex_unlock(&aubio_fftw_mutex);
99 void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
101 fvec_copy(input, s->in);
102 fftw_execute(s->pfw);
103 //fvec_copy(s->out, output);
104 s->data[0] *= s->scalers[0];
105 for (i = 1; i < s->size; i++) {
106 s->data[i] *= s->scalers[1];
108 memcpy(output->data, s->data, output->length * sizeof(smpl_t));
111 void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
113 memcpy(s->data, input->data, input->length * sizeof(smpl_t));
115 s->data[0] *= s->scalers[2];
116 for (i = 1; i < s->size; i++) {
117 s->data[i] *= s->scalers[3];
119 fftw_execute(s->pbw);
120 for (i = 0; i < s->size; i++) {
121 s->out->data[i] *= s->scalers[4];
123 fvec_copy(s->out, output);