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

Posted 2023-02-22

【中文标题】如何找到带有滞后的零交叉？

【英文标题】：How to find zero crossings with hysteresis?

【发布时间】：2014-06-10 23:34:45

【问题描述】：

在 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"

我想要一种无需显式循环遍历数组的方法：但我想不出任何方法，因为每个状态值都取决于前一个状态。没有循环可以吗？

【参考方案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(x)
    return np.where(cnt, hi[ind[cnt-1]], initial)

解释：ind 是信号低于下限或高于上限阈值的所有样本的索引，因此“开关”的位置已明确定义。使用cumsum，您可以创建某种计数器，该计数器指向最后一个明确定义的样本的索引。如果输入向量的起点在两个阈值之间，则cnt将为0，因此需要使用where函数将对应的输出设置为初始值。

Credit：这是我在某个 Matlab 论坛上的 old post 中发现的一个技巧，我将其翻译为 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()

结果：

【讨论】：

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

我问这个是因为我想了解你的函数是做什么的，所以我可以用 C 编写代码......?

【参考方案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()