过滤Python长时间序列的最有效方法

Question

我有一个很大的时间序列,比如1e10,这是由记录神经活动(即电压)引起的.在进行进一步分析之前,我想对300 Hz和7000 Hz之间的数据进行带通滤波.下面,我发布了我设计的Butterworth滤波器的代码.如何使此过滤器更快？运行需要很长时间.

更新:示例数据

这是一个链接到数据,如将在每一行data.

至于格式化,每行代表不同的记录源,每列代表一个时间点.数据以20,000Hz采样.

 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工具箱,也无法读取此格式.也许我错了,但我没有成功阅读它.

因此,我的答案基于Rössler振荡器生成的人工数据.它是一种混沌的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

这两个打印语句的输出是:

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

这个数字是20!使用命令8用于butterworth过滤器绝对不是一个好主意.我想不出任何有意义的情况.为了证明使用这种滤波器时出现的其他问题,让我们看看这些滤波器的效果.我们对数据应用带通滤波,一次是订单8,一次是订单2:

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阶的butterworth滤波器变得相当不稳定.有没有听说过共振灾难/灾难？这就是它的样子.

这些信号的功率谱密度可以绘制为:

# 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()

在这里,您看到绿线边缘更锐利,但价格是多少？人工峰值约为 300赫兹.信号完全失真.

那你该怎么办？

• 切勿使用8阶的butterworth过滤器
• 如果足够,请使用较低的顺序.
• 如果没有,请使用Parks-McGlellan或Remez-Exchange-Algorithms创建一些FIR滤波器.例如,有scipy.signal.remez.

另一个暗示:如果你关心信号的相位,你肯定应该及时过滤前后.lfilter否则会改变阶段.这种算法的实现,通常被称为filtfilt,可以在我的github存储库中找到.

还有一个编程提示:

如果你有传递参数的情况(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()