一维数据中的步长检测答案

【问题标题】：Step detection in one-dimensional data一维数据中的步长检测
【发布时间】：2018-06-08 14:23:09
【问题描述】：

在 Python 中是否存在用于检测一维数据中的步骤的现有实现？

例如在此数据中检测到一个步骤的东西：

算法out there有不少描述，但我想知道Python中是否存在适合该工作的东西？

我不确定我是否应该/如何提供这些数据，但它是：

[ 594.          568.55555556  577.22222222  624.55555556  546.66666667
552.88888889  575.55555556  592.33333333  528.88888889  576.11111111
625.          574.22222222  556.33333333  567.66666667  576.66666667
591.66666667  566.33333333  567.33333333  547.44444444  631.11111111
555.66666667  548.66666667  579.44444444  546.88888889  597.55555556
519.88888889  582.33333333  618.88888889  574.55555556  547.44444444
593.11111111  565.66666667  544.66666667  562.66666667  554.11111111
543.88888889  602.33333333  609.77777778  550.55555556  561.88888889
719.33333333  784.44444444  711.22222222  843.66666667  691.33333333
690.11111111  684.33333333  749.11111111  759.11111111  653.33333333
817.11111111  705.22222222  689.44444444  712.33333333  659.
683.88888889  713.          740.44444444  692.22222222  677.33333333
681.44444444  640.          717.55555556  717.88888889  769.22222222
690.88888889  786.          774.66666667  799.44444444  743.44444444
789.88888889  673.66666667  685.66666667  709.88888889  645.55555556
846.11111111  792.77777778  702.22222222  749.44444444  678.55555556
707.55555556  665.77777778  643.55555556  671.44444444  795.66666667
627.22222222  684.55555556  708.44444444  829.66666667  719.        ]

【问题讨论】：

这是一个非常有趣的问题，但遗憾的是，图书馆推荐与 SO 无关。不过，请给我点赞。
@errantlinguist 我稍微改变了措辞，以免直接要求模块。希望它不会被删除！
这让我想起了最近发的一个类似问题：stackoverflow.com/questions/47290732/…
"相当多的算法描述"链接到期刊文章，而不是算法本身。你想要的python算法类型很可能在scipy中找到
既然你有。标记的信号处理scipy-signal 可能是您进行自己研究的好地方。您确实需要了解自己的目标。您提出的问题远非微不足道，而且正如指出的那样不太适合 SO。

标签： python signal-processing

【解决方案1】：

与一步卷积，看看峰值分辨率是否足够好

import numpy as np
from matplotlib import pyplot as plt


d = '''594.          568.55555556  577.22222222  624.55555556  546.66666667
552.88888889  575.55555556  592.33333333  528.88888889  576.11111111
625.          574.22222222  556.33333333  567.66666667  576.66666667
591.66666667  566.33333333  567.33333333  547.44444444  631.11111111
555.66666667  548.66666667  579.44444444  546.88888889  597.55555556
519.88888889  582.33333333  618.88888889  574.55555556  547.44444444
593.11111111  565.66666667  544.66666667  562.66666667  554.11111111
543.88888889  602.33333333  609.77777778  550.55555556  561.88888889
719.33333333  784.44444444  711.22222222  843.66666667  691.33333333
690.11111111  684.33333333  749.11111111  759.11111111  653.33333333
817.11111111  705.22222222  689.44444444  712.33333333  659.
683.88888889  713.          740.44444444  692.22222222  677.33333333
681.44444444  640.          717.55555556  717.88888889  769.22222222
690.88888889  786.          774.66666667  799.44444444  743.44444444
789.88888889  673.66666667  685.66666667  709.88888889  645.55555556
846.11111111  792.77777778  702.22222222  749.44444444  678.55555556
707.55555556  665.77777778  643.55555556  671.44444444  795.66666667
627.22222222  684.55555556  708.44444444  829.66666667  719.        '''

