具有 numpy 1d 数组容差的阈值过滤器

Question

假设我有一个 1d numpy 数组，其中包含一些嘈杂的数据系列。

我想建立一个阈值来检查值什么时候是high，什么时候是low。 然而，因为数据是嘈杂的，只做是没有意义的

is_high = data > threshold

我试图设置此阈值的容差，就像许多控制系统（例如大多数供暖和空调系统）所做的那样。这个想法是，信号的状态只有在超过阈值加上容差时才会从低变为高。同样，只有在低于阈值减去容差时，信号才会从高变为低。换句话说：

def tolerance_filter(data, threshold, tolerance):
    currently_high = False  # start low
    signal_state = np.empty_like(data, dtype=np.bool)
    for i in range(data.size):
        # if we were high and we are getting too low, become low
        if currently_high and data[i] < (threshold-tolerance):
            currently_high = False
        # if we were low and are getting too high, become high
        elif not currently_high and data[i] > (threshold+tolerance):
            currently_high = True
        signal_state[i] = currently_high
    return signal_state

这个函数给出了我期望的输出。但是，我想知道是否有任何方法可以使用 numpy 或 scipy 的速度而不是原始 python for 循环来做到这一点。

有什么想法吗？ :)

---

更新：

感谢 Joe Kington 的评论指出了 滞后 术语，我找到了this other SO 问题。恐怕它很相似（重复？），还有 Bas Swinckels 的one nice working solution。

无论如何，我尝试实现 Joe Kington 建议的加速（不知道我是否做对了），并将他的解决方案、Fergal 和 Bas 的解决方案与我的幼稚方法进行了比较。结果如下（代码如下）：

Proposed function in my original question
10 loops, best of 3: 22.6 ms per loop

Proposed function by Fergal
1000 loops, best of 3: 995 µs per loop

Proposed function by Bas Swinckels in the hysteresis question
1000 loops, best of 3: 1.05 ms per loop

Proposed function by Joe Kington using Cython
Approximate time cost of compiling: 2.195411
1000 loops, best of 3: 1.35 ms per loop

答案中的所有方法都执行相似（尽管 Fergal 需要一些额外的步骤来获得布尔向量！）。有没有考虑在这里添加？ 此外，我很惊讶Cython 方法速度较慢（不过只是稍微慢了一点）。无论如何，我不得不承认，如果您不熟悉所有 numpy 函数，它可能是最快的编码方法......

这是我用来对不同选项进行基准测试的代码。非常欢迎审核和修订！ :P （Cython 代码在中间强制 SO 将所有代码保留在同一个可滚动块中。当然我把它放在不同的文件中）

# Naive approach from the original question
def tolerance_filter1(data, threshold, tolerance):
    currently_high = False  # start low
    signal_state = np.empty_like(data, dtype=np.bool)
    for i in range(data.size):
        # if we were high and we are getting too low, become low
        if currently_high and data[i] < (threshold-tolerance):
            currently_high = False
        # if we were low and are getting too high, become high
        elif not currently_high and data[i] > (threshold+tolerance):
            currently_high = True
        signal_state[i] = currently_high
    return signal_state
        
# Numpythonic approach suggested by Fergal
def tolerance_filter2(data, threshold, tolerance):
    a = np.zeros_like(data)
    a[ data < threshold-tolerance] = -1
    a[ data > threshold+tolerance] = +1
    wh = np.where(a != 0)[0]
    idx= np.diff( a[wh]) == 2
    #This variable indexes the values of data where data crosses
    #from below threshold-tol to above threshold+tol
    crossesAboveThreshold = wh[idx]
    return crossesAboveThreshold
    
