ラティス・フィルタのPythonによる実装

参考文献

「プログラム101付き 音声信号処理」2021/01/01
(著)川村 新
文献元

準備

console

pip install numpy

ソースコード

sample.py

# -*- coding: utf-8 -*-
import numpy as np

def standardization_(x):
    return (x-x.mean())/x.std()

def arc_standardization(x,original):
    return (x*original.std())+original.mean()

def Lattice_filter(x):
    a=standardization_(x)
    b=a/np.abs(a).max()
    s=b
    t=0
    MEM_SIZE=b.shape[1]
    t_out=0
    N=128#b.shape[1]*2                                            # フィルタ次数
    h=np.full((1,N+1),0.)
    y=np.full((x.shape[0],x.shape[1]),0.)
    
    f  =np.full((1,N+1),0.)
    b0 =np.full((1,N+1),0.)
    b1 =np.full((1,N+1),0.)          #// 前向き予測誤差，後ろ向き予測誤差
    rf =np.full((1,N+1),0.)
    rb =np.full((1,N+1),0.)                       #// 反射係数
    k  =np.full((1,N+1),0.)
    Eb0=np.full((1,N+1),0.)
    Eb1=np.full((1,N+1),0.)
    Ef =np.full((1,N+1),0.) #// 係数更新用変数
    
    #///////////////////////////////
    #//                           //
    #//       変数の初期設定      //
    #//                           //
    #///////////////////////////////
    _lambda=0.9                                            #// 忘却係数の設定
    add_len = b.shape[1]-1
    #//************************************************************************//
    
    #///////////////////////////////////
    #//                               //
    #//        メインループ           //
    #//                               //
    #///////////////////////////////////
    while True:                                               # メインループ
        
        if t_out > add_len:
            break                # ループ終了判定
        
        #//************************************************************************//
        
        #////////////////////////////////////////////////////
        #//                                                //
        #//              Signal Processing                 //
        #//                                                //
        #//  現在時刻tの入力 s[t] から出力 y[t] をつくる   //
        #//                                                //
        #//  ※ tは0からMEM_SIZE-1までをループ             //
        #//                                                //
        #////////////////////////////////////////////////////
        
        f[:,0]  = s[:,t]                                                #// ラティスフィルタへの入力信号
        b0[:,0] = s[:,t]
        Ef[:,0] = _lambda * Ef[:,0] + (s[:,t] * s[:,t])                 #// 入力信号の分散の推定
        Eb0[:,0]= _lambda * Ef[:,0] + (s[:,t] * s[:,t])
        for i in range(1,N,1):                                          #// 出力計算ループ
            k[:,i-1]   = _lambda  *   k[:,i-1] +  f[:,i-1] *  b1[:,i-1] #// 反射係数の分子
            if Eb1[:,i-1] != 0.0:
                rf[:,i] = -k[:,i-1] / Eb1[:,i-1]#// 前向き反射係数の更新
            if Ef[:,i-1] != 0.0:
                rb[:,i] = -k[:,i-1] /  Ef[:,i-1]#// 後向き反射係数の更新
            f  [:,i]   =  f[:,i-1] + rf[:,i] * b1[:,i-1]            #// 前向き予測誤差
            b0 [:,i]   = b1[:,i-1] + rb[:,i] *  f[:,i-1]            #// 後向き予測誤差
            Ef [:,i]   = _lambda *  Ef[:,i]   +  f[:,i] *  f[:,i]   #// 前向き反射係数の分母更新
            Eb0[:,i]   = _lambda * Eb0[:,i]   + b0[:,i] * b0[:,i]   #// 後向き反射係数の分母更新
            b1 [:,i-1] =  b0[:,i-1]                                 #// 後向き予測誤差の信号遅延
            Eb1[:,i-1] = Eb0[:,i-1]
        
        
        y[:,t]=s[:,t]-f[:,N]                            #// 観測信号から予測誤差を減算(予測値)
        
        #//************************************************************************//
        
        
        t=(t+1)%MEM_SIZE                                # 時刻 t の更新
        t_out+=1                                        # ループ終了時刻の計測
    y[:,:] = y[:,:]*np.abs(a).max()                 # 出力を整数化
    return arc_standardization(y,x)