Merge fe0e6f1 into e555bed

docs/changelog.rst

@@ -4,6 +4,9 @@ Change log
 v0.1.14 <2018-xx-xx>
 --------------------
 
+- `#62`_, `#64`_: **IMPORTANT**: Add ``use_scipy`` option to ``levdur`` and ``lpc``, and set it to True by default to improve numerical stability. The results should be same regardless to the option, but this would be technically a breaking change.
+- `#64`_: Add ``acorr``.
+
 v0.1.13 <2018-11-19>
 --------------------
 
@@ -15,7 +18,7 @@ v0.1.13 <2018-11-19>
 v0.1.12 <2018-10-27>
 --------------------
 
-- #63`_: Fix lpc2lsp bug, add lsp2lpc function. Add regression tests for those.
+- `#63`_: Fix lpc2lsp bug, add lsp2lpc function. Add regression tests for those.
 
 v0.1.11 <2018-02-05>
 --------------------
@@ -118,3 +121,5 @@ v0.1.0 <2015-09-05>
 .. _#54: https://github.com/r9y9/pysptk/pull/54
 .. _#55: https://github.com/r9y9/pysptk/issues/55
 .. _#63: https://github.com/r9y9/pysptk/pull/63
+.. _#62: https://github.com/r9y9/pysptk/issues/62
+.. _#64: https://github.com/r9y9/pysptk/pull/64

pysptk/_sptk.pxd

@@ -21,6 +21,8 @@ cdef extern from "SPTK.h":
     int _mseq "mseq"()
 
 
+    void _acorr "acorr"(double *x, int l, double *r, const int np)
+
     # Adaptive mel-generalized cepstrum analysis
     double _acep "acep"(double x, double *c, const int m, const double lambda_coef,
                         const double step, const double tau, const int pd,
@@ -32,7 +34,6 @@ cdef extern from "SPTK.h":
                           const double lambda_coef, const double step, const double tau,
                           const int pd, const double eps);
 
