如何找到带滞后的零交叉点？

Question

在numpy中，我想找出信号从低于某个阈值变为高于另一个阈值的点。这种情况通常出现在去抖动或者在噪声影响下准确检测零交叉等场景。

就像这样：

import numpy

# set up little test problem
N = 1000
values = numpy.sin(numpy.linspace(0, 20, N))
values += 0.4 * numpy.random.random(N) - 0.2
v_high = 0.3
v_low = -0.3

# find transitions from below v_low to above v_high    
transitions = numpy.zeros_like(values, dtype=numpy.bool)

state = "high"

for i in range(N):
    if values[i] > v_high:
        # previous state was low, this is a low-to-high transition
        if state == "low":
            transitions[i] = True
        state = "high"
    if values[i] < v_low:
        state = "low"

我希望能有一种方法来实现这个，而不需要显式地遍历整个数组。但是我想不出有什么办法，因为每个状态的值都依赖于前一个状态。有没有可能不使用循环来做到这一点呢？

Answer 1

这里有一个解决方案，它能达到和Bas Swinckels的回答完全相同的效果，而且在我看来更容易理解，性能也差不多（在Colab上测试时，处理1000万个元素大约需要0.4秒）：

def hyst(x, th_lo, th_hi, initial = False):
    outside_values = np.full(x.size, np.nan)
    outside_values[0] = initial
    outside_values[x < th_lo] = 0
    outside_values[x > th_hi] = 1
    outside_value_indexes = np.where(np.isnan(outside_values), 0, np.arange(x.size))
    np.maximum.accumulate(outside_value_indexes, out=outside_value_indexes)
    return outside_values[outside_value_indexes]

接下来，我们运行和其他回答相同的测试：

x = np.linspace(0,20, 1000)
y = np.sin(x)
h1 = hyst(y, -0.5, 0.5)
h2 = hyst(y, -0.5, 0.5, True)
plt.plot(x, y, x, -0.5 + h1, x, -0.5 + h2)
plt.legend(('input', 'output, start=0', 'output, start=1'))
plt.title('Thresholding with hysteresis')
plt.show()

这个方法的基本思路来源于一个关于前向填充的不同问题的回答：

1. 首先，我们创建一个叫做outside_values的数组，这个数组会把输入值映射成：如果小于某个阈值就变成0，如果大于阈值就变成1，其他情况则用NaN作为占位符。
2. 然后，我们创建一个叫做outside_values_indexes的数组，这个数组列出了所有outside_values的索引（也就是[0, 1, 2, 3, ...]）。那些映射到内部值（NaN）的索引会被替换成0。
3. 接着，我们把outside_values_indexes替换成它的累积最大值，也就是说，每个索引会被替换成它前面所有索引中的最大值。由于我们把内部值的索引设为0，所以它们不会影响最大值，而是使用所有前面外部值的最大索引（换句话说，就是最后一个看到的外部值的索引）。
4. 现在，outside_values_indexes已经把每个输入值映射到最后一个外部值的索引，我们可以用这个索引去查找outside_values，这样就完成了！

Answer 2

为了我的工作，我根据Bas Swinckels的回答做了一些修改，这样可以在使用标准阈值和反向阈值时检测到阈值穿越。

不过我对命名不太满意，可能现在应该叫th_hi2lo和th_lo2hi，而不是th_lo和th_hi？使用原来的值时，行为是一样的。

def hyst(x, th_lo, th_hi, initial = False):
    """
    x : Numpy Array
        Series to apply hysteresis to.
    th_lo : float or int
        Below this threshold the value of hyst will be False (0).
    th_hi : float or int
        Above this threshold the value of hyst will be True (1).
    """        

    if th_lo > th_hi: # If thresholds are reversed, x must be reversed as well
        x = x[::-1]
        th_lo, th_hi = th_hi, th_lo
        rev = True
    else:
        rev = False

    hi = x >= th_hi
    lo_or_hi = (x <= th_lo) | hi

    ind = np.nonzero(lo_or_hi)[0]  # Index für alle darunter oder darüber
    if not ind.size:  # prevent index error if ind is empty
        x_hyst = np.zeros_like(x, dtype=bool) | initial
    else:
        cnt = np.cumsum(lo_or_hi)  # from 0 to len(x)
        x_hyst = np.where(cnt, hi[ind[cnt-1]], initial)

    if rev:
        x_hyst = x_hyst[::-1]

    return x_hyst

接下来是对代码的测试，看看它的效果：

x = np.linspace(0,20, 1000)
y = np.sin(x)
h1 = hyst(y, -0.2, 0.2)
h2 = hyst(y, +0.5, -0.5)
plt.plot(x, y, x, -0.2 + h1*0.4, x, -0.5 + h2)
plt.legend(('input', 'output, classic, hyst(y, -0.2, +0.2)', 
            'output, reversed, hyst(y, +0.5, -0.5)'))
plt.title('Thresholding with hysteresis')
plt.show()

Answer 3

可以这样做：

def hyst(x, th_lo, th_hi, initial = False):
    hi = x >= th_hi
    lo_or_hi = (x <= th_lo) | 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(ind)
    return np.where(cnt, hi[ind[cnt-1]], initial)

解释一下：ind 是所有样本的索引，这些样本的信号要么低于下限，要么高于上限，因此“开关”的位置就很明确。通过 cumsum，你可以创建一个计数器，它指向最后一个明确的样本的索引。如果输入向量的开始部分在两个阈值之间，cnt 就会是 0，这时你需要用 where 函数将对应的输出设置为初始值。

致谢：这是我在一个老帖子中发现的技巧，原帖是在某个 Matlab 论坛上，我把它翻译成了 Numpy。这个代码有点难理解，而且需要分配各种中间数组。如果 Numpy 能提供一个专门的函数，类似于你简单的 for 循环，但用 C 实现以提高速度，那就更好了。

快速测试：

x = np.linspace(0,20, 1000)
y = np.sin(x)
h1 = hyst(y, -0.5, 0.5)
h2 = hyst(y, -0.5, 0.5, True)
plt.plot(x, y, x, -0.5 + h1, x, -0.5 + h2)
plt.legend(('input', 'output, start=0', 'output, start=1'))
plt.title('Thresholding with hysteresis')
plt.show()

结果：