dary = np.array([*map(float, d.split())])

dary -= np.average(dary)

step = np.hstack((np.ones(len(dary)), -1*np.ones(len(dary))))

dary_step = np.convolve(dary, step, mode='valid')

# get the peak of the convolution, its index

step_indx = np.argmax(dary_step)  # yes, cleaner than np.where(dary_step == dary_step.max())[0][0]

# plots

plt.plot(dary)

plt.plot(dary_step/10)

plt.plot((step_indx, step_indx), (dary_step[step_indx]/10, 0), 'r')

【讨论】：

这太不可思议了。您是否知道任何可用的在线资源或涵盖此类方法的书籍？我指的不仅仅是卷积或信号处理方法，而是应用信号处理来检测异常值、异常和状态变化，就像你刚才描述的那样？
argmax 可能优于np.where(dary_step == dary_step.max())[0][0]
conolution 是 EE 的中流砥柱，我的背景 - 根据卷积步骤的相对极性，峰值也可能为负 - 所以一般来说需要更多的案例测试跨度>
您也可以使用-2 * np.cumsum(dary) 代替np.convolve(dary, step, mode='valid')。速度有点快。
@datapug：零均值转换是一种使用argmax() 很容易找到峰值的技巧（看看没有转换会发生什么）。基本上，带有阶梯核的卷积是累积和（在这种情况下是缩放的），因此零均值变换确保累积和以零结尾。请注意，这种步骤检测方法只有在您已经知道只有一个步骤时才能可靠地工作。另见例如the DSP guide.

【解决方案2】：

我自己对此很感兴趣。我不是专家，但正如 in this answer 建议的那样，如果您先对信号进行去噪，它可能有助于区分真正的步骤和噪声。 scikit-image 中的 total variation denoising 算法可以解决问题：

import numpy as np
from matplotlib import pyplot as plt
from skimage.restoration import denoise_tv_chambolle


data = '''594.          568.55555556  577.22222222  624.55555556  546.66666667
552.88888889  575.55555556  592.33333333  528.88888889  576.11111111
625.          574.22222222  556.33333333  567.66666667  576.66666667
591.66666667  566.33333333  567.33333333  547.44444444  631.11111111
555.66666667  548.66666667  579.44444444  546.88888889  597.55555556
519.88888889  582.33333333  618.88888889  574.55555556  547.44444444
593.11111111  565.66666667  544.66666667  562.66666667  554.11111111
543.88888889  602.33333333  609.77777778  550.55555556  561.88888889
719.33333333  784.44444444  711.22222222  843.66666667  691.33333333
690.11111111  684.33333333  749.11111111  759.11111111  653.33333333
817.11111111  705.22222222  689.44444444  712.33333333  659.
683.88888889  713.          740.44444444  692.22222222  677.33333333
681.44444444  640.          717.55555556  717.88888889  769.22222222
690.88888889  786.          774.66666667  799.44444444  743.44444444
789.88888889  673.66666667  685.66666667  709.88888889  645.55555556
846.11111111  792.77777778  702.22222222  749.44444444  678.55555556
707.55555556  665.77777778  643.55555556  671.44444444  795.66666667
627.22222222  684.55555556  708.44444444  829.66666667  719.        '''


x = np.array(data.split()).astype('float')
x_std = (x - x.mean()) / x.std()
x_denoise = denoise_tv_chambolle(x_std, weight=1)  # adjust the parameters
x_step = -2*np.cumsum(x_denoise)
step_indicator = x_step == x_step.max()

n = x.shape[0]
plt.subplot(211)
plt.plot(range(n), x_std, label='standardized')
plt.plot(range(n), x_denoise, label='denoised (TV)')
plt.legend()
plt.grid()
plt.subplot(212)
plt.step(range(n), step_indicator)
plt.show()

输出：

注意：

对这些方法的一个警告是，您必须确保信号中只有一个步骤（并且必须有一个）。

【讨论】：