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"
24 #include "spectral/dct.h"
33 #error "Using aubio in double precision with fftw3 in single precision"
34 #endif /* HAVE_AUBIO_DOUBLE */
35 #else /* HAVE_FFTW3F */
36 #if !HAVE_AUBIO_DOUBLE
37 #error "Using aubio in single precision with fftw3 in double precision"
38 #endif /* HAVE_AUBIO_DOUBLE */
39 #endif /* HAVE_FFTW3F */
42 #define fftw_malloc fftwf_malloc
43 #define fftw_free fftwf_free
44 #define fftw_execute fftwf_execute
45 #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
46 #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
47 #define fftw_plan_r2r_1d fftwf_plan_r2r_1d
48 #define fftw_plan fftwf_plan
49 #define fftw_destroy_plan fftwf_destroy_plan
52 // defined in src/spectral/fft.c
53 extern pthread_mutex_t aubio_fftw_mutex;
55 extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *);
65 aubio_dct_t * new_aubio_dct (uint_t size) {
66 aubio_dct_t * s = AUBIO_NEW(aubio_dct_t);
71 s->in = new_fvec(size);
72 s->out = new_fvec(size);
73 pthread_mutex_lock(&aubio_fftw_mutex);
74 s->data = (smpl_t *)fftw_malloc(sizeof(smpl_t) * size);
75 s->pfw = fftw_plan_r2r_1d(size, s->in->data, s->data, FFTW_REDFT10,
77 s->pbw = fftw_plan_r2r_1d(size, s->data, s->out->data, FFTW_REDFT01,
79 pthread_mutex_unlock(&aubio_fftw_mutex);
80 s->scalers[0] = SQRT(1./(4.*s->size));
81 s->scalers[1] = SQRT(1./(2.*s->size));
82 s->scalers[2] = 1. / s->scalers[0];
83 s->scalers[3] = 1. / s->scalers[1];
84 s->scalers[4] = .5 / s->size;
91 void del_aubio_dct(aubio_dct_t *s) {
92 pthread_mutex_lock(&aubio_fftw_mutex);
93 fftw_destroy_plan(s->pfw);
94 fftw_destroy_plan(s->pbw);
96 pthread_mutex_unlock(&aubio_fftw_mutex);
102 void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
104 fvec_copy(input, s->in);
105 fftw_execute(s->pfw);
106 //fvec_copy(s->out, output);
107 s->data[0] *= s->scalers[0];
108 for (i = 1; i < s->size; i++) {
109 s->data[i] *= s->scalers[1];
111 memcpy(output->data, s->data, output->length * sizeof(smpl_t));
114 void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) {
116 memcpy(s->data, input->data, input->length * sizeof(smpl_t));
118 s->data[0] *= s->scalers[2];
119 for (i = 1; i < s->size; i++) {
120 s->data[i] *= s->scalers[3];
122 fftw_execute(s->pbw);
123 for (i = 0; i < s->size; i++) {
124 s->out->data[i] *= s->scalers[4];
126 fvec_copy(s->out, output);