查找曲线上斜率变化的位置

Question

我有创建如下曲线的时间和电压数据点。

time 数据是

array([  0.10810811,   0.75675676,   1.62162162,   2.59459459,
         3.56756757,   4.21621622,   4.97297297,   4.97297297,
         4.97297297,   4.97297297,   4.97297297,   4.97297297,
         4.97297297,   4.97297297,   5.08108108,   5.18918919,
         5.2972973 ,   5.51351351,   5.72972973,   5.94594595,
         6.27027027,   6.59459459,   7.13513514,   7.67567568,
         8.32432432,   9.18918919,  10.05405405,  10.91891892,
        11.78378378,  12.64864865,  13.51351351,  14.37837838,
        15.35135135,  16.32432432,  17.08108108,  18.16216216,
        19.02702703,  20.        ,  20.        ,  20.        ,
        20.        ,  20.        ,  20.        ,  20.        ,
        20.10810811,  20.21621622,  20.43243243,  20.64864865,
        20.97297297,  21.40540541,  22.05405405,  22.91891892,
        23.78378378,  24.86486486,  25.83783784,  26.7027027 ,
        27.56756757,  28.54054054,  29.51351351,  30.48648649,
        31.56756757,  32.64864865,  33.62162162,  34.59459459,
        35.67567568,  36.64864865,  37.62162162,  38.59459459,
        39.67567568,  40.75675676,  41.83783784,  42.81081081,
        43.89189189,  44.97297297,  46.05405405,  47.02702703,
        48.10810811,  49.18918919,  50.27027027,  51.35135135,
        52.43243243,  53.51351351,  54.48648649,  55.56756757,
        56.75675676,  57.72972973,  58.81081081,  59.89189189])

volts 数据是

array([ 4.11041056,  4.11041056,  4.11041056,  4.11041056,  4.11041056,
        4.11041056,  4.11041056,  4.10454545,  4.09794721,  4.09208211,
        4.08621701,  4.07961877,  4.07228739,  4.06568915,  4.05909091,
        4.05175953,  4.04516129,  4.03782991,  4.03123167,  4.02463343,
        4.01803519,  4.01217009,  4.00557185,  3.99970674,  3.99384164,
        3.98797654,  3.98284457,  3.97771261,  3.97331378,  3.96891496,
        3.96451613,  3.96085044,  3.95645161,  3.95205279,  3.9483871 ,
        3.94398827,  3.94032258,  3.93665689,  3.94325513,  3.94985337,
        3.95645161,  3.96378299,  3.97038123,  3.97624633,  3.98284457,
        3.98944282,  3.99604106,  4.0026393 ,  4.00923754,  4.01510264,
        4.02096774,  4.02609971,  4.02903226,  4.03196481,  4.03416422,
        4.0356305 ,  4.03709677,  4.03856305,  4.03929619,  4.04002933,
        4.04076246,  4.04222874,  4.04296188,  4.04296188,  4.04369501,
        4.04442815,  4.04516129,  4.04516129,  4.04589443,  4.04589443,
        4.04662757,  4.04662757,  4.0473607 ,  4.0473607 ,  4.04809384,
        4.04809384,  4.04809384,  4.04882698,  4.04882698,  4.04882698,
        4.04956012,  4.04956012,  4.04956012,  4.04956012,  4.05029326,
        4.05029326,  4.05029326,  4.05029326])

我想确定标记为 A、B、C、D 和 E 的点的位置。点 A 是斜率从零变为未定义的第一个位置。 B点是线不再垂直的位置。 C点是曲线的最小值。点 D 是曲线不再垂直的位置。 E点是斜率再次接近零的地方。下面的 Python 代码确定了点 A 和 C 的位置。

tdiff = np.diff(time)
vdiff = np.diff(volts)

# point A
idxA = np.where(vdiff < 0)[0][0]
timeA = time[idxA]
voltA = volts[idxA]

