替换 scikits.talkbox 中的 Levinson 实现

Question

模块 scikits.talkbox 包含一些用于 Levinson-Durbin 递归的 C 代码。不幸的是，这段代码在最新版本的 Python 中不起作用，我想用纯 Python 实现替换它。 （速度不是问题，只要它有效。）

损坏的 C 函数的文档字符串如下：

Levinson-Durbin recursion, to efficiently solve symmetric linear systems
with toeplitz structure.

Parameters
----------
r : array-like
    input array to invert (since the matrix is symmetric Toeplitz, the
    corresponding pxp matrix is defined by p items only). Generally the
    autocorrelation of the signal for linear prediction coefficients
    estimation. The first item must be a non zero real, and corresponds
    to the autocorelation at lag 0 for linear prediction.
order : int
    order of the recursion. For order p, you will get p+1 coefficients.
axis : int, optional
    axis over which the algorithm is applied. -1 by default.

Returns
-------
a : array-like
    the solution of the inversion (see notes).
e : array-like
    the prediction error.
k : array-like
    reflection coefficients.

Notes
-----
Levinson is a well-known algorithm to solve the Hermitian toeplitz
equation: ::

                   _          _
    -R[1] = R[0]   R[1]   ... R[p-1]    a[1]
     :      :      :          :         :
     :      :      :          _      *  :
    -R[p] = R[p-1] R[p-2] ... R[0]      a[p]

with respect to a. Using the special symmetry in the matrix, the inversion
can be done in O(p^2) instead of O(p^3).

Only double argument are supported: float and long double are internally
converted to double, and complex input are not supported at all.

我看到有一个函数scipy.linalg.solve_toeplitz 看起来像我想要的。但是，它无法指定顺序，而是将数组的元组作为输入。

我承认我在这里有点超出我的深度，并且不完全理解这段代码应该做什么。有没有一种简单的方法可以用 Numpy 的 solve_toeplitz 替换对损坏的 C 函数的调用？

Answer 1

scikits.talkbox 包还包括一个名为 py_lpc 的模块，其中包含 LPC 估计的纯 Python 实现。适应这一点并不难。

def lpc2(signal, order):
    """Compute the Linear Prediction Coefficients.

    Return the order + 1 LPC coefficients for the signal. c = lpc(x, k) will
    find the k+1 coefficients of a k order linear filter:

      xp[n] = -c[1] * x[n-2] - ... - c[k-1] * x[n-k-1]

    Such as the sum of the squared-error e[i] = xp[i] - x[i] is minimized.

    Parameters
    ----------
    signal: array_like
        input signal
    order : int
        LPC order (the output will have order + 1 items)"""

    order = int(order)

    if signal.ndim > 1:
        raise ValueError("Array of rank > 1 not supported yet")
    if order > signal.size:
        raise ValueError("Input signal must have a lenght >= lpc order")

    if order > 0:
        p = order + 1
        r = np.zeros(p, signal.dtype)
        # Number of non zero values in autocorrelation one needs for p LPC
        # coefficients
        nx = np.min([p, signal.size])
        x = np.correlate(signal, signal, 'full')
        r[:nx] = x[signal.size-1:signal.size+order]
        phi = np.dot(sp.linalg.inv(sp.linalg.toeplitz(r[:-1])), -r[1:])
        return np.concatenate(([1.], phi)), None, None
    else:
        return np.ones(1, dtype = signal.dtype), None, None

这会在不使用任何特殊 Toeplitz 属性的情况下反转矩阵，使其速度明显变慢，但无需任何 C 代码即可工作。