【问题标题】:Finding first samples greater than a threshold value efficiently in Python (and MATLAB comparison)在 Python 中有效地查找大于阈值的第一个样本(和 MATLAB 比较)
【发布时间】:2014-07-14 22:52:36
【问题描述】:

我不想找到列表或数组中大于特定threshold 的所有样本/数据点,我只想找到signal 大于threshold 的第一个样本。信号可能会多次跨越阈值。例如,如果我有一个示例信号:

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

还有一个threshold = 2,然后

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

会给我signal 中大于threshold 的所有值。 但是,每当信号变得大于阈值时,我只想获得第一个样本。因此,我正在浏览整个列表并进行布尔比较,例如

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 倍?

提前感谢您的 cmets!

【问题讨论】:

    标签: python matlab numpy threshold


    【解决方案1】:

    其他答案给你第一个 True 的位置,如果你想要一个标记第一个 True 的 bool 数组,你可以做得更快:

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

    【讨论】:

    • 这个解决方案真的比上面的其他两个选项还要快。谢谢。
    • 它确实更快 - 而且它的扩展性也更好。我怀疑在大型阵列上缓存未命中的情况更少。不错的解决方案!
    【解决方案2】:

    基于加快速度的最佳方法是选择最佳算法这一概念,您可以使用简单的边缘检测器巧妙地做到这一点:

    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)中的事情变得更快——在我的机器上大约需要 12 毫秒才能完成需要 4.5 秒的工作。

    编辑:正如@eickenberg 所指出的,卷积可以用numpy.diff(thresholded_data) 代替,这在概念上有点简单,但在这种情况下,索引将减1,所以记得把它们加回去,并且还将thresholded_data 转换为带有thresholded_data.astype(int) 的整数数组。两种方法之间没有明显的速度差异。

    【讨论】:

    • 这看起来像是要走的路。我建议使用numpy.diff 以清晰和可能的速度(建议编辑)。
    • 使用numpy.diff实际上会稍微慢一些(尽管可能还不够担心)(而且它没有被我拒绝!)。
    • @Henry Gomersall:感谢您对边缘检测器的建议。这也适用于我,并且在我的 PC 上产生的速度与 mskimm 的解决方案相似。
    • @xaneon 确保您将数组初始化为 numpy 数组 - 如果您需要多次执行此操作,列表之间的转换会减慢您的速度。
    • 感谢您的测试! np.diff 没有更快的事实可能是由于样本量。如果它总是较慢,那么这真的很有趣,因为这意味着 np.convolution 即使在特定情况下(例如这种情况)也确实有效地完成了它的(稍微更一般的)工作。如果使用np.diff,请注意在thresholded_data.astype(int) 上使用它,否则所有差异将仅计算为True-1 也是。
    【解决方案3】:

    This post 解释了为什么您的代码比 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 世界。

    【讨论】:

    • 感谢您提供有关 JIT 加速器说明的链接。很高兴知道何时使用循环。在我的机器上测试了你上面的版本,它运行良好并且只用了0.02 秒。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2023-03-23
    • 2018-03-14
    • 1970-01-01
    • 2011-04-14
    • 1970-01-01
    • 2016-11-17
    相关资源
    最近更新 更多