python - 長い時系列のPythonをフィルタリングする最も効率的な方法

Question

私は、神経活動、すなわち電圧を記録することから生じる大きな時系列、たとえば1e10を持っています。さらに分析を行う前に、300 Hz〜7000Hzのデータをバンドパスフィルター処理したいと思います。以下に、私が設計したバターワースフィルターのコードを投稿します。このフィルターを高速化するにはどうすればよいですか？実行に時間がかかりすぎます。

更新：サンプルデータ

の各行にあるようなデータへのリンクは次のとおりですdata。

フォーマットに関しては、各行は異なる記録ソースを表し、各列は特定の時点を表します。データは20000Hzでサンプリングされます。

 def butter_bandpass(lowcut,highcut,fs,order=8):
     nyq = 0.5*fs
     low = lowcut/nyq
     high = highcut/nyq

     b,a = butter(order, [low, high], btype='band')
     return b,a

 def butter_bandpass_filter(data,lowcut,highcut,fs,order=8):
     b,a = butter_bandpass(lowcut,highcut,fs,order=order)
     return lfilter(b,a,data) 

 print 'Band-pass filtering data'
 data = map(lambda channel: butter_bandpass_filter(channel,300,7000,20000),data)

アップデート

ほとんどのNumPyと同様に、SciPy関数lfilterは多次元入力を受け取ることができるため、map不要なオーバーヘッドが発生します。つまり、書き直すことができます

 data = map(lambda channel:butter_bandpass_filter(channel,300,7000,20000),data)

なので

 data = butter_bandpass_filter(data,300,7000,20000)

デフォルトlfilterでは、最後の非シングルトン軸で動作します。2D行列の場合、これは関数が各行に適用されることを意味します。これはまさに私が必要としていることです。悲しいことに、それでもまだ遅すぎます。

Answer 1

まず、あなたのデータ サンプルは独自の形式ですよね? Python 用の biosig ツールボックスを使用しても、この形式を読み取ることはできません。間違っているかもしれませんが、読めませんでした。

したがって、レスラー発振器から生成された人工データに基づいて回答します。これは混沌とした 3D オシレーターであり、非線形時系列解析の分野でよく使用されます。

import numpy as np
from scipy.signal import butter, lfilter

##############################################################
# For generating sample-data
##############################################################
from scipy.integrate import odeint

def roessler_ode(y,t,omega=1,a=0.165,b=0.2,c=10):
    dy = np.zeros((3))
    dy[0] = -1.0*(omega*y[1] + y[2]) #+ e1*(y[3]-y[0])
    dy[1] = omega * y[0] + a * y[1]
    dy[2] = b + y[2] * (y[0] - c)
    return dy

class Roessler(object):
    """A single coupled Roessler oscillators"""
    def __init__(self, y=None, omega=1.0, a=0.165,b=0.2,c=10):
        self.omega = omega
        self.a = a
        self.b = b
        self.c = c
        if y==None:
            self.y = np.random.random((3))+0.5
        else:
            self.y = y

    def ode(self,y,t):
        dy = roessler_ode(y[:],t,self.omega,self.a,self.b,self.c)
        return dy

    def integrate(self,ts):
        rv = odeint(self.ode,self.y,ts)
        self.y = rv[-1,:]
        return rv
###############################################################

次に、関数定義について説明します。

def butter_bandpass(lowcut,highcut,fs,order=8):
    nyq = 0.5*fs
    low = lowcut/nyq
    high = highcut/nyq

    b,a = butter(order, [low, high], btype='band')
    return b,a

def butter_bandpass_filter(data,lowcut,highcut,fs,order=8):
    b,a = butter_bandpass(lowcut,highcut,fs,order=order)
    return lfilter(b,a,data) 

私はそれらを変更せずに残しました。オシレーターでいくつかのデータを生成しますが、その 3 番目のコンポーネントしか取得しません。何かを除外するためにガウス ノイズを追加します。

# generate sample data
data = Roessler().integrate(np.arange(0,1000,0.1))[:,2]
data += np.random.normal(size=data.shape)

それでは、速度の問題に移りましょう。timeit-module を使用して実行時間を確認します。これらのステートメントは、フィルタリングを 100 回実行し、全体の時間を測定します。この測定は次数 2 と次数 8 で行います (はい、よりシャープなスペクトル エッジが必要です。わかっていますが、待ってください)。

# time execution
from timeit import timeit
time_order8 = timeit("butter_bandpass_filter(data,300,2000,20000,8)", "from __main__ import butter_bandpass_filter, butter_bandpass, data", number=100)
time_order2 = timeit("butter_bandpass_filter(data,300,2000,20000,2)", "from __main__ import butter_bandpass_filter, butter_bandpass, data", number=100)
print "For order 8: %.2f seconds" % time_order8
print "For order 2: %.2f seconds" % time_order2

これら 2 つの print ステートメントの出力は次のとおりです。

For order 8: 11.70 seconds
For order 2: 0.54 seconds

これは 20 倍になります。バターワース フィルターに次数 8 を使用することは、まったくお勧めできません。これが理にかなっている状況は考えられません。このようなフィルターを使用するときに発生するその他の問題を証明するために、これらのフィルターの効果を見てみましょう。次数 8 で 1 回、次数 2 で 1 回、バンドパス フィルタリングをデータに適用します。

