在Python中高效地查找首个超过阈值的样本（并与MATLAB比较）

Question

相较于在列表或数组中找出所有大于特定阈值的样本/数据点，我希望仅找出第一次信号超过阈值时的样本。信号可能会多次穿过阈值。例如，如果我有一个示例信号：

signal = [1, 2, 3, 4, 4, 3, 2, 1, 0, 3, 2, 1, 0, 0, 1, 1, 4, 8, 7, 6, 5, 0]

如果一个阈值为2，那么

signal = numpy.array(signal)
is_bigger_than_threshold = signal > threshold

会给我所有大于阈值的signal中的值。然而，我只想在信号第一次变大于阈值时获取前几个样本。因此，我要遍历整个列表，并进行布尔比较。

first_bigger_than_threshold = list()
first_bigger_than_threshold.append(False)
for i in xrange(1, len(is_bigger_than_threshold)):
    if(is_bigger_than_threshold[i] == False):
        val = False
    elif(is_bigger_than_threshold[i]):
        if(is_bigger_than_threshold[i - 1] == False):
            val = True
        elif(is_bigger_than_threshold[i - 1] == True):
            val = False
    first_bigger_than_threshold.append(val)

这给了我我想要的结果，即

[False, False, True, False, False, False, False, False, False, True, False, False, False,   
False, False, False, True, False, False, False, False, False]

在MATLAB中，我会类似地做。

for i = 2 : numel(is_bigger_than_threshold)
    if(is_bigger_than_threshold(i) == 0)
        val = 0;
    elseif(is_bigger_than_threshold(i))
        if(is_bigger_than_threshold(i - 1) == 0)
            val = 1;
        elseif(is_bigger_than_threshold(i - 1) == 1)
            val = 0;
        end
    end
    first_bigger_than_threshold(i) = val;
end % for

有没有更高效（更快）的方法来执行这个计算？

如果我在Python中生成数据，例如：

signal = [round(random.random() * 10) for i in xrange(0, 1000000)]

"计时，这些值的计算花了4.45秒。如果我在MATLAB中生成数据。"

signal = round(rand(1, 1000000) * 10);

执行该程序仅需0.92秒。为什么MATLAB执行此任务比Python快近5倍？感谢您提前的评论！

- xaneon

Answer 1

其他答案给出了第一个True的位置，如果你想要一个标记第一个True的布尔数组，可以更快地完成：

import numpy as np

signal = np.random.rand(1000000)
th = signal > 0.5
th[1:][th[:-1] & th[1:]] = False

- HYRY

---

这个解决方案比上面的两个选项都要快得多。谢谢。 - xaneon

确实更快 - 而且它也更具可扩展性。我怀疑在大型数组上会有较少的缓存未命中。不错的解决方案！ - Henry Gomersall

Answer 2

这篇文章解释了为什么你的代码比Matlab慢。

尝试使用这段代码

import numpy as np

threshold = 2
signal = np.array([1, 2, 3, 4, 4, 3, 2, 1, 0, 3, 2, 1, 0, 0, 1, 1, 4, 8, 7, 6, 5, 0])

indices_bigger_than_threshold = np.where(signal > threshold)[0] # get item
print indices_bigger_than_threshold
# [ 2  3  4  5  9 16 17 18 19 20]
non_consecutive = np.where(np.diff(indices_bigger_than_threshold) != 1)[0]+1 # +1 for selecting the next
print non_consecutive
# [4 5]
first_bigger_than_threshold1 = np.zeros_like(signal, dtype=np.bool)
first_bigger_than_threshold1[indices_bigger_than_threshold[0]] = True # retain the first
first_bigger_than_threshold1[indices_bigger_than_threshold[non_consecutive]] = True

np.where 返回满足条件的索引。

策略是获取高于threshold的索引并移除连续的索引。

顺便说一下，欢迎来到Python/Numpy世界。

- emesday

---

感谢您提供有关JIT加速器的解释链接。在使用循环时，这是一个好的知识点。我测试了您上面的版本，它可以正常工作，并且在我的机器上只需要0.02秒。 - xaneon

Answer 3

基于一个想法，即加速的最佳方式是选择最佳算法，你可以通过简单的边缘检测器来实现这一点：

import numpy

signal = numpy.array([1, 2, 3, 4, 4, 3, 2, 1, 0, 3, 2, 1, 0, 0, 1, 1, 4, 8, 7, 6, 5, 0])

thresholded_data = signal > threshold
threshold_edges = numpy.convolve([1, -1], thresholded_data, mode='same')

thresholded_edge_indices = numpy.where(threshold_edges==1)[0]

print(thresholded_edge_indices)

打印出[2 9 16]，这些索引对应于大于阈值的第一个条目。这将使Matlab和Python（使用Numpy）更快-在我的机器上约为12ms，而你需要4.5s才能完成。

编辑：正如@eickenberg所指出的，卷积可以替换为numpy.diff(thresholded_data)，这在概念上稍微简单一些，但在这种情况下，索引会偏移1，因此请记得将其添加回去，并将thresholded_data转换为整数数组，使用thresholded_data.astype(int)。两种方法之间没有明显的速度差异。

- Henry Gomersall

---

使用numpy.diff实际上会略微慢一些（虽然无疑不足以担忧）（并且这不是我拒绝的原因！）。 - Henry Gomersall

@Henry Gomersall：感谢您提供的边缘检测建议。这对我也很有效，并且在我的电脑上与mskimm的解决方案速度相似。 - xaneon

@xaneon请确保将数组初始化为numpy数组 - 如果您需要多次在列表之间进行转换，这样会减慢速度。 - Henry Gomersall

感谢您测试！np.diff不比速度更快可能是由于样本大小的原因。如果它始终比较慢，那么这真的很有趣，因为这意味着即使在特定情况下，如此，np.convolution也可以非常高效地完成其（稍微更一般的）工作。如果使用np.diff，请注意将其用于thresholded_data.astype(int)上，否则所有差异都将仅评估为“True”，包括"-1"。 - eickenberg

在重新检查后，时间差别并不实质性，并且实际上是相反的（diff用时11毫秒，而convolve用时12毫秒）。我会添加有关“astype”的注释。这并不完全出乎意料 - 长度为2的卷积并不算太艰难！ - Henry Gomersall