在找到峰值时查找时间序列数据中达到特定值的时间
Finding the time in which a specific value is reached in time-series data when peaks are found
我想在有噪声的时间序列数据中找到达到某个值的时刻。如果数据中没有峰值，我可以在 MATLAB 中执行以下操作。
代码来自 here
% create example data
d=1:100;
t=d/100;
ts = timeseries(d,t);
% define threshold
thr = 55;
data = ts.data(:);
time = ts.time(:);
ind = find(data>thr,1,'first');
time(ind) %time where data>threshold
但是当有噪音时，我不知道该怎么办。
在上图中绘制的时间序列数据中，我想找到达到 y 轴值 5 的时刻。数据实际上在 t>=100 s 时稳定为 5。但由于数据中存在噪声，我们看到一个峰值在 20 秒左右达到 5。我想知道如何检测例如 100 秒作为正确的时间而不是 20 秒。上面的代码 posted 只会给出 20 s 作为答案。我
看到 post here 解释了使用滑动 window 来查找数据何时平衡。但是，我不确定如何实施。建议将非常有帮助。
可以找到上图中绘制的示例数据here
关于如何在 Python 或 MATLAB 代码中实现的建议将非常有帮助。
编辑：
我不想在峰值 (/noise/overshoot) 出现时捕获。我想找到达到平衡的时间。例如，在 20 秒左右，曲线上升并下降到 5 以下。大约 100 秒后，曲线平衡到稳态值 5，并且从不下降或达到峰值。
假设您的数据在 data 变量中，时间索引在 time 中。那么
import numpy as np
threshold = 0.025
stable_index = np.where(np.abs(data[-1] - data) > threshold)[0][-1] + 1
print('Stabilizes after', time[stable_index], 'sec')
Stabilizes after 96.6 sec
这里 data[-1] - data 是 data 的最后一个值与所有数据值之间的差异。这里的假设是data的最后一个值代表平衡点。
np.where( * > threshold )[0] 是 data 的值的所有指数，这些指数 大于 阈值，但仍未稳定。我们只取最后一个索引。下一个是时间序列被认为是稳定的，因此 + 1.
精确的数据分析是一项严肃的工作（也是我的热情），涉及对您正在研究的系统的大量理解。以下是评论，不幸的是，我怀疑是否有一个简单的好答案可以解决您的问题——您将不得不考虑一下。数据分析基本上总是需要“讨论”。
首先针对你的数据和问题的总体情况：
当您谈论噪声时，在数据分析中这意味着统计 运行dom 波动。最常见的是高斯分布（有时还有其他分布，例如泊松分布）。高斯噪声是 a) 每个 bin 中的 运行dom 和 b) 在负方向和正方向上对称。因此，您在 20 秒左右的峰值处观察到的不是噪声。与运行dom noise相比，它具有非常不同，非常系统和扩展的特性。这是一个必须有来源的“神器”，但我们只能在这里推测。在 real-world 应用程序中，研究和删除此类工件是最昂贵的 time-consuming 任务。
查看您的数据，运行dom 噪声可以忽略不计。这是非常精确的数据。例如，在 ~150 秒及之后，没有可见的 运行dom 波动高达第四位小数。
在得出结论这不是常识上的噪音之后，它可能至少有两件事：a) 你正在研究的系统的一个特征，因此，你可以开发一个 model/formula 你可以“适合”数据。 b) 测量链中某处带宽有限的特性，因此，这里有一个 high-frequency 截止点。参见例如https://en.wikipedia.org/wiki/Ringing_artifacts。不幸的是，对于 a 和 b，都没有 catch-all 通用解决方案。而且你的问题描述（即使有代码和数据）也不足以提出理想的方法。
现在花了 ~1 小时处理数据并绘制了一些图。我相信（推测）~10s 处的极其尖锐的特征不可能是数据的“物理”属性。简直太extreme/steep了。这里发生了一些根本性的事情。我的猜测可能是某些设备刚刚打开（之前关闭）。这样，之前的数据是没有意义的，之后还有很短的时间来稳定系统。在这种情况下，除了在系统稳定在 40 秒左右之前完全丢弃数据外，别无选择。这也使您的问题变得微不足道。只需删除前 40 秒，然后最大值就很明显了。
那么您可以使用哪些技术解决方案，请不要太沮丧，因为您必须自己考虑这个问题 assemble 最适合您的情况的解决方案。我将您的数据复制到两个 numpy 数组 x 和 y 以及 运行 中 python 中的以下测试：
去掉不稳定的时间
这是简单的解决方案 -- 我更喜欢它。
plt.figure()
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(x, y, label="original")
y_cut = y
y_cut[:40] = 0
plt.plot(x, y_cut, label="cut 40s")
plt.legend()
plt.grid()
plt.show()
注意如果你有点疯狂（关于数据），请继续阅读下面的内容。
滑动window
您提到了“滑动 window”，它最适合 运行dom 噪声（您没有）或周期性波动（您也没有）。滑动 window 只是对连续的 bin 进行平均，平均掉 运行dom 波动。从数学上讲，这是一个卷积。
从技术上讲，您实际上可以这样解决您的问题（您自己尝试更大的 Nwindow 值）：
Nwindow=10
y_slide_10 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
Nwindow=20
y_slide_20 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
Nwindow=30
y_slide_30 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(x,y, label="original")
plt.plot(x,y_slide_10, label="window=10")
plt.plot(x,y_slide_20, label='window=20')
plt.plot(x,y_slide_30, label='window=30')
plt.legend()
#plt.xscale('log') # useful
plt.grid()
plt.show()
因此，从技术上讲，您可以成功抑制最初的“驼峰”。但不要忘记这是一个 hand-tuned 而不是通用的解决方案...
任何滑动 window 解决方案的另一个警告：这总是会扭曲您的时间安排。由于您根据上升或下降信号对时间间隔进行平均，因此您的复杂轨迹会及时移动 back/forth（轻微但显着）。在您的特定情况下，这不是问题，因为主要信号区域基本上没有 time-dependence（非常平坦）。
频域
这应该是灵丹妙药，但它也不适用于您的示例 well/easily。事实上这并没有更好地工作是对我的主要暗示，即前 40 秒的数据最好被丢弃....（即在科学工作中）
您可以使用快速傅立叶 t运行sform 检查 frequency-domain 中的数据。
import scipy.fft
y_fft = scipy.fft.rfft(y)
# original frequency domain plot
plt.plot(y_fft, label="original")
plt.xlabel('frequency')
plt.ylabel('signal')
plt.yscale('log')
plt.show()
频率结构代表数据的特征。峰值零是 ~100 秒后的稳定区域，驼峰与时间的（快速）变化相关。您现在可以尝试并更改频谱（--> 滤波器），但我认为频谱是如此人为，因此在这里不会产生很好的结果。与其他数据一起尝试，您可能会印象深刻！我尝试了两件事，首先将 high-frequency 区域切出（设置为零），其次，在频域中应用 sliding-window 滤波器（保留 0 处的峰值，因为无法触及。试试看你知道为什么.
# cut high-frequency by setting to zero
y_fft_2 = np.array(y_fft)
y_fft_2[50:70] = 0
# sliding window in frequency
Nwindow = 15
Start = 10
y_fft_slide = np.array(y_fft)
y_fft_slide[Start:] = np.convolve(y_fft[Start:], np.ones((Nwindow,))/Nwindow, mode='same')
# frequency-domain plot
plt.plot(y_fft, label="original")
plt.plot(y_fft_2, label="high-frequency, filter")
plt.plot(y_fft_slide, label="frequency sliding window")
plt.xlabel('frequency')
plt.ylabel('signal')
plt.yscale('log')
plt.legend()
plt.show()
将其转换回 time-domain:
# reverse FFT into time-domain for plotting
y_filtered = scipy.fft.irfft(y_fft_2)
y_filtered_slide = scipy.fft.irfft(y_fft_slide)
# time-domain plot
plt.plot(x[:500], y[:500], label="original")
plt.plot(x[:500], y_filtered[:500], label="high-f filtered")
plt.plot(x[:500], y_filtered_slide[:500], label="frequency sliding window")
# plt.xscale('log') # useful
plt.grid()
plt.legend()
plt.show()
产量
这些解决方案中存在明显的波动，这使得它们对您的目的基本上毫无用处。这使我进入最后一个练习，再次在“频率滑动 window”time-domain
上应用 sliding-window 过滤器
# extra time-domain sliding window
Nwindow=90
y_fft_90 = np.convolve(y_filtered_slide, np.ones((Nwindow,))/Nwindow, mode='same')
# final time-domain plot
plt.plot(x[:500], y[:500], label="original")
plt.plot(x[:500], y_fft_90[:500], label="frequency-sliding window, slide")
# plt.xscale('log') # useful
plt.legend()
plt.show()
我对这个结果很满意，但它仍然有很小的振荡，因此没有解决你原来的问题。
结论
多么有趣。浪费了一个小时。也许它对某人有用。甚至对你来说，娜塔莎。请不要生我的气...
如果您处理的是最终单调收敛到某个固定值的确定性数据，那么问题就非常简单了。您的最后一次观察应该最接近极限，因此您可以相对于最后一个数据点定义一个可接受的容差阈值，并从后向前扫描您的数据以找出超出阈值的位置。
一旦将随机噪声添加到图片中，事情就会变得更加糟糕，尤其是在存在序列相关的情况下。这个问题在仿真建模中很常见（见下面的（*）），被称为initial bias. It was first identified by Conway in 1963的问题，从那时起一直是一个活跃的研究领域，没有关于如何处理它的普遍接受的明确答案。与确定性情况一样，最广泛接受的答案从数据集的右侧开始解决问题，因为这是数据最有可能处于稳定状态的地方。基于这种方法的技术使用数据集的末尾来建立某种统计标准或基线，以衡量随着向数据集前端移动添加观察结果，数据开始看起来明显不同的位置。由于序列相关的存在，这变得非常复杂。
如果时间序列处于稳定状态，从 covariance stationary 的意义上讲，则数据的简单平均值是对其期望值的无偏估计，但估计平均值的标准误差在很大程度上取决于关于序列相关。正确的标准误差平方不再是 s2/n，而是 (s2/n)*W，其中 W 是正确的自相关值的加权和。 1990 年代开发了一种称为 MSER 的方法，它通过尝试确定标准误差最小化的位置来避免尝试正确估计 W 的问题。给定足够大的样本量，它将 W 视为事实上的常数，因此如果您考虑两个标准误差估计值的比率，W 的抵消和最小值出现在 s2/n 的位置被最小化。 MSER 进行如下：
- 从末尾开始，计算一半数据集的s2，建立基线。
- 现在更新 s2 的估计值，使用 Welford's online algorithm 等有效技术一次一个观测值，计算 s2/n 其中 n 是到目前为止统计的观察次数。追踪哪个 n 值产生最小的 s2/n。起泡、冲洗、重复。
- 一旦你从后向前遍历了整个数据集，产生最小s的n2/n就是数据末尾的观测数设置无法检测到的起始条件偏差。
理由 - 如果基线足够大（一半的数据），只要时间序列保持稳定状态，s2/n 应该相对稳定。由于 n 是单调递增的，因此 s2/n 应该继续减少，但受其作为估计值的可变性的限制。但是，一旦您开始获取不稳定状态的观测值，均值和方差的漂移就会使 s2/n 的分子膨胀。因此，最小值对应于没有非平稳性迹象的最后一次观察。可以在 proceedings paper. A Ruby implementation is available on BitBucket.
中找到更多详细信息
您的数据变化如此之小，以至于 MSER 断定它仍在收敛到稳定状态。因此，我建议采用第一段中概述的确定性方法。如果您以后有嘈杂的数据，我绝对建议您试一试 MSER。
(*) - 简而言之，仿真模型是一个计算机程序，因此必须将其状态设置为一组初始值。我们通常不知道系统状态在 long 运行 中会是什么样子，所以我们将其初始化为一组任意但方便的值，然后让系统“预热”。问题是模拟的初始结果不是稳态行为的典型结果，因此在您的分析中包含该数据会使它们产生偏差。解决方案是去除数据中有偏差的部分，但应该去除多少？
我想在有噪声的时间序列数据中找到达到某个值的时刻。如果数据中没有峰值，我可以在 MATLAB 中执行以下操作。
代码来自 here
% create example data
d=1:100;
t=d/100;
ts = timeseries(d,t);
% define threshold
thr = 55;
data = ts.data(:);
time = ts.time(:);
ind = find(data>thr,1,'first');
time(ind) %time where data>threshold
但是当有噪音时，我不知道该怎么办。
在上图中绘制的时间序列数据中，我想找到达到 y 轴值 5 的时刻。数据实际上在 t>=100 s 时稳定为 5。但由于数据中存在噪声，我们看到一个峰值在 20 秒左右达到 5。我想知道如何检测例如 100 秒作为正确的时间而不是 20 秒。上面的代码 posted 只会给出 20 s 作为答案。我 看到 post here 解释了使用滑动 window 来查找数据何时平衡。但是，我不确定如何实施。建议将非常有帮助。
可以找到上图中绘制的示例数据here
关于如何在 Python 或 MATLAB 代码中实现的建议将非常有帮助。
编辑： 我不想在峰值 (/noise/overshoot) 出现时捕获。我想找到达到平衡的时间。例如，在 20 秒左右，曲线上升并下降到 5 以下。大约 100 秒后，曲线平衡到稳态值 5，并且从不下降或达到峰值。
假设您的数据在 data 变量中，时间索引在 time 中。那么
import numpy as np
threshold = 0.025
stable_index = np.where(np.abs(data[-1] - data) > threshold)[0][-1] + 1
print('Stabilizes after', time[stable_index], 'sec')
Stabilizes after 96.6 sec
这里 data[-1] - data 是 data 的最后一个值与所有数据值之间的差异。这里的假设是data的最后一个值代表平衡点。
np.where( * > threshold )[0] 是 data 的值的所有指数，这些指数 大于 阈值，但仍未稳定。我们只取最后一个索引。下一个是时间序列被认为是稳定的，因此 + 1.
精确的数据分析是一项严肃的工作（也是我的热情），涉及对您正在研究的系统的大量理解。以下是评论，不幸的是，我怀疑是否有一个简单的好答案可以解决您的问题——您将不得不考虑一下。数据分析基本上总是需要“讨论”。
首先针对你的数据和问题的总体情况：
当您谈论噪声时，在数据分析中这意味着统计 运行dom 波动。最常见的是高斯分布（有时还有其他分布，例如泊松分布）。高斯噪声是 a) 每个 bin 中的 运行dom 和 b) 在负方向和正方向上对称。因此，您在 20 秒左右的峰值处观察到的不是噪声。与运行dom noise相比，它具有非常不同，非常系统和扩展的特性。这是一个必须有来源的“神器”，但我们只能在这里推测。在 real-world 应用程序中，研究和删除此类工件是最昂贵的 time-consuming 任务。
查看您的数据，运行dom 噪声可以忽略不计。这是非常精确的数据。例如，在 ~150 秒及之后，没有可见的 运行dom 波动高达第四位小数。
在得出结论这不是常识上的噪音之后，它可能至少有两件事：a) 你正在研究的系统的一个特征，因此，你可以开发一个 model/formula 你可以“适合”数据。 b) 测量链中某处带宽有限的特性，因此，这里有一个 high-frequency 截止点。参见例如https://en.wikipedia.org/wiki/Ringing_artifacts。不幸的是，对于 a 和 b，都没有 catch-all 通用解决方案。而且你的问题描述（即使有代码和数据）也不足以提出理想的方法。
现在花了 ~1 小时处理数据并绘制了一些图。我相信（推测）~10s 处的极其尖锐的特征不可能是数据的“物理”属性。简直太extreme/steep了。这里发生了一些根本性的事情。我的猜测可能是某些设备刚刚打开（之前关闭）。这样，之前的数据是没有意义的，之后还有很短的时间来稳定系统。在这种情况下，除了在系统稳定在 40 秒左右之前完全丢弃数据外，别无选择。这也使您的问题变得微不足道。只需删除前 40 秒，然后最大值就很明显了。
那么您可以使用哪些技术解决方案，请不要太沮丧，因为您必须自己考虑这个问题 assemble 最适合您的情况的解决方案。我将您的数据复制到两个 numpy 数组 x 和 y 以及 运行 中 python 中的以下测试：
去掉不稳定的时间
这是简单的解决方案 -- 我更喜欢它。
plt.figure()
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(x, y, label="original")
y_cut = y
y_cut[:40] = 0
plt.plot(x, y_cut, label="cut 40s")
plt.legend()
plt.grid()
plt.show()
注意如果你有点疯狂（关于数据），请继续阅读下面的内容。
滑动window
您提到了“滑动 window”，它最适合 运行dom 噪声（您没有）或周期性波动（您也没有）。滑动 window 只是对连续的 bin 进行平均，平均掉 运行dom 波动。从数学上讲，这是一个卷积。
从技术上讲，您实际上可以这样解决您的问题（您自己尝试更大的 Nwindow 值）：
Nwindow=10
y_slide_10 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
Nwindow=20
y_slide_20 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
Nwindow=30
y_slide_30 = np.convolve(y, np.ones((Nwindow,))/Nwindow, mode='same')
plt.xlabel('time')
plt.ylabel('signal')
plt.plot(x,y, label="original")
plt.plot(x,y_slide_10, label="window=10")
plt.plot(x,y_slide_20, label='window=20')
plt.plot(x,y_slide_30, label='window=30')
plt.legend()
#plt.xscale('log') # useful
plt.grid()
plt.show()
因此，从技术上讲，您可以成功抑制最初的“驼峰”。但不要忘记这是一个 hand-tuned 而不是通用的解决方案...
任何滑动 window 解决方案的另一个警告：这总是会扭曲您的时间安排。由于您根据上升或下降信号对时间间隔进行平均，因此您的复杂轨迹会及时移动 back/forth（轻微但显着）。在您的特定情况下，这不是问题，因为主要信号区域基本上没有 time-dependence（非常平坦）。
频域
这应该是灵丹妙药，但它也不适用于您的示例 well/easily。事实上这并没有更好地工作是对我的主要暗示，即前 40 秒的数据最好被丢弃....（即在科学工作中）
您可以使用快速傅立叶 t运行sform 检查 frequency-domain 中的数据。
import scipy.fft
y_fft = scipy.fft.rfft(y)
# original frequency domain plot
plt.plot(y_fft, label="original")
plt.xlabel('frequency')
plt.ylabel('signal')
plt.yscale('log')
plt.show()
频率结构代表数据的特征。峰值零是 ~100 秒后的稳定区域，驼峰与时间的（快速）变化相关。您现在可以尝试并更改频谱（--> 滤波器），但我认为频谱是如此人为，因此在这里不会产生很好的结果。与其他数据一起尝试，您可能会印象深刻！我尝试了两件事，首先将 high-frequency 区域切出（设置为零），其次，在频域中应用 sliding-window 滤波器（保留 0 处的峰值，因为无法触及。试试看你知道为什么.
# cut high-frequency by setting to zero
y_fft_2 = np.array(y_fft)
y_fft_2[50:70] = 0
# sliding window in frequency
Nwindow = 15
Start = 10
y_fft_slide = np.array(y_fft)
y_fft_slide[Start:] = np.convolve(y_fft[Start:], np.ones((Nwindow,))/Nwindow, mode='same')
# frequency-domain plot
plt.plot(y_fft, label="original")
plt.plot(y_fft_2, label="high-frequency, filter")
plt.plot(y_fft_slide, label="frequency sliding window")
plt.xlabel('frequency')
plt.ylabel('signal')
plt.yscale('log')
plt.legend()
plt.show()
将其转换回 time-domain:
# reverse FFT into time-domain for plotting
y_filtered = scipy.fft.irfft(y_fft_2)
y_filtered_slide = scipy.fft.irfft(y_fft_slide)
# time-domain plot
plt.plot(x[:500], y[:500], label="original")
plt.plot(x[:500], y_filtered[:500], label="high-f filtered")
plt.plot(x[:500], y_filtered_slide[:500], label="frequency sliding window")
# plt.xscale('log') # useful
plt.grid()
plt.legend()
plt.show()
产量
这些解决方案中存在明显的波动，这使得它们对您的目的基本上毫无用处。这使我进入最后一个练习，再次在“频率滑动 window”time-domain
上应用 sliding-window 过滤器# extra time-domain sliding window
Nwindow=90
y_fft_90 = np.convolve(y_filtered_slide, np.ones((Nwindow,))/Nwindow, mode='same')
# final time-domain plot
plt.plot(x[:500], y[:500], label="original")
plt.plot(x[:500], y_fft_90[:500], label="frequency-sliding window, slide")
# plt.xscale('log') # useful
plt.legend()
plt.show()
我对这个结果很满意，但它仍然有很小的振荡，因此没有解决你原来的问题。
结论
多么有趣。浪费了一个小时。也许它对某人有用。甚至对你来说，娜塔莎。请不要生我的气...
如果您处理的是最终单调收敛到某个固定值的确定性数据，那么问题就非常简单了。您的最后一次观察应该最接近极限，因此您可以相对于最后一个数据点定义一个可接受的容差阈值，并从后向前扫描您的数据以找出超出阈值的位置。
一旦将随机噪声添加到图片中，事情就会变得更加糟糕，尤其是在存在序列相关的情况下。这个问题在仿真建模中很常见（见下面的（*）），被称为initial bias. It was first identified by Conway in 1963的问题，从那时起一直是一个活跃的研究领域，没有关于如何处理它的普遍接受的明确答案。与确定性情况一样，最广泛接受的答案从数据集的右侧开始解决问题，因为这是数据最有可能处于稳定状态的地方。基于这种方法的技术使用数据集的末尾来建立某种统计标准或基线，以衡量随着向数据集前端移动添加观察结果，数据开始看起来明显不同的位置。由于序列相关的存在，这变得非常复杂。
如果时间序列处于稳定状态，从 covariance stationary 的意义上讲，则数据的简单平均值是对其期望值的无偏估计，但估计平均值的标准误差在很大程度上取决于关于序列相关。正确的标准误差平方不再是 s2/n，而是 (s2/n)*W，其中 W 是正确的自相关值的加权和。 1990 年代开发了一种称为 MSER 的方法，它通过尝试确定标准误差最小化的位置来避免尝试正确估计 W 的问题。给定足够大的样本量，它将 W 视为事实上的常数，因此如果您考虑两个标准误差估计值的比率，W 的抵消和最小值出现在 s2/n 的位置被最小化。 MSER 进行如下：
- 从末尾开始，计算一半数据集的s2，建立基线。
- 现在更新 s2 的估计值，使用 Welford's online algorithm 等有效技术一次一个观测值，计算 s2/n 其中 n 是到目前为止统计的观察次数。追踪哪个 n 值产生最小的 s2/n。起泡、冲洗、重复。
- 一旦你从后向前遍历了整个数据集，产生最小s的n2/n就是数据末尾的观测数设置无法检测到的起始条件偏差。
理由 - 如果基线足够大（一半的数据），只要时间序列保持稳定状态，s2/n 应该相对稳定。由于 n 是单调递增的，因此 s2/n 应该继续减少，但受其作为估计值的可变性的限制。但是，一旦您开始获取不稳定状态的观测值，均值和方差的漂移就会使 s2/n 的分子膨胀。因此，最小值对应于没有非平稳性迹象的最后一次观察。可以在 proceedings paper. A Ruby implementation is available on BitBucket.
中找到更多详细信息您的数据变化如此之小，以至于 MSER 断定它仍在收敛到稳定状态。因此，我建议采用第一段中概述的确定性方法。如果您以后有嘈杂的数据，我绝对建议您试一试 MSER。
(*) - 简而言之，仿真模型是一个计算机程序，因此必须将其状态设置为一组初始值。我们通常不知道系统状态在 long 运行 中会是什么样子，所以我们将其初始化为一组任意但方便的值，然后让系统“预热”。问题是模拟的初始结果不是稳态行为的典型结果，因此在您的分析中包含该数据会使它们产生偏差。解决方案是去除数据中有偏差的部分，但应该去除多少？