Merge branch 'master' into feature/pytest
authorPaul Brossier <piem@piem.org>
Sat, 17 Nov 2018 15:38:29 +0000 (16:38 +0100)
committerPaul Brossier <piem@piem.org>
Sat, 17 Nov 2018 15:38:29 +0000 (16:38 +0100)
24 files changed:
Makefile
doc/py_utils.rst
python/ext/aubiomodule.c
python/ext/py-filterbank.c
python/ext/py-musicutils.c
python/ext/py-musicutils.h
python/lib/gen_code.py
python/tests/test_fft.py
python/tests/test_filterbank_mel.py
python/tests/test_hztomel.py [new file with mode: 0755]
python/tests/test_source_channels.py [new file with mode: 0755]
src/io/source_avcodec.c
src/mathutils.c
src/mathutils.h
src/musicutils.c [new file with mode: 0644]
src/musicutils.h
src/spectral/fft.c
src/spectral/filterbank.c
src/spectral/filterbank.h
src/spectral/filterbank_mel.c
src/spectral/filterbank_mel.h
src/spectral/mfcc.c
src/spectral/mfcc.h
tests/src/spectral/test-filterbank.c

index 92257cb..380d313 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -232,6 +232,11 @@ test_python_only_clean: test_python_only \
        uninstall_python \
        check_clean_python
 
+coverage_cycle: coverage_zero_counters coverage_report
+
+coverage_zero_counters:
+       lcov --zerocounters --directory .
+
 coverage: export CFLAGS=--coverage
 coverage: export LDFLAGS=--coverage
 coverage: export PYTHONPATH=$(PWD)/python/lib
@@ -244,6 +249,8 @@ coverage: force_uninstall_python deps_python \
        lcov --capture --no-external --directory . --output-file build/coverage_python.info
        lcov -a build/coverage_python.info -a build/coverage_lib.info -o build/coverage.info
 
+# make sure we don't build the doc, which builds a temporary python module
+coverage_report: export WAFOPTS += --disable-docs
 coverage_report: coverage
        genhtml build/coverage.info --output-directory lcov_html
        mkdir -p gcovr_html/
index 6d19c0b..94debf8 100644 (file)
@@ -28,6 +28,12 @@ Frequency conversion
 
 .. autofunction:: miditofreq
 
+.. python/ext/py-musicutils.h
+
+.. autofunction:: meltohz
+
+.. autofunction:: hztomel
+
 .. python/ext/aubiomodule.c
 
 .. autofunction:: bintomidi
index 0fa8ac7..0345064 100644 (file)
@@ -372,6 +372,10 @@ static PyMethodDef aubio_methods[] = {
   {"window", Py_aubio_window, METH_VARARGS, Py_aubio_window_doc},
   {"shift", Py_aubio_shift, METH_VARARGS, Py_aubio_shift_doc},
   {"ishift", Py_aubio_ishift, METH_VARARGS, Py_aubio_ishift_doc},
+  {"hztomel", Py_aubio_hztomel, METH_VARARGS|METH_KEYWORDS, Py_aubio_hztomel_doc},
+  {"meltohz", Py_aubio_meltohz, METH_VARARGS|METH_KEYWORDS, Py_aubio_meltohz_doc},
+  {"hztomel_htk", Py_aubio_hztomel_htk, METH_VARARGS, Py_aubio_hztomel_htk_doc},
+  {"meltohz_htk", Py_aubio_meltohz_htk, METH_VARARGS, Py_aubio_meltohz_htk_doc},
   {NULL, NULL, 0, NULL} /* Sentinel */
 };
 
index a6647d5..24f9f71 100644 (file)
@@ -138,8 +138,14 @@ Py_filterbank_set_triangle_bands (Py_filterbank * self, PyObject *args)
   err = aubio_filterbank_set_triangle_bands (self->o,
       &(self->freqs), samplerate);
   if (err > 0) {
-    PyErr_SetString (PyExc_ValueError,
-        "error when running set_triangle_bands");
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError, "error running set_triangle_bands");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
     return NULL;
   }
   Py_RETURN_NONE;
@@ -150,15 +156,79 @@ Py_filterbank_set_mel_coeffs_slaney (Py_filterbank * self, PyObject *args)
 {
   uint_t err = 0;
 
-  uint_t samplerate;
-  if (!PyArg_ParseTuple (args, "I", &samplerate)) {
+  smpl_t samplerate;
+  if (!PyArg_ParseTuple (args, AUBIO_NPY_SMPL_CHR, &samplerate)) {
     return NULL;
   }
 
   err = aubio_filterbank_set_mel_coeffs_slaney (self->o, samplerate);
   if (err > 0) {
-    PyErr_SetString (PyExc_ValueError,
-        "error when running set_mel_coeffs_slaney");
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError, "error running set_mel_coeffs_slaney");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
+    return NULL;
+  }
+  Py_RETURN_NONE;
+}
+
+static PyObject *
+Py_filterbank_set_mel_coeffs (Py_filterbank * self, PyObject *args)
+{
+  uint_t err = 0;
+
+  smpl_t samplerate;
+  smpl_t freq_min;
+  smpl_t freq_max;
+  if (!PyArg_ParseTuple (args, AUBIO_NPY_SMPL_CHR AUBIO_NPY_SMPL_CHR
+        AUBIO_NPY_SMPL_CHR, &samplerate, &freq_min, &freq_max)) {
+    return NULL;
+  }
+
+  err = aubio_filterbank_set_mel_coeffs (self->o, samplerate,
+      freq_min, freq_max);
+  if (err > 0) {
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError, "error running set_mel_coeffs");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
+    return NULL;
+  }
+  Py_RETURN_NONE;
+}
+
+static PyObject *
+Py_filterbank_set_mel_coeffs_htk (Py_filterbank * self, PyObject *args)
+{
+  uint_t err = 0;
+
+  smpl_t samplerate;
+  smpl_t freq_min;
+  smpl_t freq_max;
+  if (!PyArg_ParseTuple (args, AUBIO_NPY_SMPL_CHR AUBIO_NPY_SMPL_CHR
+        AUBIO_NPY_SMPL_CHR, &samplerate, &freq_min, &freq_max)) {
+    return NULL;
+  }
+
+  err = aubio_filterbank_set_mel_coeffs_htk (self->o, samplerate,
+      freq_min, freq_max);
+  if (err > 0) {
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError, "error running set_mel_coeffs_htk");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
     return NULL;
   }
   Py_RETURN_NONE;
@@ -195,15 +265,69 @@ Py_filterbank_get_coeffs (Py_filterbank * self, PyObject *unused)
       aubio_filterbank_get_coeffs (self->o) );
 }
 
+static PyObject *
+Py_filterbank_set_power(Py_filterbank *self, PyObject *args)
+{
+  uint_t power;
+
+  if (!PyArg_ParseTuple (args, "I", &power)) {
+    return NULL;
+  }
+  if(aubio_filterbank_set_power (self->o, power)) {
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError,
+          "error running filterbank.set_power");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
+    return NULL;
+  }
+  Py_RETURN_NONE;
+}
+
+static PyObject *
+Py_filterbank_set_norm(Py_filterbank *self, PyObject *args)
+{
+  uint_t playing;
+
+  if (!PyArg_ParseTuple (args, "I", &playing)) {
+    return NULL;
+  }
+  if(aubio_filterbank_set_norm (self->o, playing)) {
+    if (PyErr_Occurred() == NULL) {
+      PyErr_SetString (PyExc_ValueError,
+          "error running filterbank.set_power");
+    } else {
+      // change the RuntimeError into ValueError
+      PyObject *type, *value, *traceback;
+      PyErr_Fetch(&type, &value, &traceback);
+      PyErr_Restore(PyExc_ValueError, value, traceback);
+    }
+    return NULL;
+  }
+  Py_RETURN_NONE;
+}
+
 static PyMethodDef Py_filterbank_methods[] = {
   {"set_triangle_bands", (PyCFunction) Py_filterbank_set_triangle_bands,
     METH_VARARGS, "set coefficients of filterbanks"},
   {"set_mel_coeffs_slaney", (PyCFunction) Py_filterbank_set_mel_coeffs_slaney,
     METH_VARARGS, "set coefficients of filterbank as in Auditory Toolbox"},
+  {"set_mel_coeffs", (PyCFunction) Py_filterbank_set_mel_coeffs,
+    METH_VARARGS, "set coefficients of filterbank to linearly spaced mel scale"},
+  {"set_mel_coeffs_htk", (PyCFunction) Py_filterbank_set_mel_coeffs_htk,
+    METH_VARARGS, "set coefficients of filterbank to linearly spaced mel scale"},
   {"get_coeffs", (PyCFunction) Py_filterbank_get_coeffs,
     METH_NOARGS, "get coefficients of filterbank"},
   {"set_coeffs", (PyCFunction) Py_filterbank_set_coeffs,
     METH_VARARGS, "set coefficients of filterbank"},
+  {"set_power", (PyCFunction) Py_filterbank_set_power,
+    METH_VARARGS, "set power applied to filterbank input spectrum"},
+  {"set_norm", (PyCFunction) Py_filterbank_set_norm,
+    METH_VARARGS, "set norm applied to filterbank input spectrum"},
   {NULL}
 };
 