# Approach suggested by Bas Swinckels and borrowed
# from the hysteresis question
def tolerance_filter3(data, threshold, tolerance, initial=False):
    hi = data >= threshold+tolerance
    lo_or_hi = (data <= threshold-tolerance) | hi
    ind = np.nonzero(lo_or_hi)[0]
    if not ind.size: # prevent index error if ind is empty
        return np.zeros_like(x, dtype=bool) | initial
    cnt = np.cumsum(lo_or_hi) # from 0 to len(x)
    return np.where(cnt, hi[ind[cnt-1]], initial)
    
#########################################################
## IN A DIFFERENT FILE (tolerance_filter_cython.pyx)
## So that StackOverflow shows a single scrollable code block :)

import numpy as np
import cython

@cython.boundscheck(False)
def tolerance_filter(data, float threshold, float tolerance):
    cdef bint currently_high = 0  # start low
    signal_state = np.empty_like(data, dtype=int)
    cdef double[:] data_view = data
    cdef long[:] signal_state_view = signal_state
    cdef int i = 0
    cdef int l = len(data)
    low = np.empty_like(data, dtype=bool)
    high = np.empty_like(data, dtype=bool)
    low = data < (threshold-tolerance)
    high = data > (threshold+tolerance)
    
    for i in range(l):
        # if we were high and we are getting too low, become low
        if currently_high and low[i]:
            currently_high = False
        # if we were low and are getting too high, become high
        elif not currently_high and high[i]:
            currently_high = True
        signal_state_view[i] = currently_high
    return signal_state

##################################################################
# BACK TO THE PYTHON FILE

import numpy as np
from time import clock
from datetime import datetime
from IPython import get_ipython
ipython = get_ipython()
time = np.arange(0,1000,0.01)
data = np.sin(time*3) + np.cos(time/7)*8 + np.random.normal(size=time.shape)*2
threshold, tolerance = 0, 4

print "Proposed function in my original question"
ipython.magic("timeit tolerance_filter1(data, threshold, tolerance)")

print "\nProposed function by Fergal"
ipython.magic("timeit tolerance_filter2(data, threshold, tolerance)")

print "\nProposed function by Bas Swinckels in the hysteresis question"
ipython.magic("timeit tolerance_filter3(data, threshold, tolerance)")

print "\nProposed function by Joe Kington using Cython"
start = datetime.now()
import pyximport; pyximport.install()
import tolerance_filter_cython
print "Approximate time cost of compiling: {}".format((datetime.now()-start).total_seconds())
tolerance_filter4 = tolerance_filter_cython.tolerance_filter
ipython.magic("timeit tolerance_filter4(data, threshold, tolerance)")

Answer 1

我认为看到cython extensions 如此简单和类似 Python 有时会令人惊讶。这是您的代码转换为 cython。它可以从 Python 调用，但应该会给你 C++ 的速度。

def tolerance_filter(data, float threshold, float tolerance):
    cdef bint currently_high = 0  # start low
    signal_state = np.empty_like(data, dtype=int)
    cdef float[:] data_view = data
    cdef int[:] signal_state_view = signa_state
    cdef int i = 0
    cdef int l = len(data)
    for i in range(l):
        # if we were high and we are getting too low, become low
        if currently_high and data[i] < (threshold-tolerance):
            currently_high = False
        # if we were low and are getting too high, become high
        elif not currently_high and data[i] > (threshold+tolerance):
            currently_high = True
        signal_state_view[i] = currently_high

有几点需要注意：

• 注意函数开头typed memory views的使用

• 故意使函数尽可能接近您的原始代码。但是，可以通过关闭范围检查（请参阅 Cython 文档）并在循环外部计算上限和下限有效阈值来加快速度。

Answer 2

我不确定这是否比您的解决方案更好，但它是否更 numpythonic。

a = np.zeros_like(data)
a[ data < threshold-tol] = -1
a[ data > threshold+tol] = +1
wh = np.where(a != 0)
idx= np.diff( a[wh]) == 2
#This variable indexes the values of data where data crosses
#from below threshold-tol to above threshold+tol
crossesAboveThreshold = wh[idx]