-
     # Mel-generalized cepstrum analysis
     int _mcep "mcep"(double *xw, const int flng, double *mc, const int m, const double a,
                      const int itr1, const int itr2, const double dd, const int etype,

pysptk/_sptk.pyx

@@ -29,6 +29,13 @@ def mseq():
     return _mseq()
 
 
+def acorr(np.ndarray[np.float64_t, ndim=1, mode="c"] x not None, order):
+    cdef np.ndarray[np.float64_t, ndim = 1, mode = "c"] r
+    r = np.zeros(order + 1)
+    _acorr(&x[0], len(x), &r[0], order)
+    return r
+
+
 ### Adaptive mel-generalized cepstrum analysis ###
 
 def acep(x, np.ndarray[np.float64_t, ndim=1, mode="c"] c not None,

pysptk/sptk.py

@@ -10,6 +10,7 @@
     gexp
     glog
     mseq
+    acorr
 
 Adaptive cepstrum analysis
 --------------------------
@@ -262,6 +263,32 @@ def mseq():
     return _sptk.mseq()
 
 
+@apply_along_last_axis
+@automatic_type_conversion
+def acorr(x, order):
+    """Autocorrelation
+
+    Parameters
+    ----------
+    x : array
+        Input frame
+
+    order : int
+        Order of sequence
+
+    Returns
+    -------
+    r : array
+        Autocorrelation
+
+    See Also
+    --------
+    pysptk.sptk.levdur
+    pysptk.sptk.c2acr
+    """
+    return _sptk.acorr(x, order)
+
+
 ### Adaptive mel-generalized cepstrum analysis ###
 
 def acep(x, c, lambda_coef=0.98, step=0.1, tau=0.9, pd=4, eps=1.0e-6):
@@ -819,7 +846,7 @@ def fftcep(logsp,
 
 @apply_along_last_axis
 @automatic_type_conversion
-def lpc(windowed, order=25, min_det=1.0e-6):
+def lpc(windowed, order=25, min_det=1.0e-6, use_scipy=True):
     """Linear prediction analysis
 
     Parameters
@@ -834,6 +861,10 @@ def lpc(windowed, order=25, min_det=1.0e-6):
         Mimimum value of the determinant of normal matrix.
         Default is 1.0e-6.
 
+    use_scipy : bool
+        Use scipy's solve_toeplitz implementation or not. Default is True.
+        This would be more numerically stable than SPTK.
+
     Returns
     -------
     a : array, shape (``order + 1``)
@@ -858,7 +889,11 @@ def lpc(windowed, order=25, min_det=1.0e-6):
     pysptk.sptk.lspdf
 
     """
-    return _sptk.lpc(windowed, order, min_det)
+    if use_scipy:
+        r = _sptk.acorr(windowed, order)
+        return levdur(r, use_scipy=True)
+    else:
+        return _sptk.lpc(windowed, order, min_det)
 
 
 ### MFCC ###
@@ -2714,7 +2749,7 @@ def lspcheck(lsp):
 
 @apply_along_last_axis
 @automatic_type_conversion
-def levdur(r, eps=0.0):
+def levdur(r, eps=0.0, use_scipy=True):
     """Solve an Autocorrelation Normal Equation Using Levinson-Durbin Method
 
     Parameters
@@ -2725,6 +2760,10 @@ def levdur(r, eps=0.0):
     eps : float
         Singular check (eps(if -1., 0.0 is assumed))
 
+    use_scipy : bool
+        Use scipy's solve_toeplitz implementation or not. Default is True.
+        This would be more numerically stable than SPTK.
+
     Returns
     -------
     a : array
@@ -2739,4 +2778,14 @@ def levdur(r, eps=0.0):
     --------
     pysptk.sptk.c2acr
     """
-    return _sptk.levdur(r, eps)
+    if use_scipy:
+        from scipy.linalg import solve_toeplitz
+        a = np.empty_like(r)
+        # Rx = r where R is a toeplitz matrix
+        a[1:] = -solve_toeplitz(r[:-1], r[1:])
+        # http://www.seas.ucla.edu/~ingrid/ee213a/speech/vlad_present.pdf
+        G = np.sqrt(r[0] + (a[1:] * r[1:]).sum())
+        a[0] = G
+        return a
+    else:
+        return _sptk.levdur(r, eps)

setup.py

@@ -127,6 +127,7 @@ def run(self):
     cmdclass=cmdclass,
     install_requires=[
         'numpy >= 1.8.0',
+        'scipy',
         'six',
         'decorator'
     ],

tests/test_mgcep.py

@@ -220,22 +220,27 @@ def test_lpc():
     x = windowed_dummy_data(1024)
 
     def __test(order):
-        a = pysptk.lpc(x, order)
-        assert np.all(np.isfinite(a))
+        a1 = pysptk.lpc(x, order, use_scipy=False)
+        a2 = pysptk.lpc(x, order, use_scipy=True)
+        a3 = pysptk.levdur(pysptk.acorr(x, order))
 
-    for order in [15, 20, 25]:
+        assert np.all(np.isfinite(a1))
+        assert np.allclose(a1, a2)
+        assert np.allclose(a1, a3)
+
+    for order in [15, 20, 25, 30, 40]:
         yield __test, order
 
 
 def test_lpc_invalid_args():
     x = windowed_dummy_data(1024)
 
     def __test_min_det(min_det):
-        pysptk.lpc(x, min_det=min_det)
+        pysptk.lpc(x, min_det=min_det, use_scipy=False)
 
     yield raises(ValueError)(__test_min_det), -1.0
 
 
 @raises(RuntimeError)
 def test_lpc_failure():
-    pysptk.lpc(np.zeros(256), 40)
+    pysptk.lpc(np.zeros(256), 40, use_scipy=False)

tests/test_synthesis.py

@@ -240,7 +240,10 @@ def __test_synthesis_levdur(filt):
         windowed = __dummy_windowed_frames(
             source, frame_len=512, hopsize=hopsize)
         c = pysptk.mcep(windowed, filt.order)
-        lpc = pysptk.levdur(pysptk.c2acr(c))
+        a = pysptk.c2acr(c)
+        lpc = pysptk.levdur(a)
+        lpc2 = pysptk.levdur(a, use_scipy=True)
+        assert np.allclose(lpc, lpc2)
 
         # make sure lpc has loggain
         lpc[:, 0] = np.log(lpc[:, 0])