data_bp8 = butter_bandpass_filter(data,300,2000,20000,8)
data_bp2 = butter_bandpass_filter(data,300,2000,20000,2)

それでは、いくつかのプロットを実行してみましょう。まず、単純な線 (x 軸は気にしませんでした)

# plot signals
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(data, label="raw")
plt.plot(data_bp8, label="order 8")
plt.plot(data_bp2, label="order 2")
plt.legend()

これにより、次のことがわかります。

緑の線はどうしたの？変ですね。その理由は、次数 8 のバターワース フィルターがかなり不安定なものになるためです。共鳴災害/大惨事について聞いたことがありますか? これは、それがどのように見えるかです。

これらの信号のパワー スペクトル密度は、次のようにプロットできます。

# plot power spectral densities
plt.figure(2)
plt.psd(data, Fs=200000, label="raw")
plt.psd(data_bp8, Fs=20000, label="order 8")
plt.psd(data_bp2, Fs=20000, label="order 2")
plt.legend()

plt.show()

ここでは、緑の線がより鋭いエッジを持っていることがわかりますが、価格はいくらですか? 人工的なピーク 300Hz。信号は完全に歪んでいます。

それで、あなたは何をしますか？

• 次数 8 のバターワース フィルターを使用しない
• 十分な場合は、より低い次数を使用してください。
• そうでない場合は、Parks-McGllellan または Remez-Exchange-Algorithms を使用して FIR フィルターを作成します。たとえば、scipy.signal.remez があります。

もう 1 つのヒント: 信号の位相が気になる場合は、時間の前後に確実にフィルタリングする必要があります。lfilterそれ以外の場合は位相をシフトします。一般に と呼ばれるこのようなアルゴリズムの実装は、filtfilt私の github リポジトリ にあります。

もう 1 つのプログラミングのヒント:

パススルー パラメーターの状況がある場合( の 4 つのパラメーターがbutter_bandpass_filterにのみ渡される場合、とbutter_bandpassを使用できます。*args**kwargs

def butter_bandpass(lowcut,highcut,fs,order=8):
    nyq = 0.5*fs
    low = lowcut/nyq
    high = highcut/nyq

    b,a = butter(order, [low, high], btype='band')
    return b,a

def butter_bandpass_filter(data, *args, **kwargs):
    b,a = butter_bandpass(*args, **kwargs)
    return lfilter(b,a,data) 

これにより、コードの冗長性が減り、コードのエラーが発生しにくくなります。

最後に、簡単にコピーして貼り付けて試してみるための完全なスクリプトを次に示します。

import numpy as np
from scipy.signal import butter, lfilter

##############################################################
# For generating sample-data
##############################################################
from scipy.integrate import odeint

def roessler_ode(y,t,omega=1,a=0.165,b=0.2,c=10):
    dy = np.zeros((3))
    dy[0] = -1.0*(omega*y[1] + y[2]) #+ e1*(y[3]-y[0])
    dy[1] = omega * y[0] + a * y[1]
    dy[2] = b + y[2] * (y[0] - c)
    return dy

class Roessler(object):
    """A single coupled Roessler oscillators"""
    def __init__(self, y=None, omega=1.0, a=0.165,b=0.2,c=10):
        self.omega = omega
        self.a = a
        self.b = b
        self.c = c
        if y==None:
            self.y = np.random.random((3))+0.5
        else:
            self.y = y

    def ode(self,y,t):
        dy = roessler_ode(y[:],t,self.omega,self.a,self.b,self.c)
        return dy

    def integrate(self,ts):
        rv = odeint(self.ode,self.y,ts)
        self.y = rv[-1,:]
        return rv
###############################################################


def butter_bandpass(lowcut,highcut,fs,order=8):
    nyq = 0.5*fs
    low = lowcut/nyq
    high = highcut/nyq

    b,a = butter(order, [low, high], btype='band')
    return b,a

def butter_bandpass_filter(data,lowcut,highcut,fs,order=8):
    b,a = butter_bandpass(lowcut,highcut,fs,order=order)
    return lfilter(b,a,data) 

# generate sample data
data = Roessler().integrate(np.arange(0,1000,0.1))[:,2]
data += np.random.normal(size=data.shape)

# time execution
from timeit import timeit
time_order8 = timeit("butter_bandpass_filter(data,300,2000,20000,8)", "from __main__ import butter_bandpass_filter, butter_bandpass, data", number=100)
time_order2 = timeit("butter_bandpass_filter(data,300,2000,20000,2)", "from __main__ import butter_bandpass_filter, butter_bandpass, data", number=100)
print "For order 8: %.2f seconds" % time_order8
print "For order 2: %.2f seconds" % time_order2

data_bp8 = butter_bandpass_filter(data,300,2000,20000,8)
data_bp2 = butter_bandpass_filter(data,300,2000,20000,2)

# plot signals
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(data, label="raw")
plt.plot(data_bp8, label="order 8")
plt.plot(data_bp2, label="order 2")
plt.legend()

# plot power spectral densities
plt.figure(2)
plt.psd(data, Fs=200000, label="raw")
plt.psd(data_bp8, Fs=20000, label="order 8")
plt.psd(data_bp2, Fs=20000, label="order 2")
plt.legend()

plt.show()