scikits.talkboxでのLevinson実装の置き換え

Question

モジュールscikits.talkboxには、Levinson-Durbin再帰のためのCコードが含まれています。残念ながら、このコードはPythonの最近のバージョンでは動作しません。純粋なPython実装に置き換えたいと思います。 （速度は、それが作品として、問題ではありません。）

壊れたC関数のドキュメンテーション文字列を読み取ります。私は私が欲しいもののように見える機能scipy.linalg.solve_toeplitzがあることを見

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. 

。ただし、順序を指定する方法はなく、配列のタプルを入力として受け取ります。

私はここで少し深く考えていますが、このコードが何をすべきかを十分に理解していません。 Numpyの壊れたC関数の呼び出しをNumpyのsolve_toeplitzに置き換える簡単な方法はありますか？

Answer 1

scikits.talkboxパッケージには、py_lpcというモジュールが含まれています。このモジュールには、純粋なPythonのLPC推定の実装が含まれています。これを適応させることはあまり難しくありません。

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コードなしで動作します。