From 60583a3817b8c4120861ec7df2e9836146efd54d Mon Sep 17 00:00:00 2001 From: Paul Brossier Date: Tue, 5 Sep 2017 12:56:39 +0200 Subject: [PATCH] src/spectral/dct.h: add dct type II using ooura --- python/lib/gen_code.py | 2 + src/aubio.h | 1 + src/spectral/dct.h | 85 +++++++++++++++++++++++++++++++++++++++++++ src/spectral/dct_ooura.c | 85 +++++++++++++++++++++++++++++++++++++++++++ tests/src/spectral/test-dct.c | 32 ++++++++++++++++ 5 files changed, 205 insertions(+) create mode 100644 src/spectral/dct.h create mode 100644 src/spectral/dct_ooura.c create mode 100644 tests/src/spectral/test-dct.c diff --git a/python/lib/gen_code.py b/python/lib/gen_code.py index d3911974..0498aa20 100644 --- a/python/lib/gen_code.py +++ b/python/lib/gen_code.py @@ -2,6 +2,7 @@ aubiodefvalue = { # we have some clean up to do 'buf_size': 'Py_default_vector_length', 'win_s': 'Py_default_vector_length', + 'size': 'Py_default_vector_length', # and here too 'hop_size': 'Py_default_vector_length / 2', 'hop_s': 'Py_default_vector_length / 2', @@ -82,6 +83,7 @@ objoutsize = { 'tempo': '1', 'filterbank': 'self->n_filters', 'tss': 'self->buf_size', + 'dct': 'self->size', } objinputsize = { diff --git a/src/aubio.h b/src/aubio.h index 0d021361..73deb24c 100644 --- a/src/aubio.h +++ b/src/aubio.h @@ -182,6 +182,7 @@ extern "C" #include "temporal/a_weighting.h" #include "temporal/c_weighting.h" #include "spectral/fft.h" +#include "spectral/dct.h" #include "spectral/phasevoc.h" #include "spectral/filterbank.h" #include "spectral/filterbank_mel.h" diff --git a/src/spectral/dct.h b/src/spectral/dct.h new file mode 100644 index 00000000..1da1a0d5 --- /dev/null +++ b/src/spectral/dct.h @@ -0,0 +1,85 @@ +/* + Copyright (C) 2017 Paul Brossier + + This file is part of aubio. + + aubio is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + aubio is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with aubio. If not, see . + +*/ + +/** \file + + Discrete Cosine Transform + + Functions aubio_dct_do() and aubio_dct_rdo() are equivalent to MATLAB/Octave + dct() and idct() functions, as well as scipy.fftpack.dct(x, norm='ortho') and + scipy.fftpack.idct(x, norm='ortho') + + \example spectral/test-dct.c + +*/ + +#ifndef AUBIO_DCT_H +#define AUBIO_DCT_H + +#ifdef __cplusplus +extern "C" { +#endif + +/** DCT object + + This object computes forward and backward DCT type 2 with orthonormal + scaling. + +*/ +typedef struct _aubio_dct_t aubio_dct_t; + +/** create new DCT computation object + + \param size length of the DCT + +*/ +aubio_dct_t * new_aubio_dct(uint_t size); + +/** compute forward DCT + + \param s dct object as returned by new_aubio_dct + \param input input signal + \param output transformed input array + +*/ +void aubio_dct_do (aubio_dct_t *s, const fvec_t * input, fvec_t * output); + +/** compute backward DCT + + \param s dct object as returned by new_aubio_dct + \param input input signal + \param output transformed input array + +*/ +void aubio_dct_rdo (aubio_dct_t *s, const fvec_t * input, fvec_t * output); + + +/** delete DCT object + + \param s dct object as returned by new_aubio_dct + +*/ +void del_aubio_dct (aubio_dct_t *s); + +#ifdef __cplusplus +} +#endif + +#endif /* AUBIO_DCT_H */ diff --git a/src/spectral/dct_ooura.c b/src/spectral/dct_ooura.c new file mode 100644 index 00000000..1fafc0e0 --- /dev/null +++ b/src/spectral/dct_ooura.c @@ -0,0 +1,85 @@ +/* + Copyright (C) 2017 Paul Brossier + + This file is part of aubio. + + aubio is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + aubio is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with aubio. If not, see . + +*/ + +#include "aubio_priv.h" +#include "fvec.h" +#include "cvec.h" +#include "spectral/dct.h" + +extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *); + +struct _aubio_dct_t { + uint_t size; + fvec_t *input; + smpl_t *w; + int *ip; +}; + +aubio_dct_t * new_aubio_dct (uint_t size) { + aubio_dct_t * s = AUBIO_NEW(aubio_dct_t); + if (aubio_is_power_of_two(size) != 1) { + AUBIO_ERR("dct: can only create with sizes power of two, requested %d\n", + size); + goto beach; + } + s->size = size; + s->input = new_fvec(s->size); + s->w = AUBIO_ARRAY(smpl_t, s->size * 5 / 4); + s->ip = AUBIO_ARRAY(int, 3 + (1 << (int)FLOOR(LOG(s->size/2) / LOG(2))) / 2); + s->ip[0] = 0; + return s; +beach: + AUBIO_FREE(s); + return NULL; +} + +void del_aubio_dct(aubio_dct_t *s) { + del_fvec(s->input); + AUBIO_FREE(s->ip); + AUBIO_FREE(s->w); + AUBIO_FREE(s); +} + +void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + uint_t i = 0; + fvec_copy(input, s->input); + aubio_ooura_ddct(s->size, -1, s->input->data, s->ip, s->w); + // apply orthonormal scaling + s->input->data[0] *= 2.* SQRT(1./(4.*s->input->length)); + for (i = 1; i < s->input->length; i++) { + s->input->data[i] *= 2. * SQRT(1./(2.*s->input->length)); + } + fvec_copy(s->input, output); +} + +void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + uint_t i = 0; + fvec_copy(input, s->input); + s->input->data[0] /= 2.* SQRT(1./(4.*s->input->length)); + for (i = 1; i < s->input->length; i++) { + s->input->data[i] /= 2. * SQRT(1./(2.*s->input->length)); + } + s->input->data[0] *= .5; + aubio_ooura_ddct(s->size, 1, s->input->data, s->ip, s->w); + for (i = 0; i < s->input->length; i++) { + s->input->data[i] *= 2./s->input->length; + } + fvec_copy(s->input, output); +} diff --git a/tests/src/spectral/test-dct.c b/tests/src/spectral/test-dct.c new file mode 100644 index 00000000..dc3a9149 --- /dev/null +++ b/tests/src/spectral/test-dct.c @@ -0,0 +1,32 @@ +#include + +int main (void) +{ + int return_code = 0; + uint_t win_s = 32; // window size + uint_t i, n_iters = 10; // number of iterations + // create dct object + aubio_dct_t * dct = new_aubio_dct(win_s); + + fvec_t * in = new_fvec (win_s); // input buffer + fvec_t * dctout = new_fvec (win_s); // output buffer + + if (!dct || !in || !dctout) { + return_code = 1; + return return_code; + } + + in->data[0] = 1.; + for (i = 0; i < n_iters; i++) { + // execute stft + aubio_dct_do (dct, in, dctout); + aubio_dct_rdo (dct, dctout, in); + } + fvec_print(dctout); + fvec_print(in); + del_fvec(dctout); + del_fvec(in); + + del_aubio_dct(dct); + return return_code; +} -- 2.11.0