在Python中高效计算事件/刺激触发的平均值

Question

我想要高效地计算事件或刺激触发的平均值。假设我有一个 signal，比如说：

signal = [random.random() for i in xrange(0, 1000)]

这个信号有 n_signal 个数据点。

n_signal = len(signal) 

我知道这个信号的采样频率是：

Fs = 25000 # Hz

在这种情况下，我知道信号的总时间是：

T_sec = n_signal / float(Fs)  

在特定的时间点，会发生某些事件，比如：

t_events = [0.01, 0.017, 0.018, 0.022, 0.034, 0.0345, 0.03456]

现在我想找到这些事件发生之前的一段时间的信号，比如：

t_bef = 0.001

以及这些事件发生之后的一段时间，比如：

t_aft = 0.002

一旦我得到了所有这些信号片段，我想对它们进行平均。在过去，我会创建信号的时间向量：

t_signal = numpy.linspace(0, T_sec, n_signal)

然后在 t_signal 中查找所有 t_events 的索引，比如使用 numpy.searchsorted (链接)。

因为我知道信号的采样率，所以可以更快地完成这些操作，比如：

indices = [int(i * Fs) for i in t_events]

这样可以节省 t_signal 的内存，而且我不需要遍历整个信号来找到我的索引。

接下来，我会确定 t_bef 和 t_aft 对应的数据样本数量：

nsamples_t_bef = int(t_bef * Fs)
nsamples_t_aft = int(t_aft * Fs)

然后我会把信号片段保存在一个 list 中：

signal_chunks = list()
for i in xrange(0, len(t_events)):
    signal_chunks.append(signal[indices[i] - nsamples_t_bef : indices[i] + nsamples_t_aft])

最后，我会对这些片段进行平均：

event_triggered_average = numpy.mean(signal_chunks, axis = 0)

如果我对时间向量感兴趣，我会用以下方法计算：

t_event_triggered_average = numpy.linspace(-t_signal[nsamples_t_bef], t_signal[nsamples_t_aft], nsamples_t_bef + nsamples_t_aft)

现在我的问题是：有没有更高效的计算方法？如果我有一个数据点很多、事件也很多的信号，这样的计算可能会花费一些时间。使用 list 来保存这些片段是最好的数据结构吗？你知道如何更快地获取这些数据片段吗？也许可以使用缓冲区？感谢你们的评论和建议。

最小工作示例

import numpy
import random

random.seed(0)
signal = [random.random() for i in xrange(0, 1000)]

# sampling rate
Fs = 25000 # Hz

# total time of the signal
n_signal = len(signal) 
T_sec = n_signal / float(Fs)  

# time of events of interest
t_events = [0.01, 0.017, 0.018, 0.022, 0.034, 0.0345, 0.03456]

# and their corresponding indices
indices = [int(i * Fs) for i in t_events]

# define the time window of interest around each event
t_bef = 0.001
t_aft = 0.002

# and the corresponding index offset
nsamples_t_bef = int(t_bef * Fs)
nsamples_t_aft = int(t_aft * Fs)

# vector of signal times
t_signal = numpy.linspace(0, T_sec, n_signal)

signal_chunks = list()
for i in xrange(0, len(t_events)):
    signal_chunks.append(signal[indices[i] - nsamples_t_bef : indices[i] + nsamples_t_aft])

# average signal value across chunks
event_triggered_average = numpy.mean(signal_chunks, axis = 0)

# not sure what's going on here
t_event_triggered_average = numpy.linspace(-t_signal[nsamples_t_bef],
                                            t_signal[nsamples_t_aft],
                                            nsamples_t_bef + nsamples_t_aft)

Answer 1

因为你的信号是在一个规则的网格上定义的，所以你可以通过一些简单的数学运算来找到你需要的所有样本的索引。然后，你可以通过一次索引操作来构建包含这些样本的数组。

import numpy as np

# Making some test data
n_signal = 1000
signal = np.random.rand(n_signal)
Fs = 25000 # Hz
t_events = np.array([0.01, 0.017, 0.018, 0.022, 0.034, 0.0345, 0.03456])

# Preferences
t_bef = 0.001
t_aft = 0.002

# The number of samples in a chunk
nsamples = int((t_bef+t_aft) * Fs)

# Create a vector from 0 up to nsamples
sample_idx = np.arange(nsamples)

# Calculate the index of the first sample for each chunk
# Require integers, because it will be used for indexing
start_idx = ((t_events - t_bef) * Fs).astype(int)

# Use broadcasting to create an array with indices
# Each row contains consecutive indices for each chunk
idx = start_idx[:, None] + sample_idx[None, :]

# Get all the chunks using fancy indexing
signal_chunks = signal[idx]

# Calculate the average like you did earlier
event_triggered_average = signal_chunks.mean(axis=0)

需要注意的是，包含 .astype(int) 的那一行并不是四舍五入到最近的整数，而是向零的方向取整。