# point C
idxC = volts.idxmin()
timeC = time[idxC]
voltC = volts[idxC]

如何确定曲线上由点 B、D 和 E 表示的其他位置？

Answer 1

您正在寻找标记斜率变化为零或无穷大的任何位置的点。我们实际上不需要在任何地方计算斜率：y<sub>n</sub> - y<sub>n-1</sub> == 0 和 y<sub>n+1</sub> - y<sub>n</sub> != 0，反之亦然，或者 x 相同。

我们可以获取x 的差异。如果两个连续元素之一为零，则差异的差异将是该点的差异或负差异。因此，我们只想找到并标记diff(x) == diff(diff(x)) 和diff(x) != 0 的所有点，并根据数组之间的大小差异进行适当调整。我们还想要y 的所有点。

在 numpy 术语中，这可以写成如下

def masks(vec):
    d = np.diff(vec)
    dd = np.diff(d)

    # Mask of locations where graph goes to vertical or horizontal, depending on vec
    to_mask = ((d[:-1] != 0) & (d[:-1] == -dd))
    # Mask of locations where graph comes from vertical or horizontal, depending on vec
    from_mask = ((d[1:] != 0) & (d[1:] == dd))
    return to_mask, from_mask

to_vert_mask, from_vert_mask = masks(time)
to_horiz_mask, from_horiz_mask = masks(volts)

请记住，掩码是根据二阶差分计算的，因此它们比输入短两个元素。掩码中的元素对应于输入数组中的元素，在前缘和后缘有一个元素的边框（因此下面的索引[1:-1]）。您可以使用np.nonzero 将掩码转换为索引，也可以直接使用掩码作为索引获取 x 和 y 值：

def apply_mask(mask, x, y):
    return x[1:-1][mask], y[1:-1][mask]

to_vert_t, to_vert_v = apply_mask(to_vert_mask, time, volts)
from_vert_t, from_vert_v = apply_mask(from_vert_mask, time, volts)
to_horiz_t, to_horiz_v = apply_mask(to_horiz_mask, time, volts)
from_horiz_t, from_horiz_v = apply_mask(from_horiz_mask, time, volts)

plt.plot(time, volts, 'b-')
plt.plot(to_vert_t, to_vert_v, 'r^', label='Plot goes vertical')
plt.plot(from_vert_t, from_vert_v, 'kv', label='Plot stops being vertical')
plt.plot(to_horiz_t, to_horiz_v, 'r>', label='Plot goes horizontal')
plt.plot(from_horiz_t, from_horiz_v, 'k<', label='Plot stops being horizontal')
plt.legend()
plt.show()

这是结果图：

请注意，由于分类是单独进行的，因此“A 点”被正确识别为垂直开始和水平结束的点。问题是根据这些标准，“E 点”似乎无法解决。放大显示所有增殖点都正确识别水平线段：

您可以通过完全丢弃from_horiz 并仅丢弃to_horiz 中的最后一个值来选择“正确”版本的“E 点”：

to_horiz_t, to_horiz_v = apply_mask(to_horiz_mask, time, volts)
to_horiz_t, to_horiz_v = to_horiz_t[-1], to_horiz_v[-1]

plt.plot(time, volts, 'b-')
plt.plot(*apply_mask(to_vert_mask, time, volts), 'r^', label='Plot goes vertical')
plt.plot(*apply_mask(from_vert_mask, time, volts), 'kv', label='Plot stops being vertical')
plt.plot(to_horiz_t, to_horiz_v, 'r>', label='Plot goes horizontal')
plt.legend()
plt.show()

我将其用作apply_mask 结果的明星扩展展示。结果图是：

这几乎正是您要寻找的情节。丢弃from_horiz 也会使“A 点”仅被识别为垂直下降，这很好。

正如to_horiz 中的多个值所示，此方法对数据中的噪声非常敏感。您的数据非常平滑，但这种方法不太可能适用于未经过滤的原始测量。