如何使用滞后效应找到零交叉点？

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"

我希望找到一种不需要显式循环数组的方法来实现这个，但是我想不出任何方法，因为每个状态值都依赖于前一个状态。有没有可能在不使用循环的情况下完成？

- Alex I

Answer 1

这可以这样做：

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()

结果：

- Bas Swinckels

---

@SpaceDog 你可以这样做，但如果你正在使用C语言，最好编写一个类似于原始问题中的简单循环。我的答案中的技巧在Python中更快，因为它使用矢量化的numpy代码，而不是缓慢的Python循环。矢量化代码必须多次通过数据，而C中的简单循环可以一次完成所有操作。 - Bas Swinckels

我问这个是因为我想了解你的函数是做什么的，这样我才能用C语言编写代码... - Duck

谢谢！这个解决方案以一种创造性的方式使用了索引，让我想起了 APL。 - Tobia

这个解决方案非常巧妙，我花了一些时间才理解它的工作原理。关键是要明白有两层间接引用在进行 - cnt-1 是一个索引数组，指向 ind，而 ind 本身是一个索引数组，指向 lo_or_hi。cnt 可以被看作是“迄今为止我们已经遇到的外部值的数量”，由于 ind 是外部值的索引列表，所以 cnt-1 实际上给出了最后一个外部值的索引，以便“锁定”。非常令人费解。 - Etienne Dechamps

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()

- AMTQ

Answer 3

这里有一个解决方案，它实现了与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，完成！

- Etienne Dechamps