index a0d450a..3223945 100644 (file)
@@ -181,3 +181,57 @@ Py_aubio_ishift(PyObject *self, PyObject *args)
   //Py_RETURN_NONE;
   return (PyObject *) PyAubio_CFvecToArray(&vec);
 }
+
+PyObject*
+Py_aubio_hztomel(PyObject *self, PyObject *args, PyObject *kwds)
+{
+  smpl_t v;
+  PyObject *htk = NULL;
+  static char *kwlist[] = {"f", "htk", NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, AUBIO_NPY_SMPL_CHR "|O",
+        kwlist, &v, &htk))
+  {
+    return NULL;
+  }
+  if (htk != NULL && PyObject_IsTrue(htk) == 1)
+    return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_hztomel_htk(v));
+  else
+    return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_hztomel(v));
+}
+
+PyObject*
+Py_aubio_meltohz(PyObject *self, PyObject *args, PyObject *kwds)
+{
+  smpl_t v;
+  PyObject *htk = NULL;
+  static char *kwlist[] = {"m", "htk", NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, AUBIO_NPY_SMPL_CHR "|O",
+        kwlist, &v, &htk))
+  {
+    return NULL;
+  }
+  if (htk != NULL && PyObject_IsTrue(htk) == 1)
+    return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_meltohz_htk(v));
+  else
+    return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_meltohz(v));
+}
+
+PyObject*
+Py_aubio_hztomel_htk(PyObject *self, PyObject *args)
+{
+  smpl_t v;
+  if (!PyArg_ParseTuple(args, AUBIO_NPY_SMPL_CHR, &v)) {
+    return NULL;
+  }
+  return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_hztomel_htk(v));
+}
+
+PyObject*
+Py_aubio_meltohz_htk(PyObject *self, PyObject *args)
+{
+  smpl_t v;
+  if (!PyArg_ParseTuple(args, AUBIO_NPY_SMPL_CHR, &v)) {
+    return NULL;
+  }
+  return Py_BuildValue(AUBIO_NPY_SMPL_CHR, aubio_meltohz_htk(v));
+}
index 02415e4..27698d1 100644 (file)
@@ -300,4 +300,72 @@ static char Py_aubio_ishift_doc[] = ""
 "";
 PyObject * Py_aubio_ishift(PyObject *self, PyObject *args);
 
+static char Py_aubio_hztomel_doc[] = ""
+"hztomel(f, htk=False)\n"
+"\n"
+"Convert a scalar from frequency to mel scale.\n"
+"\n"
+"Parameters\n"
+"----------\n"
+"m : float\n"
+"   input frequency, in Hz\n"
+"htk : bool\n"
+"   if `True`, use Htk mel scale instead of Slaney.\n"
+"\n"
+"Returns\n"
+"-------\n"
+"float\n"
+"   output mel\n"
+"\n"
+"See Also\n"
+"--------\n"
+"meltohz\n"
+"";
+PyObject * Py_aubio_hztomel(PyObject *self, PyObject *args);
+
+static char Py_aubio_meltohz_doc[] = ""
+"meltohz(m, htk=False)\n"
+"\n"
+"Convert a scalar from mel scale to frequency.\n"
+"\n"
+"Parameters\n"
+"----------\n"
+"m : float\n"
+"   input mel\n"
+"htk : bool\n"
+"   if `True`, use Htk mel scale instead of Slaney.\n"
+"\n"
+"Returns\n"
+"-------\n"
+"float\n"
+"   output frequency, in Hz\n"
+"\n"
+"See Also\n"
+"--------\n"
+"hztomel\n"
+"";
+PyObject * Py_aubio_meltohz(PyObject *self, PyObject *args);
+
+static char Py_aubio_hztomel_htk_doc[] = ""
+"hztomel_htk(m)\n"
+"\n"
+"Same as `hztomel(m, htk=True)`\n"
+"\n"
+"See Also\n"
+"--------\n"
+"hztomel\n"
+"";
+PyObject * Py_aubio_hztomel_htk(PyObject *self, PyObject *args);
+
+static char Py_aubio_meltohz_htk_doc[] = ""
+"meltohz_htk(m)\n"
+"\n"
+"Same as `meltohz(m, htk=True)`\n"
+"\n"
+"See Also\n"
+"--------\n"
+"meltohz\n"
+"";
+PyObject * Py_aubio_meltohz_htk(PyObject *self, PyObject *args);
+
 #endif /* PY_AUBIO_MUSICUTILS_H */
