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

Question

我有一个很大的时间序列，比如 1e10，它是记录神经事件(即电压)的结果。在进行进一步分析之前，我想对 300 Hz 和 7000 Hz 之间的数据进行带通滤波。下面，我发布了我设计的 Butterworth 滤波器的代码。 如何让这个过滤器更快？运行时间太长。

更新:样本数据

Here是指向数据的链接，例如 data 的每一行。

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

 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 在最后一个非单例轴上运行。对于二维矩阵，这意味着函数应用于每一行，这正是我所需要的。遗憾的是，它仍然太慢了。

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) 

我没有改变它们。我用我的振荡器生成了一些数据，但我只取了它的第三个分量。我添加了一些高斯噪声以便过滤掉一些东西。

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

现在让我们来谈谈速度问题。我使用 timeit 模块来检查执行时间。这些语句执行过滤 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 绝对不是一个好主意。我想不出在任何情况下这都有意义。为了证明使用这种过滤器时出现的其他问题，让我们看看这些过滤器的效果。我们对数据应用带通滤波，一次是 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 阶巴特沃斯滤波器变得相当不稳定。听说过resonance disaster / catastrophe ？这是它的样子。

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

# 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 阶巴特沃斯滤波器
• 如果足够，使用低阶。
• 如果没有，请使用 Parks-McGellan 或 Remez-Exchange-Algorithms 创建一些 FIR 滤波器。例如，有 scipy.signal.remez。

另一个提示:如果您关心信号的相位，您绝对应该及时向前和向后过滤。 lfilter 确实改变了相位。这种算法的实现，通常称为 filtfilt，可以在 my github repository 找到。 .

另一个编程提示:

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