From: Paul Brossier Date: Tue, 7 Aug 2018 12:05:49 +0000 (+0200) Subject: Merge branch 'master' into dct X-Git-Tag: 0.4.7~20^2~20 X-Git-Url: https://git.aubio.org/?p=aubio.git;a=commitdiff_plain;h=189c6ae6ed15ad2875bb0331110f681949fe2868;hp=fb4ab893f8d305b2aa969c766222a4c22f717294 Merge branch 'master' into dct --- diff --git a/python/lib/gen_code.py b/python/lib/gen_code.py index d3911974..b29a59fd 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 = { @@ -190,6 +192,11 @@ class MappedObject(object): out += self.gen_init() out += self.gen_del() out += self.gen_do() + if len(self.prototypes['rdo']): + self.do_proto = self.prototypes['rdo'][0] + self.do_inputs = [get_params_types_names(self.do_proto)[1]] + self.do_outputs = get_params_types_names(self.do_proto)[2:] + out += self.gen_do(method='rdo') out += self.gen_memberdef() out += self.gen_set() out += self.gen_get() @@ -370,12 +377,12 @@ Py_{shortname}_del (Py_{shortname} * self, PyObject * unused) """.format(del_fn = del_fn) return out - def gen_do(self): + def gen_do(self, method = 'do'): out = """ // do {shortname} static PyObject* -Py_{shortname}_do (Py_{shortname} * self, PyObject * args) -{{""".format(**self.__dict__) +Pyaubio_{shortname}_{method} (Py_{shortname} * self, PyObject * args) +{{""".format(**self.__dict__, method = method) input_params = self.do_inputs output_params = self.do_outputs #print input_params @@ -515,6 +522,12 @@ static PyMethodDef Py_{shortname}_methods[] = {{""".format(**self.__dict__) out += """ {{"{shortname}", (PyCFunction) Py{name}, METH_NOARGS, ""}},""".format(name = name, shortname = shortname) + for m in self.prototypes['rdo']: + name = get_name(m) + shortname = name.replace('aubio_%s_' % self.shortname, '') + out += """ + {{"{shortname}", (PyCFunction) Py{name}, + METH_VARARGS, ""}},""".format(name = name, shortname = shortname) out += """ {NULL} /* sentinel */ }; @@ -540,7 +553,7 @@ PyTypeObject Py_{shortname}Type = {{ 0, 0, 0, - (ternaryfunc)Py_{shortname}_do, + (ternaryfunc)Pyaubio_{shortname}_do, 0, 0, 0, diff --git a/python/lib/gen_external.py b/python/lib/gen_external.py index 8095cb60..1e35bfd0 100644 --- a/python/lib/gen_external.py +++ b/python/lib/gen_external.py @@ -181,7 +181,7 @@ def generate_lib_from_c_declarations(cpp_objects, c_declarations): if o[:6] == 'aubio_': shortname = o[6:-2] # without aubio_ prefix and _t suffix - lib[shortname] = {'struct': [], 'new': [], 'del': [], 'do': [], 'get': [], 'set': [], 'other': []} + lib[shortname] = {'struct': [], 'new': [], 'del': [], 'do': [], 'rdo': [], 'get': [], 'set': [], 'other': []} lib[shortname]['longname'] = o lib[shortname]['shortname'] = shortname @@ -195,6 +195,8 @@ def generate_lib_from_c_declarations(cpp_objects, c_declarations): lib[shortname]['struct'].append(fn) elif '_do' in fn: lib[shortname]['do'].append(fn) + elif '_rdo' in fn: + lib[shortname]['rdo'].append(fn) elif 'new_' in fn: lib[shortname]['new'].append(fn) elif 'del_' in fn: diff --git a/python/tests/test_dct.py b/python/tests/test_dct.py new file mode 100755 index 00000000..07c336d9 --- /dev/null +++ b/python/tests/test_dct.py @@ -0,0 +1,68 @@ +#! /usr/bin/env python + + +import numpy as np +from numpy.testing import TestCase, assert_almost_equal +import aubio + +precomputed_arange = [ 9.89949512, -6.44232273, 0., -0.67345482, 0., + -0.20090288, 0., -0.05070186] + +precomputed_some_ones = [ 4.28539848, 0.2469689, -0.14625292, -0.58121818, + -0.83483052, -0.75921834, -0.35168475, 0.24087936, + 0.78539824, 1.06532764, 0.97632152, 0.57164496, 0.03688532, + -0.39446154, -0.54619485, -0.37771079] + +class aubio_dct(TestCase): + + def test_init(self): + """ test that aubio.dct() is created with expected size """ + a_dct = aubio.dct() + self.assertEqual(a_dct.size, 1024) + + def test_arange(self): + """ test that dct(arange(8)) is computed correctly + + >>> from scipy.fftpack import dct + >>> a_in = np.arange(8).astype('float32') + >>> precomputed = dct(a_in, norm='ortho') + """ + N = len(precomputed_arange) + a_dct = aubio.dct(8) + a_in = np.arange(8).astype('float32') + a_expected = aubio.fvec(precomputed_arange) + assert_almost_equal(a_dct(a_in), a_expected, decimal=6) + + def test_some_ones(self): + """ test that dct(somevector) is computed correctly """ + a_dct = aubio.dct(16) + a_in = np.ones(16).astype('float32') + a_in[1] = 0 + a_in[3] = np.pi + a_expected = aubio.fvec(precomputed_some_ones) + assert_almost_equal(a_dct(a_in), a_expected, decimal=6) + + def test_reconstruction(self): + """ test that some_ones vector can be recontructed """ + a_dct = aubio.dct(16) + a_in = np.ones(16).astype('float32') + a_in[1] = 0 + a_in[3] = np.pi + a_dct_in = a_dct(a_in) + a_dct_reconstructed = a_dct.rdo(a_dct_in) + assert_almost_equal(a_dct_reconstructed, a_in, decimal=6) + + def test_negative_size(self): + """ test that creation fails with a negative size """ + with self.assertRaises(ValueError): + aubio.dct(-1) + + def test_wrong_size(self): + """ test that creation fails with a non power-of-two size """ + # supports for non 2** fft sizes only when compiled with fftw3 + size = 13 + try: + with self.assertRaises(RuntimeError): + aubio.dct(size) + except AssertionError: + self.skipTest('creating aubio.dct with size %d did not fail' % size) 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/aubio_priv.h b/src/aubio_priv.h index 65a3d8a3..1e9881e1 100644 --- a/src/aubio_priv.h +++ b/src/aubio_priv.h @@ -85,6 +85,8 @@ #ifndef HAVE_AUBIO_DOUBLE #define aubio_vDSP_mmov vDSP_mmov #define aubio_vDSP_vmul vDSP_vmul +#define aubio_vDSP_vsmul vDSP_vsmul +#define aubio_vDSP_vsadd vDSP_vsadd #define aubio_vDSP_vfill vDSP_vfill #define aubio_vDSP_meanv vDSP_meanv #define aubio_vDSP_sve vDSP_sve @@ -97,6 +99,8 @@ #else /* HAVE_AUBIO_DOUBLE */ #define aubio_vDSP_mmov vDSP_mmovD #define aubio_vDSP_vmul vDSP_vmulD +#define aubio_vDSP_vsmul vDSP_vsmulD +#define aubio_vDSP_vsadd vDSP_vsaddD #define aubio_vDSP_vfill vDSP_vfillD #define aubio_vDSP_meanv vDSP_meanvD #define aubio_vDSP_sve vDSP_sveD 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_accelerate.c b/src/spectral/dct_accelerate.c new file mode 100644 index 00000000..8bf5f8e2 --- /dev/null +++ b/src/spectral/dct_accelerate.c @@ -0,0 +1,99 @@ +/* + 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" + +#if defined(HAVE_ACCELERATE) + +#if HAVE_AUBIO_DOUBLE +#warning "no double-precision dct with accelerate" +#endif + +struct _aubio_dct_t { + uint_t size; + fvec_t *tmp; + vDSP_DFT_Setup setup; + vDSP_DFT_Setup setupInv; +}; + +aubio_dct_t * new_aubio_dct (uint_t size) { + aubio_dct_t * s = AUBIO_NEW(aubio_dct_t); + + if ((sint_t)size < 16 || !aubio_is_power_of_two(size)) { + AUBIO_ERR("dct: can only create with sizes greater than 16 and" + "that are powers of two, requested %d\n", size); + goto beach; + } + + s->setup = vDSP_DCT_CreateSetup(NULL, (vDSP_Length)size, vDSP_DCT_II); + s->setupInv = vDSP_DCT_CreateSetup(NULL, (vDSP_Length)size, vDSP_DCT_III); + if (s->setup == NULL || s->setupInv == NULL) { + goto beach; + } + + s->size = size; + + return s; + +beach: + del_aubio_dct(s); + return NULL; +} + +void del_aubio_dct(aubio_dct_t *s) { + if (s->setup) vDSP_DFT_DestroySetup(s->setup); + if (s->setupInv) vDSP_DFT_DestroySetup(s->setupInv); + AUBIO_FREE(s); +} + +void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + + vDSP_DCT_Execute(s->setup, (const float *)input->data, (float *)output->data); + + // apply orthonormal scaling + output->data[0] *= SQRT(1./s->size); + smpl_t scaler = SQRT(2./s->size); + + aubio_vDSP_vsmul(output->data + 1, 1, &scaler, output->data + 1, 1, + output->length - 1); + +} + +void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + + output->data[0] = input->data[0] / SQRT(1./s->size); + smpl_t scaler = 1./SQRT(2./s->size); + + aubio_vDSP_vsmul(input->data + 1, 1, &scaler, output->data + 1, 1, + output->length - 1); + + vDSP_DCT_Execute(s->setupInv, (const float *)output->data, + (float *)output->data); + + scaler = 2./s->size; + + aubio_vDSP_vsmul(output->data, 1, &scaler, output->data, 1, output->length); + +} + +#endif //defined(HAVE_ACCELERATE) diff --git a/src/spectral/dct_fftw.c b/src/spectral/dct_fftw.c new file mode 100644 index 00000000..7ea1f2a1 --- /dev/null +++ b/src/spectral/dct_fftw.c @@ -0,0 +1,129 @@ +/* + 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" + +#ifdef HAVE_FFTW3 + +#include +#include + +#ifdef HAVE_FFTW3F +#if HAVE_AUBIO_DOUBLE +#error "Using aubio in double precision with fftw3 in single precision" +#endif /* HAVE_AUBIO_DOUBLE */ +#else /* HAVE_FFTW3F */ +#if !HAVE_AUBIO_DOUBLE +#error "Using aubio in single precision with fftw3 in double precision" +#endif /* HAVE_AUBIO_DOUBLE */ +#endif /* HAVE_FFTW3F */ + +#ifdef HAVE_FFTW3F +#define fftw_malloc fftwf_malloc +#define fftw_free fftwf_free +#define fftw_execute fftwf_execute +#define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d +#define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d +#define fftw_plan_r2r_1d fftwf_plan_r2r_1d +#define fftw_plan fftwf_plan +#define fftw_destroy_plan fftwf_destroy_plan +#endif + +// defined in src/spectral/fft.c +extern pthread_mutex_t aubio_fftw_mutex; + +extern void aubio_ooura_ddct(int, int, smpl_t *, int *, smpl_t *); + +struct _aubio_dct_t { + uint_t size; + fvec_t *in, *out; + smpl_t *data; + fftw_plan pfw, pbw; + smpl_t scalers[5]; +}; + +aubio_dct_t * new_aubio_dct (uint_t size) { + aubio_dct_t * s = AUBIO_NEW(aubio_dct_t); + if (!s) { + goto beach; + } + s->size = size; + s->in = new_fvec(size); + s->out = new_fvec(size); + pthread_mutex_lock(&aubio_fftw_mutex); + s->data = (smpl_t *)fftw_malloc(sizeof(smpl_t) * size); + s->pfw = fftw_plan_r2r_1d(size, s->in->data, s->data, FFTW_REDFT10, + FFTW_ESTIMATE); + s->pbw = fftw_plan_r2r_1d(size, s->data, s->out->data, FFTW_REDFT01, + FFTW_ESTIMATE); + pthread_mutex_unlock(&aubio_fftw_mutex); + s->scalers[0] = SQRT(1./(4.*s->size)); + s->scalers[1] = SQRT(1./(2.*s->size)); + s->scalers[2] = 1. / s->scalers[0]; + s->scalers[3] = 1. / s->scalers[1]; + s->scalers[4] = .5 / s->size; + return s; +beach: + AUBIO_FREE(s); + return NULL; +} + +void del_aubio_dct(aubio_dct_t *s) { + pthread_mutex_lock(&aubio_fftw_mutex); + fftw_destroy_plan(s->pfw); + fftw_destroy_plan(s->pbw); + fftw_free(s->data); + pthread_mutex_unlock(&aubio_fftw_mutex); + del_fvec(s->in); + del_fvec(s->out); + AUBIO_FREE(s); +} + +void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + uint_t i; + fvec_copy(input, s->in); + fftw_execute(s->pfw); + //fvec_copy(s->out, output); + s->data[0] *= s->scalers[0]; + for (i = 1; i < s->size; i++) { + s->data[i] *= s->scalers[1]; + } + memcpy(output->data, s->data, output->length * sizeof(smpl_t)); +} + +void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + uint_t i; + memcpy(s->data, input->data, input->length * sizeof(smpl_t)); + //s->data[0] *= .5; + s->data[0] *= s->scalers[2]; + for (i = 1; i < s->size; i++) { + s->data[i] *= s->scalers[3]; + } + fftw_execute(s->pbw); + for (i = 0; i < s->size; i++) { + s->out->data[i] *= s->scalers[4]; + } + fvec_copy(s->out, output); +} + +#endif //HAVE_FFTW3 diff --git a/src/spectral/dct_ipp.c b/src/spectral/dct_ipp.c new file mode 100644 index 00000000..cfaba388 --- /dev/null +++ b/src/spectral/dct_ipp.c @@ -0,0 +1,144 @@ +/* + 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" + +#if defined(HAVE_INTEL_IPP) + +#if !HAVE_AUBIO_DOUBLE +#define aubio_IppFloat Ipp32f +#define aubio_ippsDCTFwdSpec IppsDCTFwdSpec_32f +#define aubio_ippsDCTInvSpec IppsDCTInvSpec_32f +#define aubio_ippsDCTFwdGetSize ippsDCTFwdGetSize_32f +#define aubio_ippsDCTInvGetSize ippsDCTInvGetSize_32f +#define aubio_ippsDCTFwdInit ippsDCTFwdInit_32f +#define aubio_ippsDCTInvInit ippsDCTInvInit_32f +#define aubio_ippsDCTFwd ippsDCTFwd_32f +#define aubio_ippsDCTInv ippsDCTInv_32f +#else /* HAVE_AUBIO_DOUBLE */ +#define aubio_IppFloat Ipp64f +#define aubio_ippsDCTFwdSpec IppsDCTFwdSpec_64f +#define aubio_ippsDCTInvSpec IppsDCTInvSpec_64f +#define aubio_ippsDCTFwdGetSize ippsDCTFwdGetSize_64f +#define aubio_ippsDCTInvGetSize ippsDCTInvGetSize_64f +#define aubio_ippsDCTFwdInit ippsDCTFwdInit_64f +#define aubio_ippsDCTInvInit ippsDCTInvInit_64f +#define aubio_ippsDCTFwd ippsDCTFwd_64f +#define aubio_ippsDCTInv ippsDCTInv_64f +#endif + +struct _aubio_dct_t { + uint_t size; + Ipp8u* pSpec; + Ipp8u* pSpecBuffer; + Ipp8u* pBuffer; + aubio_ippsDCTFwdSpec* pFwdDCTSpec; + aubio_ippsDCTInvSpec* pInvDCTSpec; +}; + +aubio_dct_t * new_aubio_dct (uint_t size) { + aubio_dct_t * s = AUBIO_NEW(aubio_dct_t); + + const IppHintAlgorithm qualityHint = ippAlgHintAccurate; // ippAlgHintFast; + int pSpecSize, pSpecBufferSize, pBufferSize; + IppStatus status; + + if ((sint_t)size <= 1) { + AUBIO_ERR("dct: can only create with sizes greater than 1, requested %d\n", + size); + goto beach; + } + + status = aubio_ippsDCTFwdGetSize(size, qualityHint, &pSpecSize, + &pSpecBufferSize, &pBufferSize); + if (status != ippStsNoErr) { + AUBIO_ERR("dct: failed to initialize dct. IPP error: %d\n", status); + goto beach; + } + + //AUBIO_INF("dct: fwd initialized with %d %d %d\n", pSpecSize, pSpecBufferSize, + // pBufferSize); + + s->pSpec = ippsMalloc_8u(pSpecSize); + if (pSpecSize > 0) { + s->pSpecBuffer = ippsMalloc_8u(pSpecBufferSize); + } else { + s->pSpecBuffer = NULL; + } + s->pBuffer = ippsMalloc_8u(pBufferSize); + + status = aubio_ippsDCTInvGetSize(size, qualityHint, &pSpecSize, + &pSpecBufferSize, &pBufferSize); + if (status != ippStsNoErr) { + AUBIO_ERR("dct: failed to initialize dct. IPP error: %d\n", status); + goto beach; + } + + //AUBIO_INF("dct: inv initialized with %d %d %d\n", pSpecSize, pSpecBufferSize, + // pBufferSize); + + status = aubio_ippsDCTFwdInit(&(s->pFwdDCTSpec), size, qualityHint, s->pSpec, + s->pSpecBuffer); + if (status != ippStsNoErr) { + AUBIO_ERR("dct: failed to initialize fwd dct. IPP error: %d\n", status); + goto beach; + } + + status = aubio_ippsDCTInvInit(&(s->pInvDCTSpec), size, qualityHint, s->pSpec, + s->pSpecBuffer); + if (status != ippStsNoErr) { + AUBIO_ERR("dct: failed to initialize inv dct. IPP error: %d\n", status); + goto beach; + } + + s->size = size; + + return s; + +beach: + del_aubio_dct(s); + return NULL; +} + +void del_aubio_dct(aubio_dct_t *s) { + ippFree(s->pSpec); + ippFree(s->pSpecBuffer); + ippFree(s->pBuffer); + AUBIO_FREE(s); +} + +void aubio_dct_do(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + + aubio_ippsDCTFwd((const aubio_IppFloat*)input->data, + (aubio_IppFloat*)output->data, s->pFwdDCTSpec, s->pBuffer); + +} + +void aubio_dct_rdo(aubio_dct_t *s, const fvec_t *input, fvec_t *output) { + + aubio_ippsDCTInv((const aubio_IppFloat*)input->data, + (aubio_IppFloat*)output->data, s->pInvDCTSpec, s->pBuffer); + +} + +#endif //defined(HAVE_INTEL_IPP) diff --git a/src/spectral/dct_ooura.c b/src/spectral/dct_ooura.c new file mode 100644 index 00000000..38f1949e --- /dev/null +++ b/src/spectral/dct_ooura.c @@ -0,0 +1,95 @@ +/* + 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" + +#if !defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3) && !defined(HAVE_INTEL_IPP) + +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; + smpl_t scalers[5]; +}; + +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; + s->scalers[0] = 2. * SQRT(1./(4.*s->size)); + s->scalers[1] = 2. * SQRT(1./(2.*s->size)); + s->scalers[2] = 1. / s->scalers[0]; + s->scalers[3] = 1. / s->scalers[1]; + s->scalers[4] = 2. / s->size; + 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] *= s->scalers[0]; + for (i = 1; i < s->input->length; i++) { + s->input->data[i] *= s->scalers[1]; + } + 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] *= s->scalers[2]; + for (i = 1; i < s->input->length; i++) { + s->input->data[i] *= s->scalers[3]; + } + 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] *= s->scalers[4]; + } + fvec_copy(s->input, output); +} + +#endif //!defined(HAVE_ACCELERATE) && !defined(HAVE_FFTW3) diff --git a/tests/src/spectral/test-dct.c b/tests/src/spectral/test-dct.c new file mode 100644 index 00000000..5ea4d1dd --- /dev/null +++ b/tests/src/spectral/test-dct.c @@ -0,0 +1,31 @@ +#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++) { + 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; +}