index 1266336..a29a23f 100644 (file)
@@ -462,22 +462,35 @@ Pyaubio_{shortname}_{method}  (Py_{shortname} * self, PyObject * args)
 // {shortname} setters
 """.format(**self.__dict__)
         for set_param in self.prototypes['set']:
-            params = get_params_types_names(set_param)[1]
-            paramtype = params['type']
+            params = get_params_types_names(set_param)[1:]
+            param = self.shortname.split('_set_')[-1]
+            paramdecls = "".join(["""
+   {0} {1};""".format(p['type'], p['name']) for p in params])
             method_name = get_name(set_param)
             param = method_name.split('aubio_'+self.shortname+'_set_')[-1]
-            pyparamtype = pyargparse_chars[paramtype]
+            refs = ", ".join(["&%s" % p['name'] for p in params])
+            paramlist = ", ".join(["%s" % p['name'] for p in params])
+            if len(params):
+                paramlist = "," + paramlist
+            pyparamtypes = ''.join([pyargparse_chars[p['type']] for p in params])
             out += """
 static PyObject *
 Pyaubio_{shortname}_set_{param} (Py_{shortname} *self, PyObject *args)
 {{
   uint_t err = 0;
-  {paramtype} {param};
+  {paramdecls}
+""".format(param = param, paramdecls = paramdecls, **self.__dict__)
 
-  if (!PyArg_ParseTuple (args, "{pyparamtype}", &{param})) {{
+            if len(refs) and len(pyparamtypes):
+                out += """
+
+  if (!PyArg_ParseTuple (args, "{pyparamtypes}", {refs})) {{
     return NULL;
   }}
-  err = aubio_{shortname}_set_{param} (self->o, {param});
+""".format(pyparamtypes = pyparamtypes, refs = refs)
+
+            out += """
+  err = aubio_{shortname}_set_{param} (self->o {paramlist});
 
   if (err > 0) {{
     if (PyErr_Occurred() == NULL) {{
@@ -492,7 +505,8 @@ Pyaubio_{shortname}_set_{param} (Py_{shortname} *self, PyObject *args)
   }}
   Py_RETURN_NONE;
 }}
-""".format(param = param, paramtype = paramtype, pyparamtype = pyparamtype, **self.__dict__)
+""".format(param = param, refs = refs, paramdecls = paramdecls,
+        pyparamtypes = pyparamtypes, paramlist = paramlist, **self.__dict__)
         return out
 
     def gen_get(self):
index 5b34411..abe95d3 100755 (executable)
@@ -141,6 +141,37 @@ class aubio_fft_test_case(TestCase):
         assert_almost_equal ( r[0], impulse, decimal = 6)
         assert_almost_equal ( r[1:], 0)
 
+class aubio_fft_odd_sizes(TestCase):
+
+    def test_reconstruct_with_odd_size(self):
+        win_s = 29
+        self.recontruct(win_s, 'odd sizes not supported')
+
+    def test_reconstruct_with_radix15(self):
+        win_s = 2 ** 4 * 15
+        self.recontruct(win_s, 'radix 15 supported')
+
+    def test_reconstruct_with_radix5(self):
+        win_s = 2 ** 4 * 5
+        self.recontruct(win_s, 'radix 5 supported')
+
+    def test_reconstruct_with_radix3(self):
+        win_s = 2 ** 4 * 3
+        self.recontruct(win_s, 'radix 3 supported')
+
+    def recontruct(self, win_s, skipMessage):
+        try:
+            f = fft(win_s)
+        except RuntimeError:
+            self.skipTest(skipMessage)
+        input_signal = fvec(win_s)
+        input_signal[win_s//2] = 1
+        c = f(input_signal)
+        output_signal = f.rdo(c)
+        assert_almost_equal(input_signal, output_signal)
+
+class aubio_fft_wrong_params(TestCase):
+
     def test_large_input_timegrain(self):
         win_s = 1024
         f = fft(win_s)
@@ -169,22 +200,11 @@ class aubio_fft_test_case(TestCase):
         with self.assertRaises(ValueError):
             f.rdo(s)
 
-class aubio_fft_wrong_params(TestCase):
-
     def test_wrong_buf_size(self):
         win_s = -1
         with self.assertRaises(ValueError):
             fft(win_s)
 
-    def test_buf_size_not_power_of_two(self):
-        # when compiled with fftw3, aubio supports non power of two fft sizes
-        win_s = 320
-        try:
-            with self.assertRaises(RuntimeError):
-                fft(win_s)
-        except AssertionError:
-            self.skipTest('creating aubio.fft with size %d did not fail' % win_s)
-
     def test_buf_size_too_small(self):
         win_s = 1
         with self.assertRaises(RuntimeError):
index 90161fa..0771264 100755 (executable)
@@ -109,6 +109,63 @@ class aubio_filterbank_mel_test_case(TestCase):
             f.set_triangle_bands(fvec(freq_list), samplerate)
 
 
+    def test_triangle_freqs_without_norm(self):
+        """make sure set_triangle_bands works without """
+        samplerate = 22050
+        freq_list = fvec([0, 100, 1000, 10000])
+        f = filterbank(len(freq_list) - 2, 1024)
+        f.set_norm(0)
+        f.set_triangle_bands(freq_list, samplerate)
+        expected = f.get_coeffs()
+        f.set_norm(1)
+        f.set_triangle_bands(fvec(freq_list), samplerate)
+        assert_almost_equal(f.get_coeffs().T,
+                expected.T * 2. / (freq_list[2:] - freq_list[:-2]))
+
+    def test_triangle_freqs_wrong_norm(self):
+        f = filterbank(10, 1024)
+        with self.assertRaises(ValueError):
+            f.set_norm(-1)
+
+    def test_triangle_freqs_with_power(self):
+        f = filterbank(9, 1024)
+        freqs = fvec([40, 80, 200, 400, 800, 1600, 3200, 6400, 12800, 15000,
+            24000])
+        f.set_power(2)
+        f.set_triangle_bands(freqs, 48000)
+        spec = cvec(1024)
+        spec.norm[:] = .1
+        expected = fvec([0.02070313, 0.02138672, 0.02127604, 0.02135417,
+            0.02133301, 0.02133301, 0.02133311, 0.02133334, 0.02133345])
+        expected /= 100.
+        assert_almost_equal(f(spec), expected)
+
+    def test_mel_coeffs(self):
+        f = filterbank(40, 1024)
+        f.set_mel_coeffs(44100, 0, 44100 / 2)
+
+    def test_zero_fmax(self):
+        f = filterbank(40, 1024)
+        f.set_mel_coeffs(44100, 0, 0)
+
+    def test_wrong_mel_coeffs(self):
+        f = filterbank(40, 1024)
+        with self.assertRaises(ValueError):
+            f.set_mel_coeffs_slaney(0)
+        with self.assertRaises(ValueError):
+            f.set_mel_coeffs(44100, 0, -44100 / 2)
+        with self.assertRaises(ValueError):
+            f.set_mel_coeffs(44100, -0.1, 44100 / 2)
+        with self.assertRaises(ValueError):
+            f.set_mel_coeffs(-44100, 0.1, 44100 / 2)
+        with self.assertRaises(ValueError):
+            f.set_mel_coeffs_htk(-1, 0, 0)
+
+    def test_mel_coeffs_htk(self):
+        f = filterbank(40, 1024)
+        f.set_mel_coeffs_htk(44100, 0, 44100 / 2)
+
+
 if __name__ == '__main__':
     from unittest import main
     main()
diff --git a/python/tests/test_hztomel.py b/python/tests/test_hztomel.py
new file mode 100755 (executable)
index 0000000..c7c402d
--- /dev/null
@@ -0,0 +1,109 @@
+#! /usr/bin/env python
+
+from unittest import main
+from numpy.testing import TestCase
+from numpy.testing import assert_equal, assert_almost_equal
+import numpy as np
+import aubio
+
+from aubio import hztomel, meltohz
+from aubio import hztomel_htk, meltohz_htk
+
+
+class aubio_hztomel_test_case(TestCase):
+
+    def test_hztomel(self):
+        assert_equal(hztomel(0.), 0.)
+        assert_almost_equal(hztomel(400. / 3.), 2., decimal=5)
+        assert_almost_equal(hztomel(1000. / 3), 5.)
+        assert_equal(hztomel(200.), 3.)
+        assert_almost_equal(hztomel(1000.), 15)
+        assert_almost_equal(hztomel(6400), 42)
+        assert_almost_equal(hztomel(40960), 69)
+
+        for m in np.linspace(0, 1000, 100):
+            assert_almost_equal(hztomel(meltohz(m)) - m, 0, decimal=3)
+
+    def test_meltohz(self):
+        assert_equal(meltohz(0.), 0.)
+        assert_almost_equal(meltohz(2), 400. / 3., decimal=4)
+        assert_equal(meltohz(3.), 200.)
+        assert_almost_equal(meltohz(5), 1000. / 3., decimal=4)
+        assert_almost_equal(meltohz(15), 1000., decimal=4)
+        assert_almost_equal(meltohz(42), 6400., decimal=2)
+        assert_almost_equal(meltohz(69), 40960., decimal=1)
+
+        for f in np.linspace(0, 20000, 1000):
+            assert_almost_equal(meltohz(hztomel(f)) - f, 0, decimal=1)
+
+    def test_meltohz_negative(self):
+        # TODO add assert_warns
+        assert_equal(meltohz(-1), 0)
+
+    def test_hztomel_negative(self):
+        # TODO add assert_warns
+        assert_equal(hztomel(-1), 0)
+
+
+class aubio_hztomel_htk_test_case(TestCase):
+
+    def test_meltohz(self):
+        assert_equal(meltohz(0, htk=True), 0)
+        assert_almost_equal(meltohz(2595, htk=True), 6300., decimal=1)
+
+    def test_hztomel(self):
+        assert_equal(hztomel(0, htk=True), 0)
+        assert_almost_equal(hztomel(3428.7, htk=True), 2000., decimal=1)
+        assert_almost_equal(hztomel(6300, htk=True), 2595., decimal=1)
+
+    def test_meltohz_negative(self):
+        # TODO add assert_warns
+        assert_equal(meltohz(-1, htk=True), 0)
+        assert_almost_equal(meltohz(2000, htk=True), 3428.7, decimal=1)
+        assert_almost_equal(meltohz(1000, htk=True), 1000., decimal=1)
+
+    def test_hztomel_negative(self):
+        # TODO add assert_warns
+        assert_equal(hztomel(-1, htk=True), 0)
+        assert_almost_equal(hztomel(1000, htk=True), 1000., decimal=1)
+
+    def test_hztomel_htk(self):
+        for f in np.linspace(0, 20000, 1000):
+            assert_almost_equal(meltohz_htk(hztomel_htk(f)) - f, 0, decimal=1)
+        for f in np.linspace(0, 20000, 1000):
+            assert_almost_equal(hztomel_htk(meltohz_htk(f)) - f, 0, decimal=1)
+
+
+class aubio_hztomel_wrong_values(TestCase):
+    """ more tests to cover all branches """
+
+    def test_hztomel_wrong_values(self):
+        with self.assertRaises(TypeError):
+            hztomel('s')
+
+    def test_meltohz_wrong_values(self):
+        with self.assertRaises(TypeError):
+            meltohz(bytes('ad'))
+
+    def test_meltohz_no_arg(self):
+        with self.assertRaises(TypeError):
+            meltohz()
+
+    def test_meltohz_htk_no_arg(self):
+        with self.assertRaises(TypeError):
+            meltohz_htk()
+
+    def test_hztomel_htk_wrong_values(self):
+        with self.assertRaises(TypeError):
+            hztomel_htk('0')
+
+    def test_hztomel_htk_false(self):
+        assert hztomel(120, htk=False) == hztomel(120)
+
+    def test_meltohz_htk_false(self):
+        assert meltohz(12, htk=False) == meltohz(12)
+
+
+if __name__ == '__main__':
+    from unittest import main
+    main()
diff --git a/python/tests/test_source_channels.py b/python/tests/test_source_channels.py
new file mode 100755 (executable)
index 0000000..eee696d
--- /dev/null
@@ -0,0 +1,93 @@
+#! /usr/bin/env python
+
+"""A brute force test using `sink` to create and write samples to a stereo
+file, then `source` to check the correct content is read from the files."""
+
+import os.path
+import unittest
+import aubio
+import numpy as np
+from numpy.testing import assert_equal
+from utils import get_tmp_sink_path
+
+class aubio_source_test_case(unittest.TestCase):
+
+    def test_read_from_mono(self):
+        out = get_tmp_sink_path()
+        samplerate = 44100
+        hop_size = 256
+        blocks = 10
+        channels = 1
+        write_samples = np.ones([channels, hop_size], dtype=aubio.float_type)
+        write_samples *= .5
+        self.check_write_and_read(samplerate, channels, hop_size, blocks,
+                write_samples)
+
+    def test_read_from_stereo(self):
+        out = get_tmp_sink_path()
+        samplerate = 44100
+        hop_size = 256
+        blocks = 10
+        channels = 1
+        write_samples = np.ones([channels, hop_size], dtype=aubio.float_type)
+        write_samples *= .5
+        self.check_write_and_read(samplerate, channels, hop_size, blocks,
+                write_samples)
+
+    def test_read_from_half_stereo(self):
+        samplerate = 16000
+        channels = 2
+        hop_size = 512
+        blocks = 10
+        write_samples = np.ones([channels, hop_size], dtype=aubio.float_type)
+        write_samples *= .5
+        write_samples[1, :] = 0
+        self.check_write_and_read(samplerate, channels, hop_size, blocks,
+                write_samples)
+
+    def test_read_from_cancelling_channels(self):
+        samplerate = 16000
+        channels = 2
+        hop_size = 512
+        blocks = 10
+        write_samples = np.ones([channels, hop_size], dtype=aubio.float_type)
+        write_samples *= .5
+        write_samples[1] *= -1
+        self.check_write_and_read(samplerate, channels, hop_size, blocks,
+                write_samples)
+
+    def test_read_from_strange_three_channels(self):
+        samplerate = 8000
+        channels = 3
+        hop_size = 123
+        blocks = 10
+        write_samples = np.ones([channels, hop_size], dtype=aubio.float_type)
+        write_samples *= .5
+        write_samples[1, :] = 0
+        self.check_write_and_read(samplerate, channels, hop_size, blocks,
+                write_samples)
+
+    def check_write_and_read(self, samplerate, channels,
+            hop_size, blocks, write_samples):
+        expected_mono = np.sum(write_samples, axis=0)/write_samples.shape[0]
+        out = get_tmp_sink_path()
+        snk = aubio.sink(out, samplerate, channels=channels)
+        for i in range(blocks):
+            snk.do_multi(write_samples, hop_size)
+        # close the sink before reading from it
+        snk.close()
+
+        src = aubio.source(out, samplerate, hop_size)
+        for i in range(blocks):
+            read_samples, read = src.do_multi()
+            assert_equal (read_samples, write_samples)
+            assert_equal (read, hop_size)
+
+        src.seek(0)
+        for i in range(blocks):
+            read_samples, read = src()
+            assert_equal (read, hop_size)
+            assert_equal (read_samples, expected_mono)
+
+if __name__ == '__main__':
+    unittest.main()
index f0be674..3b358bc 100644 (file)
@@ -89,12 +89,10 @@ struct _aubio_source_avcodec_t {
   uint_t read_index;
   sint_t selected_stream;
   uint_t eof;
-  uint_t multi;
 };
 
 // create or re-create the context when _do or _do_multi is called
-void aubio_source_avcodec_reset_resampler(aubio_source_avcodec_t * s,
-    uint_t multi);
+void aubio_source_avcodec_reset_resampler(aubio_source_avcodec_t * s);
 // actually read a frame
 void aubio_source_avcodec_readframe(aubio_source_avcodec_t *s,
     uint_t * read_samples);
@@ -284,13 +282,11 @@ aubio_source_avcodec_t * new_aubio_source_avcodec(const char_t * path,
   s->avCodecCtx = avCodecCtx;
   s->avFrame = avFrame;
 
-  // default to mono output
-  aubio_source_avcodec_reset_resampler(s, 0);
+  aubio_source_avcodec_reset_resampler(s);
 
   if (s->avr == NULL) goto beach;
 
   s->eof = 0;
-  s->multi = 0;
 
   //av_log_set_level(AV_LOG_QUIET);
 
@@ -303,21 +299,17 @@ beach:
   return NULL;
 }
 
-void aubio_source_avcodec_reset_resampler(aubio_source_avcodec_t * s,
-    uint_t multi)
+void aubio_source_avcodec_reset_resampler(aubio_source_avcodec_t * s)
 {
   // create or reset resampler to/from mono/multi-channel
-  if ( (multi != s->multi) || (s->avr == NULL) ) {
+  if ( s->avr == NULL ) {
     int err;
     int64_t input_layout = av_get_default_channel_layout(s->input_channels);
-    uint_t output_channels = multi ? s->input_channels : 1;
-    int64_t output_layout = av_get_default_channel_layout(output_channels);
+    int64_t output_layout = av_get_default_channel_layout(s->input_channels);
 #ifdef HAVE_AVRESAMPLE
     AVAudioResampleContext *avr = avresample_alloc_context();
-    AVAudioResampleContext *oldavr = s->avr;
 #elif defined(HAVE_SWRESAMPLE)
     SwrContext *avr = swr_alloc();
-    SwrContext *oldavr = s->avr;
 #endif /* HAVE_AVRESAMPLE || HAVE_SWRESAMPLE */
 
     av_opt_set_int(avr, "in_channel_layout",  input_layout,              0);
@@ -345,16 +337,6 @@ void aubio_source_avcodec_reset_resampler(aubio_source_avcodec_t * s,
       return;
     }
     s->avr = avr;
-    if (oldavr != NULL) {
-#ifdef HAVE_AVRESAMPLE
-      avresample_close( oldavr );
-#elif defined(HAVE_SWRESAMPLE)
-      swr_close ( oldavr );
-#endif /* HAVE_AVRESAMPLE || HAVE_SWRESAMPLE */
-      av_free ( oldavr );
-      oldavr = NULL;
-    }
-    s->multi = multi;
   }
 }
 
@@ -495,15 +477,18 @@ beach:
 
 void aubio_source_avcodec_do(aubio_source_avcodec_t * s, fvec_t * read_data,
     uint_t * read) {
-  uint_t i;
+  uint_t i, j;
   uint_t end = 0;
   uint_t total_wrote = 0;
-  // switch from multi
-  if (s->multi == 1) aubio_source_avcodec_reset_resampler(s, 0);
   while (total_wrote < s->hop_size) {
     end = MIN(s->read_samples - s->read_index, s->hop_size - total_wrote);
     for (i = 0; i < end; i++) {
-      read_data->data[i + total_wrote] = s->output[i + s->read_index];
+      read_data->data[i + total_wrote] = 0.;
+      for (j = 0; j < s->input_channels; j++) {
+        read_data->data[i + total_wrote] +=
+          s->output[(i + s->read_index) * s->input_channels + j];
+      }
+      read_data->data[i + total_wrote] *= 1./s->input_channels;
     }
     total_wrote += end;
     if (total_wrote < s->hop_size) {
@@ -531,8 +516,6 @@ void aubio_source_avcodec_do_multi(aubio_source_avcodec_t * s,
   uint_t i,j;
   uint_t end = 0;
   uint_t total_wrote = 0;
-  // switch from mono
-  if (s->multi == 0) aubio_source_avcodec_reset_resampler(s, 1);
   while (total_wrote < s->hop_size) {
     end = MIN(s->read_samples - s->read_index, s->hop_size - total_wrote);
     for (j = 0; j < read_data->height; j++) {
index b9de9e5..35755fe 100644 (file)
@@ -387,6 +387,15 @@ fvec_add (fvec_t * o, smpl_t val)
   }
 }
 
+void
+fvec_mul (fvec_t *o, smpl_t val)
+{
+  uint_t j;
+  for (j = 0; j < o->length; j++) {
+    o->data[j] *= val;
+  }
+}
+
 void fvec_adapt_thres(fvec_t * vec, fvec_t * tmp,
     uint_t post, uint_t pre) {
   uint_t length = vec->length, j;
index 56e52ea..4336d7e 100644 (file)
@@ -193,6 +193,14 @@ void fvec_alpha_normalise (fvec_t * v, smpl_t p);
 */
 void fvec_add (fvec_t * v, smpl_t c);
 
+/** multiply each elements of a vector by a scalar
+
+  \param v vector to add constant to
+  \param s constant to scale v with
+
+*/
+void fvec_mul (fvec_t * v, smpl_t s);
+
 /** remove the minimum value of the vector to each elements
 
   \param v vector to remove minimum from
diff --git a/src/musicutils.c b/src/musicutils.c
new file mode 100644 (file)
index 0000000..14ef849
--- /dev/null
@@ -0,0 +1,85 @@
+/*
+  Copyright (C) 2018 Paul Brossier <piem@aubio.org>
+
+  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 <http://www.gnu.org/licenses/>.
+
+*/
+
+#include "aubio_priv.h"
+#include "musicutils.h"
+
+smpl_t
+aubio_hztomel (smpl_t freq)
+{
+  const smpl_t lin_space = 3./200.;
+  const smpl_t split_hz = 1000.;
+  const smpl_t split_mel = split_hz * lin_space;
+  const smpl_t log_space = 27./LOG(6400/1000.);
+  if (freq < 0) {
+    AUBIO_WRN("hztomel: input frequency should be >= 0\n");
+    return 0;
+  }
+  if (freq < split_hz)
+  {
+    return freq * lin_space;
+  } else {
+    return split_mel + log_space * LOG (freq / split_hz);
+  }
+
+}
+
+smpl_t
+aubio_meltohz (smpl_t mel)
+{
+  const smpl_t lin_space = 200./3.;
+  const smpl_t split_hz = 1000.;
+  const smpl_t split_mel = split_hz / lin_space;
+  const smpl_t logSpacing = POW(6400/1000., 1/27.);
+  if (mel < 0) {
+    AUBIO_WRN("meltohz: input mel should be >= 0\n");
+    return 0;
+  }
+  if (mel < split_mel) {
+    return lin_space * mel;
+  } else {
+    return split_hz * POW(logSpacing, mel - split_mel);
+  }
+}
+
+smpl_t
+aubio_hztomel_htk (smpl_t freq)
+{
+  const smpl_t split_hz = 700.;
+  const smpl_t log_space = 1127.;
+  if (freq < 0) {
+    AUBIO_WRN("hztomel_htk: input frequency should be >= 0\n");
+    return 0;
+  }
+  return log_space * LOG (1 + freq / split_hz);
+}
+
+smpl_t
+aubio_meltohz_htk (smpl_t mel)
+{
+  const smpl_t split_hz = 700.;
+  const smpl_t log_space = 1./1127.;
+  if (mel < 0) {
+    AUBIO_WRN("meltohz_htk: input frequency should be >= 0\n");
+    return 0;
+  }
+  return split_hz * ( EXP ( mel * log_space) - 1.);
+}
+
index d9638ac..a659495 100644 (file)
@@ -86,6 +86,105 @@ smpl_t aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize);
 /** convert frequency (Hz) to frequency bin */
 smpl_t aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize);
 
+/** convert frequency (Hz) to mel
+
+  \param freq input frequency, in Hz
+
+  \return output mel
+
+  Converts a scalar from the frequency domain to the mel scale using Slaney
+  Auditory Toolbox's implementation:
+
+  If \f$ f < 1000 \f$, \f$ m = 3 f / 200 \f$.
+
+  If \f$ f >= 1000 \f$, \f$ m = 1000 + 27 \frac{{ln}(f) - ln(1000))}
+  {{ln}(6400) - ln(1000)}
+  \f$
+
+  See also
+  --------
+
+  aubio_meltohz(), aubio_hztomel_htk().
+
+*/
+smpl_t aubio_hztomel (smpl_t freq);
+
+/** convert mel to frequency (Hz)
+
+  \param mel input mel
+
+  \return output frequency, in Hz
+
+  Converts a scalar from the mel scale to the frequency domain using Slaney
+  Auditory Toolbox's implementation:
+
+  If \f$ f < 1000 \f$, \f$ f = 200 m/3 \f$.
+
+  If \f$ f \geq 1000 \f$, \f$ f = 1000 + \left(\frac{6400}{1000}\right)
+  ^{\frac{m - 1000}{27}} \f$
+
+  See also
+  --------
+
+  aubio_hztomel(), aubio_meltohz_htk().
+
+  References
+  ----------
+
+  Malcolm Slaney, *Auditory Toolbox Version 2, Technical Report #1998-010*
+  https://engineering.purdue.edu/~malcolm/interval/1998-010/
+
+*/
+smpl_t aubio_meltohz (smpl_t mel);
+
+/** convert frequency (Hz) to mel
+
+  \param freq input frequency, in Hz
+
+  \return output mel
+
+  Converts a scalar from the frequency domain to the mel scale, using the
+  equation defined by O'Shaughnessy, as implemented in the HTK speech
+  recognition toolkit:
+
+  \f$ m = 1127 + ln(1 + \frac{f}{700}) \f$
+
+  See also
+  --------
+
+  aubio_meltohz_htk(), aubio_hztomel().
+
+  References
+  ----------
+
+  Douglas O'Shaughnessy (1987). *Speech communication: human and machine*.
+  Addison-Wesley. p. 150. ISBN 978-0-201-16520-3.
+
+  HTK Speech Recognition Toolkit: http://htk.eng.cam.ac.uk/
+
+ */
+smpl_t aubio_hztomel_htk (smpl_t freq);
+
+/** convert mel to frequency (Hz)
+
+  \param mel input mel
+
+  \return output frequency, in Hz
+
+  Converts a scalar from the mel scale to the frequency domain, using the
+  equation defined by O'Shaughnessy, as implemented in the HTK speech
+  recognition toolkit:
+
+  \f$ f = 700 * {e}^\left(\frac{f}{1127} - 1\right) \f$
+
+  See also
+  --------
+
+  aubio_hztomel_htk(), aubio_meltohz().
+
+*/
+smpl_t aubio_meltohz_htk (smpl_t mel);
+
 /** convert frequency (Hz) to midi value (0-128) */
 smpl_t aubio_freqtomidi (smpl_t freq);
 
index 0ca467c..598fe22 100644 (file)
@@ -89,11 +89,12 @@ pthread_mutex_t aubio_fftw_mutex = PTHREAD_MUTEX_INITIALIZER;
 #define aubio_vDSP_zvphas              vDSP_zvphas
 #define aubio_vDSP_vsadd               vDSP_vsadd
 #define aubio_vDSP_vsmul               vDSP_vsmul
-#define aubio_vDSP_create_fftsetup     vDSP_create_fftsetup
-#define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetup
 #define aubio_DSPComplex               DSPComplex
 #define aubio_DSPSplitComplex          DSPSplitComplex
-#define aubio_FFTSetup                 FFTSetup
+#define aubio_vDSP_DFT_Setup           vDSP_DFT_Setup
+#define aubio_vDSP_DFT_zrop_CreateSetup vDSP_DFT_zrop_CreateSetup
+#define aubio_vDSP_DFT_Execute         vDSP_DFT_Execute
+#define aubio_vDSP_DFT_DestroySetup    vDSP_DFT_DestroySetup
 #define aubio_vvsqrt                   vvsqrtf
 #else
 #define aubio_vDSP_ctoz                vDSP_ctozD
@@ -103,11 +104,12 @@ pthread_mutex_t aubio_fftw_mutex = PTHREAD_MUTEX_INITIALIZER;
 #define aubio_vDSP_zvphas              vDSP_zvphasD
 #define aubio_vDSP_vsadd               vDSP_vsaddD
 #define aubio_vDSP_vsmul               vDSP_vsmulD
-#define aubio_vDSP_create_fftsetup     vDSP_create_fftsetupD
-#define aubio_vDSP_destroy_fftsetup    vDSP_destroy_fftsetupD
 #define aubio_DSPComplex               DSPDoubleComplex
 #define aubio_DSPSplitComplex          DSPDoubleSplitComplex
-#define aubio_FFTSetup                 FFTSetupD
+#define aubio_vDSP_DFT_Setup           vDSP_DFT_SetupD
+#define aubio_vDSP_DFT_zrop_CreateSetup vDSP_DFT_zrop_CreateSetupD
+#define aubio_vDSP_DFT_Execute         vDSP_DFT_ExecuteD
+#define aubio_vDSP_DFT_DestroySetup    vDSP_DFT_DestroySetupD
 #define aubio_vvsqrt                   vvsqrt
 #endif /* HAVE_AUBIO_DOUBLE */
 
@@ -152,8 +154,8 @@ struct _aubio_fft_t {
   fft_data_t * specdata; /* complex spectral data */
 
 #elif defined HAVE_ACCELERATE  // using ACCELERATE
-  int log2fftsize;
-  aubio_FFTSetup fftSetup;
+  aubio_vDSP_DFT_Setup fftSetupFwd;
+  aubio_vDSP_DFT_Setup fftSetupBwd;
   aubio_DSPSplitComplex spec;
   smpl_t *in, *out;
 
@@ -210,15 +212,32 @@ aubio_fft_t * new_aubio_fft (uint_t winsize) {
   }
 
 #elif defined HAVE_ACCELERATE  // using ACCELERATE
+  {
+    uint_t radix = winsize;
+    uint_t order = 0;
+    while ((radix / 2) * 2 == radix) {
+      radix /= 2;
+      order++;
+    }
+    if (order < 4 || (radix != 1 && radix != 3 && radix != 5 && radix != 15)) {
+      AUBIO_ERR("fft: vDSP/Accelerate supports FFT with sizes = "
+          "f * 2 ** n, where n > 4 and f in [1, 3, 5, 15], but requested %d. "
+          "Use the closest power of two, or try recompiling aubio with "
+          "--enable-fftw3.\n", winsize);
+      goto beach;
+    }
+  }
   s->winsize = winsize;
   s->fft_size = winsize;
   s->compspec = new_fvec(winsize);
-  s->log2fftsize = aubio_power_of_two_order(s->fft_size);
   s->in = AUBIO_ARRAY(smpl_t, s->fft_size);
   s->out = AUBIO_ARRAY(smpl_t, s->fft_size);
   s->spec.realp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
   s->spec.imagp = AUBIO_ARRAY(smpl_t, s->fft_size/2);
-  s->fftSetup = aubio_vDSP_create_fftsetup(s->log2fftsize, FFT_RADIX2);
+  s->fftSetupFwd = aubio_vDSP_DFT_zrop_CreateSetup(NULL,
+      s->fft_size, vDSP_DFT_FORWARD);
+  s->fftSetupBwd = aubio_vDSP_DFT_zrop_CreateSetup(s->fftSetupFwd,
+      s->fft_size, vDSP_DFT_INVERSE);
 
 #elif defined HAVE_INTEL_IPP  // using Intel IPP
   const IppHintAlgorithm qualityHint = ippAlgHintAccurate; // OR ippAlgHintFast;
@@ -292,7 +311,8 @@ void del_aubio_fft(aubio_fft_t * s) {
 #elif defined HAVE_ACCELERATE // using ACCELERATE
   AUBIO_FREE(s->spec.realp);
   AUBIO_FREE(s->spec.imagp);
-  aubio_vDSP_destroy_fftsetup(s->fftSetup);
+  aubio_vDSP_DFT_DestroySetup(s->fftSetupBwd);
+  aubio_vDSP_DFT_DestroySetup(s->fftSetupFwd);
 
 #elif defined HAVE_INTEL_IPP  // using Intel IPP
   ippFree(s->memSpec);
@@ -350,7 +370,8 @@ void aubio_fft_do_complex(aubio_fft_t * s, const fvec_t * input, fvec_t * compsp
   // convert real data to even/odd format used in vDSP
   aubio_vDSP_ctoz((aubio_DSPComplex*)s->in, 2, &s->spec, 1, s->fft_size/2);
   // compute the FFT
-  aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_FORWARD);
+  aubio_vDSP_DFT_Execute(s->fftSetupFwd, s->spec.realp, s->spec.imagp,
+      s->spec.realp, s->spec.imagp);
   // convert from vDSP complex split to [ r0, r1, ..., rN, iN-1, .., i2, i1]
   compspec->data[0] = s->spec.realp[0];
   compspec->data[s->fft_size / 2] = s->spec.imagp[0];
@@ -418,7 +439,8 @@ void aubio_fft_rdo_complex(aubio_fft_t * s, const fvec_t * compspec, fvec_t * ou
   // convert to split complex format used in vDSP
   aubio_vDSP_ctoz((aubio_DSPComplex*)s->out, 2, &s->spec, 1, s->fft_size/2);
   // compute the FFT
-  aubio_vDSP_fft_zrip(s->fftSetup, &s->spec, 1, s->log2fftsize, FFT_INVERSE);
+  aubio_vDSP_DFT_Execute(s->fftSetupBwd, s->spec.realp, s->spec.imagp,
+      s->spec.realp, s->spec.imagp);
   // convert result to real output
   aubio_vDSP_ztoc(&s->spec, 1, (aubio_DSPComplex*)output->data, 2, s->fft_size/2);
   // apply scaling
@@ -493,11 +515,22 @@ void aubio_fft_get_phas(const fvec_t * compspec, cvec_t * spectrum) {
         compspec->data[i]);
   }
 #endif
-  if (compspec->data[compspec->length/2] < 0) {
-    spectrum->phas[spectrum->length - 1] = PI;
+#ifdef HAVE_FFTW3
+  // for even length only, make sure last element is 0 or PI
+  if (2 * (compspec->length / 2) == compspec->length) {
+#endif
+    if (compspec->data[compspec->length/2] < 0) {
+      spectrum->phas[spectrum->length - 1] = PI;
+    } else {
+      spectrum->phas[spectrum->length - 1] = 0.;
+    }
+#ifdef HAVE_FFTW3
   } else {
-    spectrum->phas[spectrum->length - 1] = 0.;
+    i = spectrum->length - 1;
+    spectrum->phas[i] = ATAN2(compspec->data[compspec->length-i],
+        compspec->data[i]);
   }
+#endif
 }
 
 void aubio_fft_get_norm(const fvec_t * compspec, cvec_t * spectrum) {
@@ -507,8 +540,19 @@ void aubio_fft_get_norm(const fvec_t * compspec, cvec_t * spectrum) {
     spectrum->norm[i] = SQRT(SQR(compspec->data[i])
         + SQR(compspec->data[compspec->length - i]) );
   }
-  spectrum->norm[spectrum->length-1] =
-    ABS(compspec->data[compspec->length/2]);
+#ifdef HAVE_FFTW3
+  // for even length, make sure last element is > 0
+  if (2 * (compspec->length / 2) == compspec->length) {
+#endif
+    spectrum->norm[spectrum->length-1] =
+      ABS(compspec->data[compspec->length/2]);
+#ifdef HAVE_FFTW3
+  } else {
+    i = spectrum->length - 1;
+    spectrum->norm[i] = SQRT(SQR(compspec->data[i])
+        + SQR(compspec->data[compspec->length - i]) );
+  }
+#endif
 }
 
 void aubio_fft_get_imag(const cvec_t * spectrum, fvec_t * compspec) {
index 0323700..5497a4f 100644 (file)
@@ -23,6 +23,7 @@
 #include "fvec.h"
 #include "fmat.h"
 #include "cvec.h"
+#include "vecutils.h"
 #include "spectral/filterbank.h"
 #include "mathutils.h"
 
@@ -32,6 +33,8 @@ struct _aubio_filterbank_t
   uint_t win_s;
   uint_t n_filters;
   fmat_t *filters;
+  smpl_t norm;
+  smpl_t power;
 };
 
 aubio_filterbank_t *
@@ -45,6 +48,10 @@ new_aubio_filterbank (uint_t n_filters, uint_t win_s)
   /* allocate filter tables, a matrix of length win_s and of height n_filters */
   fb->filters = new_fmat (n_filters, win_s / 2 + 1);
 
+  fb->norm = 1;
+
+  fb->power = 1;
+
   return fb;
 }
 
@@ -67,6 +74,8 @@ aubio_filterbank_do (aubio_filterbank_t * f, const cvec_t * in, fvec_t * out)
   tmp.length = in->length;
   tmp.data = in->norm;
 
+  if (f->power != 1.) fvec_pow(&tmp, f->power);
+
   fmat_vecmul(f->filters, &tmp, out);
 
   return;
@@ -84,3 +93,30 @@ aubio_filterbank_set_coeffs (aubio_filterbank_t * f, const fmat_t * filter_coeff
   fmat_copy(filter_coeffs, f->filters);
   return 0;
 }
+
+uint_t
+aubio_filterbank_set_norm (aubio_filterbank_t *f, smpl_t norm)
+{
+  if (norm != 0 && norm != 1) return AUBIO_FAIL;
+  f->norm = norm;
+  return AUBIO_OK;
+}
+
+smpl_t
+aubio_filterbank_get_norm (aubio_filterbank_t *f)
+{
+  return f->norm;
+}
+
+uint_t
+aubio_filterbank_set_power (aubio_filterbank_t *f, smpl_t power)
+{
+  f->power = power;
+  return AUBIO_OK;
+}
+
+smpl_t
+aubio_filterbank_get_power (aubio_filterbank_t *f)
+{
+  return f->norm;
+}
index 769b5e7..876b14d 100644 (file)
@@ -83,6 +83,47 @@ fmat_t *aubio_filterbank_get_coeffs (const aubio_filterbank_t * f);
  */
 uint_t aubio_filterbank_set_coeffs (aubio_filterbank_t * f, const fmat_t * filters);
 
+/** set norm parameter
+
+  \param f filterbank object, as returned by new_aubio_filterbank()
+  \param norm `1` to norm the filters, `0` otherwise.
+
+  If set to `0`, the filters will not be normalized. If set to `1`,
+  each filter will be normalized to one. Defaults to `1`.
+
+  This function should be called *before* setting the filters with one of
+  aubio_filterbank_set_triangle_bands(), aubio_filterbank_set_mel_coeffs(),
+  aubio_filterbank_set_mel_coeffs_htk(), or
+  aubio_filterbank_set_mel_coeffs_slaney().
+
+ */
+uint_t aubio_filterbank_set_norm (aubio_filterbank_t *f, smpl_t norm);
+
+/** get norm parameter
+
+  \param f filterbank object, as returned by new_aubio_filterbank()
+  \returns `1` if norm is set, `0` otherwise. Defaults to `1`.
+
+ */
+smpl_t aubio_filterbank_get_norm (aubio_filterbank_t *f);
+
+/** set power parameter
+
+  \param f filterbank object, as returned by new_aubio_filterbank()
+  \param power Raise norm of the input spectrum norm to this power before
+  computing filterbank.  Defaults to `1`.
+
+ */
+uint_t aubio_filterbank_set_power (aubio_filterbank_t *f, smpl_t power);
+
+/** get power parameter
+
+  \param f filterbank object, as returned by new_aubio_filterbank()
+  \return current power parameter. Defaults to `1`.
+
+ */
+smpl_t aubio_filterbank_get_power (aubio_filterbank_t *f);
+
 #ifdef __cplusplus
 }
 #endif
index 8a8d85a..37713a3 100644 (file)
@@ -90,9 +90,13 @@ aubio_filterbank_set_triangle_bands (aubio_filterbank_t * fb,
   }
 
   /* compute triangle heights so that each triangle has unit area */
-  for (fn = 0; fn < n_filters; fn++) {
-    triangle_heights->data[fn] =
-        2. / (upper_freqs->data[fn] - lower_freqs->data[fn]);
+  if (aubio_filterbank_get_norm(fb)) {
+    for (fn = 0; fn < n_filters; fn++) {
+      triangle_heights->data[fn] =
+          2. / (upper_freqs->data[fn] - lower_freqs->data[fn]);
+    }
+  } else {
+    fvec_ones (triangle_heights);
   }
 
   /* fill fft_freqs lookup table, which assigns the frequency in hz to each bin */
@@ -117,9 +121,8 @@ aubio_filterbank_set_triangle_bands (aubio_filterbank_t * fb,
     }
 
     /* compute positive slope step size */
-    riseInc =
-        triangle_heights->data[fn] /
-        (center_freqs->data[fn] - lower_freqs->data[fn]);
+    riseInc = triangle_heights->data[fn]
+      / (center_freqs->data[fn] - lower_freqs->data[fn]);
 
     /* compute coefficients in positive slope */
     for (; bin < win_s - 1; bin++) {
@@ -133,9 +136,8 @@ aubio_filterbank_set_triangle_bands (aubio_filterbank_t * fb,
     }
 
     /* compute negative slope step size */
-    downInc =
-        triangle_heights->data[fn] /
-        (upper_freqs->data[fn] - center_freqs->data[fn]);
+    downInc = triangle_heights->data[fn]
+      / (upper_freqs->data[fn] - center_freqs->data[fn]);
 
     /* compute coefficents in negative slope */
     for (; bin < win_s - 1; bin++) {
@@ -168,23 +170,27 @@ uint_t
 aubio_filterbank_set_mel_coeffs_slaney (aubio_filterbank_t * fb,
     smpl_t samplerate)
 {
-  uint_t retval;
-
   /* Malcolm Slaney parameters */
-  smpl_t lowestFrequency = 133.3333;
-  smpl_t linearSpacing = 66.66666666;
-  smpl_t logSpacing = 1.0711703;
+  const smpl_t lowestFrequency = 133.3333;
+  const smpl_t linearSpacing = 66.66666666;
+  const smpl_t logSpacing = 1.0711703;
 
-  uint_t linearFilters = 13;
-  uint_t logFilters = 27;
-  uint_t n_filters = linearFilters + logFilters;
-
-  uint_t fn;                    /* filter counter */
+  const uint_t linearFilters = 13;
+  const uint_t logFilters = 27;
+  const uint_t n_filters = linearFilters + logFilters;
 
+  uint_t fn, retval;
   smpl_t lastlinearCF;
 
   /* buffers to compute filter frequencies */
-  fvec_t *freqs = new_fvec (n_filters + 2);
+  fvec_t *freqs;
+
+  if (samplerate <= 0) {
+    AUBIO_ERR("filterbank: set_mel_coeffs_slaney samplerate should be > 0\n");
+    return AUBIO_FAIL;
+  }
+
+  freqs = new_fvec (n_filters + 2);
 
   /* first step: fill all the linear filter frequencies */
   for (fn = 0; fn < linearFilters; fn++) {
@@ -206,3 +212,87 @@ aubio_filterbank_set_mel_coeffs_slaney (aubio_filterbank_t * fb,
 
   return retval;
 }
+
+static uint_t aubio_filterbank_check_freqs (aubio_filterbank_t *fb UNUSED,
+    smpl_t samplerate, smpl_t *freq_min, smpl_t *freq_max)
+{
+  if (samplerate <= 0) {
+    AUBIO_ERR("filterbank: set_mel_coeffs samplerate should be > 0\n");
+    return AUBIO_FAIL;
+  }
+  if (*freq_max < 0) {
+    AUBIO_ERR("filterbank: set_mel_coeffs freq_max should be > 0\n");
+    return AUBIO_FAIL;
+  } else if (*freq_max == 0) {
+    *freq_max = samplerate / 2.;
+  }
+  if (*freq_min < 0) {
+    AUBIO_ERR("filterbank: set_mel_coeffs freq_min should be > 0\n");
+    return AUBIO_FAIL;
+  }
+  return AUBIO_OK;
+}
+
+uint_t
+aubio_filterbank_set_mel_coeffs (aubio_filterbank_t * fb, smpl_t samplerate,
+    smpl_t freq_min, smpl_t freq_max)
+{
+  uint_t m, retval;
+  smpl_t start = freq_min, end = freq_max, step;
+  fvec_t *freqs;
+  fmat_t *coeffs = aubio_filterbank_get_coeffs(fb);
+  uint_t n_bands = coeffs->height;
+
+  if (aubio_filterbank_check_freqs(fb, samplerate, &start, &end)) {
+    return AUBIO_FAIL;
+  }
+
+  start = aubio_hztomel(start);
+  end = aubio_hztomel(end);
+
+  freqs = new_fvec(n_bands + 2);
+  step = (end - start) / (n_bands + 1);
+
+  for (m = 0; m < n_bands + 2; m++)
+  {
+    freqs->data[m] = MIN(aubio_meltohz(start + step * m), samplerate/2.);
+  }
+
+  retval = aubio_filterbank_set_triangle_bands (fb, freqs, samplerate);
+
+  /* destroy vector used to store frequency limits */
+  del_fvec (freqs);
+  return retval;
+}
+
+uint_t
+aubio_filterbank_set_mel_coeffs_htk (aubio_filterbank_t * fb, smpl_t samplerate,
+    smpl_t freq_min, smpl_t freq_max)
+{
+  uint_t m, retval;
+  smpl_t start = freq_min, end = freq_max, step;
+  fvec_t *freqs;
+  fmat_t *coeffs = aubio_filterbank_get_coeffs(fb);
+  uint_t n_bands = coeffs->height;
+
+  if (aubio_filterbank_check_freqs(fb, samplerate, &start, &end)) {
+    return AUBIO_FAIL;
+  }
+
+  start = aubio_hztomel_htk(start);
+  end = aubio_hztomel_htk(end);
+
+  freqs = new_fvec (n_bands + 2);
+  step = (end - start) / (n_bands + 1);
+
+  for (m = 0; m < n_bands + 2; m++)
+  {
+    freqs->data[m] = MIN(aubio_meltohz_htk(start + step * m), samplerate/2.);
+  }
+
+  retval = aubio_filterbank_set_triangle_bands (fb, freqs, samplerate);
+
+  /* destroy vector used to store frequency limits */
+  del_fvec (freqs);
+  return retval;
+}
index ffa4948..06bbf6c 100644 (file)
@@ -55,17 +55,63 @@ uint_t aubio_filterbank_set_triangle_bands (aubio_filterbank_t * fb,
 /** filterbank initialization for Mel filters using Slaney's coefficients
 
   \param fb filterbank object
-  \param samplerate audio sampling rate
+  \param samplerate audio sampling rate, in Hz
+
+  The filter coefficients are built to match exactly Malcolm Slaney's Auditory
+  Toolbox implementation (see file mfcc.m). The number of filters should be 40.
 
-  The filter coefficients are built according to Malcolm Slaney's Auditory
-  Toolbox, available online at the following address (see file mfcc.m):
+  References
+  ----------
 
+  Malcolm Slaney, *Auditory Toolbox Version 2, Technical Report #1998-010*
   https://engineering.purdue.edu/~malcolm/interval/1998-010/
 
 */
 uint_t aubio_filterbank_set_mel_coeffs_slaney (aubio_filterbank_t * fb,
     smpl_t samplerate);
 
+/** Mel filterbank initialization
+
+  \param fb filterbank object
+  \param samplerate audio sampling rate
+  \param fmin start frequency, in Hz
+  \param fmax end frequency, in Hz
+
+  The filterbank will be initialized with bands linearly spaced in the mel
+  scale, from `fmin` to `fmax`.
+
+  References
+  ----------
+
+  Malcolm Slaney, *Auditory Toolbox Version 2, Technical Report #1998-010*
+  https://engineering.purdue.edu/~malcolm/interval/1998-010/
+
+*/
+uint_t aubio_filterbank_set_mel_coeffs(aubio_filterbank_t * fb,
+    smpl_t samplerate, smpl_t fmin, smpl_t fmax);
+
+/** Mel filterbank initialization
+
+  \param fb filterbank object
+  \param samplerate audio sampling rate
+  \param fmin start frequency, in Hz
+  \param fmax end frequency, in Hz
+
+  The bank of filters will be initalized to to cover linearly spaced bands in
+  the Htk mel scale, from `fmin` to `fmax`.
+
+  References
+  ----------
+
+  Douglas O'Shaughnessy (1987). *Speech communication: human and machine*.
+  Addison-Wesley. p. 150. ISBN 978-0-201-16520-3.
+
+  HTK Speech Recognition Toolkit: http://htk.eng.cam.ac.uk/
+
+*/
+uint_t aubio_filterbank_set_mel_coeffs_htk(aubio_filterbank_t * fb,
+    smpl_t samplerate, smpl_t fmin, smpl_t fmax);
+
 #ifdef __cplusplus
 }
 #endif
index 32f0beb..fea7c7c 100644 (file)
@@ -51,6 +51,7 @@ struct _aubio_mfcc_t
   aubio_dct_t *dct;
   fvec_t *output;
 #endif
+  smpl_t scale;
 };
 
 
@@ -74,7 +75,11 @@ new_aubio_mfcc (uint_t win_s, uint_t n_filters, uint_t n_coefs,
 
   /* filterbank allocation */
   mfcc->fb = new_aubio_filterbank (n_filters, mfcc->win_s);
-  aubio_filterbank_set_mel_coeffs_slaney (mfcc->fb, samplerate);
+  if (n_filters == 40)
+    aubio_filterbank_set_mel_coeffs_slaney (mfcc->fb, samplerate);
+  else
+    aubio_filterbank_set_mel_coeffs(mfcc->fb, samplerate,
+        0, samplerate/2.);
 
   /* allocating buffers */
   mfcc->in_dct = new_fvec (n_filters);
@@ -97,6 +102,8 @@ new_aubio_mfcc (uint_t win_s, uint_t n_filters, uint_t n_coefs,
   mfcc->output = new_fvec (n_filters);
 #endif
 
+  mfcc->scale = 1.;
+
   return mfcc;
 }
 
@@ -127,14 +134,14 @@ aubio_mfcc_do (aubio_mfcc_t * mf, const cvec_t * in, fvec_t * out)
 #ifndef HAVE_SLOW_DCT
   fvec_t tmp;
 #endif
+
   /* compute filterbank */
   aubio_filterbank_do (mf->fb, in, mf->in_dct);
 
   /* compute log10 */
   fvec_log10 (mf->in_dct);
 
-  /* raise power */
-  //fvec_pow (mf->in_dct, 3.);
+  if (mf->scale != 1) fvec_mul (mf->in_dct, mf->scale);
 
   /* compute mfccs */
 #if defined(HAVE_SLOW_DCT)
@@ -150,3 +157,43 @@ aubio_mfcc_do (aubio_mfcc_t * mf, const cvec_t * in, fvec_t * out)
 
   return;
 }
+
+uint_t aubio_mfcc_set_power (aubio_mfcc_t *mf, smpl_t power)
+{
+  return aubio_filterbank_set_power(mf->fb, power);
+}
+
+uint_t aubio_mfcc_get_power (aubio_mfcc_t *mf)
+{
+  return aubio_filterbank_get_power(mf->fb);
+}
+
+uint_t aubio_mfcc_set_scale (aubio_mfcc_t *mf, smpl_t scale)
+{
+  mf->scale = scale;
+  return AUBIO_OK;
+}
+
+uint_t aubio_mfcc_get_scale (aubio_mfcc_t *mf)
+{
+  return mf->scale;
+}
+
+uint_t aubio_mfcc_set_mel_coeffs (aubio_mfcc_t *mf, smpl_t freq_min,
+    smpl_t freq_max)
+{
+  return aubio_filterbank_set_mel_coeffs(mf->fb, mf->samplerate,
+      freq_min, freq_max);
+}
+
+uint_t aubio_mfcc_set_mel_coeffs_htk (aubio_mfcc_t *mf, smpl_t freq_min,
+    smpl_t freq_max)
+{
+  return aubio_filterbank_set_mel_coeffs_htk(mf->fb, mf->samplerate,
+      freq_min, freq_max);
+}
+
+uint_t aubio_mfcc_set_mel_coeffs_slaney (aubio_mfcc_t *mf)
+{
+  return aubio_filterbank_set_mel_coeffs_slaney (mf->fb, mf->samplerate);
+}
index ba2c7f2..ffe6cb6 100644 (file)
@@ -73,6 +73,100 @@ void del_aubio_mfcc (aubio_mfcc_t * mf);
 */
 void aubio_mfcc_do (aubio_mfcc_t * mf, const cvec_t * in, fvec_t * out);
 
+/** set power parameter
+
+  \param mf mfcc object, as returned by new_aubio_mfcc()
+  \param power Raise norm of the input spectrum norm to this power before
+  computing filterbank.  Defaults to `1`.
+
+  See aubio_filterbank_set_power().
+
+ */
+uint_t aubio_mfcc_set_power (aubio_mfcc_t *mf, smpl_t power);
+
+/** get power parameter
+
+  \param mf mfcc object, as returned by new_aubio_mfcc()
+  \return current power parameter. Defaults to `1`.
+
+  See aubio_filterbank_get_power().
+
+ */
+uint_t aubio_mfcc_get_power (aubio_mfcc_t *mf);
+
+/** set scaling parameter
+
+  \param mf mfcc object, as returned by new_aubio_mfcc()
+  \param scale Scaling value to apply.
+
+  Scales the output of the filterbank after taking its logarithm and before
+  computing the DCT. Defaults to `1`.
+
+*/
+uint_t aubio_mfcc_set_scale (aubio_mfcc_t *mf, smpl_t scale);
+
+/** get scaling parameter
+
+  \param mf mfcc object, as returned by new_aubio_mfcc()
+  \return current scaling parameter. Defaults to `1`.
+
+ */
+uint_t aubio_mfcc_get_scale (aubio_mfcc_t *mf);
+
+/** Mel filterbank initialization
+
+  \param mf mfcc object
+  \param fmin start frequency, in Hz
+  \param fmax end frequency, in Hz
+
+  The filterbank will be initialized with bands linearly spaced in the mel
+  scale, from `fmin` to `fmax`.
+
+  See also
+  --------
+
+  aubio_filterbank_set_mel_coeffs()
+
+*/
+uint_t aubio_mfcc_set_mel_coeffs (aubio_mfcc_t *mf,
+        smpl_t fmin, smpl_t fmax);
+
+/** Mel filterbank initialization
+
+  \param mf mfcc object
+  \param fmin start frequency, in Hz
+  \param fmax end frequency, in Hz
+
+  The bank of filters will be initalized to to cover linearly spaced bands in
+  the Htk mel scale, from `fmin` to `fmax`.
+
+  See also
+  --------
+
+  aubio_filterbank_set_mel_coeffs_htk()
+
+*/
+uint_t aubio_mfcc_set_mel_coeffs_htk (aubio_mfcc_t *mf,
+        smpl_t fmin, smpl_t fmax);
+
+/** Mel filterbank initialization (Auditory Toolbox's parameters)
+
+  \param mf mfcc object
+  \param samplerate audio sampling rate, in Hz
+
+  The filter coefficients are built to match exactly Malcolm Slaney's Auditory
+  Toolbox implementation. The number of filters should be 40.
+
+  This is the default filterbank when `mf` was created with `n_filters = 40`.
+
+  See also
+  --------
+
+  aubio_filterbank_set_mel_coeffs_slaney()
+
+*/
+uint_t aubio_mfcc_set_mel_coeffs_slaney (aubio_mfcc_t *mf);
+
 #ifdef __cplusplus
 }
 #endif
index 543226d..8b487c2 100644 (file)
@@ -11,6 +11,15 @@ int main (void)
   // create filterbank object
   aubio_filterbank_t *o = new_aubio_filterbank (n_filters, win_s);
 
+  smpl_t power = aubio_filterbank_get_power(o);
+  smpl_t norm = aubio_filterbank_get_norm(o);
+  if (aubio_filterbank_set_power(o, power)) {
+    return 1;
+  }
+  if (aubio_filterbank_set_norm(o, norm)) {
+    return 1;
+  }
+
   // apply filterbank ten times
   uint_t n = 10;
   while (n) {