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

Question

我有一个很大的时间序列，比如 1e10，它是由记录神经活动（即电压）产生的。在进行进一步分析之前，我想对 300 Hz 到 7000 Hz 之间的数据进行带通滤波。下面，我贴出我设计的巴特沃斯滤波器的代码。如何使此过滤器更快？运行时间太长。

更新：示例数据

这是一个指向数据的链接，例如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，在最后一个非单轴上运行。对于 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) 

我让它们保持不变。我用我的振荡器生成了一些数据，但我只取了它的第三部分。我添加了一些高斯噪声以便过滤掉一些东西。

# 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 阶绝对不是一个好主意。我想不出任何情况下这是有意义的。为了证明使用这种过滤器时出现的其他问题，让我们看看这些过滤器的效果。我们对我们的数据应用带通滤波，一次是 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 阶巴特沃斯滤波器变得相当不稳定。听说过共振灾难/灾难吗？这就是它的样子。

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

# 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-McGlellan 或 Remez-Exchange-Algorithms 创建一些 FIR 滤波器。例如，有 scipy.signal.remez。

另一个提示：如果您关心信号的相位，则绝对应该及时向前和向后过滤。lfilter否则会改变阶段。这种算法的实现，通常称为filtfilt，可以在我的 github 存储库中找到。

另一个编程提示：

如果有传参的情况（四个参数butter_bandpass_filter只传给butter_bandpass，可以使用*argsand **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()