Merge pull request #44 from r9y9/mcepalpha

Add mcepalpha

NEWS.md

@@ -2,6 +2,7 @@
 
 ## v0.1.6
 
+- Add `mcepalpha`. [#43]
 - Add `mc2e`. [#42]
 - Add `sp2mc` and `mc2sp`. [#41]
 
@@ -61,3 +62,4 @@
 [#39]: https://github.com/r9y9/pysptk/pull/39
 [#41]: https://github.com/r9y9/pysptk/pull/41
 [#42]: https://github.com/r9y9/pysptk/pull/42
+[#43]: https://github.com/r9y9/pysptk/issues/43

pysptk/util.py

@@ -12,11 +12,20 @@
 
     example_audio_file
 
+
+Mel-cepstrum analysis
+---------------------
+
+.. autosummary::
+    :toctree: generated/
+
+    mcepalpha
 """
 
 from __future__ import division, print_function, absolute_import
 
 import pkg_resources
+import numpy as np
 
 # 16kHz, 16bit example audio from cmu_us_awb_arctic
 # see COPYING for the license of the audio file.
@@ -63,3 +72,68 @@ def example_audio_file():
     """
 
     return pkg_resources.resource_filename(__name__, EXAMPLE_AUDIO)
+
+
+def mcepalpha(fs, start=0.0, stop=1.0, step=0.001, num_points=1000):
+    """Compute appropriate frequency warping parameter given a sampling frequency
+
+    It would be useful to determine alpha parameter in mel-cepstrum analysis.
+
+    The code is traslated from https://bitbucket.org/happyalu/mcep_alpha_calc.
+
+    Parameters
+    ----------
+    fs : int
+        Sampling frequency
+
+    start : float
+        start value that will be passed to numpy.arange. Default is 0.0.
+
+    stop : float
+        stop value that will be passed to numpy.arange. Default is 1.0.
+
+    step : float
+        step value that will be passed to numpy.arange. Default is 0.001.
+
+    num_points : int
+        Number of points used in approximating mel-scale vectors in fixed-
+        length.
+
+    Returns
+    -------
+    alpha : float
+        frequency warping paramter (offen denoted by alpha)
+
+    See Also
+    --------
+    pysptk.sptk.mcep
+    pysptk.sptk.mgcep
+
+    """
+    alpha_candidates = np.arange(start, stop, step)
+    mel = _melscale_vector(fs, num_points)
+    distances = [rms_distance(mel, _warping_vector(alpha, num_points)) for
+                 alpha in alpha_candidates]
+    return alpha_candidates[np.argmin(distances)]
+
+
+def _melscale_vector(fs, length):
+    step = (fs / 2.0) / length
+    melscalev = 1000.0 / np.log(2) * np.log(1 +
+                                            step * np.arange(0, length) / 1000.0)
+    return melscalev / melscalev[-1]
+
+
+def _warping_vector(alpha, length):
+    step = np.pi / length
+    omega = step * np.arange(0, length)
+    num = (1 - alpha * alpha) * np.sin(omega)
+    den = (1 + alpha * alpha) * np.cos(omega) - 2 * alpha
+    warpfreq = np.arctan(num / den)
+    warpfreq[warpfreq < 0] += np.pi
+    return warpfreq / warpfreq[-1]
+
+
+def rms_distance(v1, v2):
+    d = v1 - v2
+    return np.sum(np.abs(d * d)) / len(v1)

tests/test_utils.py

@@ -3,6 +3,7 @@
 import numpy as np
 import pysptk
 from nose.tools import raises
+from pysptk.util import mcepalpha
 
 
 def test_assert_gamma():
@@ -58,3 +59,12 @@ def test_example_audio_file():
     from os.path import exists
     path = pysptk.util.example_audio_file()
     assert exists(path)
+
+
+def test_mcepalpha():
+    assert np.isclose(mcepalpha(8000), 0.312)
+    assert np.isclose(mcepalpha(11025), 0.357)
+    assert np.isclose(mcepalpha(16000), 0.41)
+    assert np.isclose(mcepalpha(22050), 0.455)
+    assert np.isclose(mcepalpha(44100),  0.544)
+    assert np.isclose(mcepalpha(48000), 0.554)