在找到峰值时查找时间序列数据中达到特定值的时间

如何解决在找到峰值时查找时间序列数据中达到特定值的时间

我想在有噪声的时间序列数据中找到达到某个值的时刻。如果数据中没有峰值，我可以在 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 s 左右达到 5。我想知道如何检测例如 100 秒作为正确的时间而不是 20 s 。上面发布的代码只会给出 20 秒作为答案。一世 看到一篇文章 here 解释了使用滑动窗口来查找数据何时达到平衡。但是，我不确定如何实现相同的功能。建议会很有帮助。

可以找到上图中绘制的示例数据here

有关如何在 Python 或 MATLAB 代码中实现的建议将非常有帮助。

编辑： 我不想在峰值（/噪声/过冲）发生时捕获。我想找到达到平衡的时间。例如，大约 20 秒后，曲线上升和下降到 5 以下。大约 100 秒后，曲线平衡到稳态值 5，不再下降或达到峰值。

解决方法

精确的数据分析是一项严肃的工作（也是我的热情所在），它涉及对您正在研究的系统的大量理解。以下是评论，不幸的是，我怀疑您的问题是否有一个简单的好答案——您将不得不考虑它。数据分析基本上总是需要“讨论”。

首先是您的数据和一般问题：

1. 当您谈到噪声时，在数据分析中这意味着统计随机波动。最常见的是高斯分布（有时还有其他分布，例如 Poission）。高斯噪声 a) 在每个 bin 中是随机的，并且 b) 在正负方向上对称。因此，您在大约 20 秒的峰值中观察到的不是噪声。与随机噪声相比，它具有非常不同、非常系统和扩展的特性。这是一件必有渊源的“神器”，但我们只能在这里推测。在实际应用中，研究和移除此类工件是最昂贵和耗时的任务。

2. 查看您的数据，随机噪声可以忽略不计。这是非常精确的数据。例如，在大约 150 秒之后，直到第四个十进制数都没有可见的随机波动。

3. 在得出这不是常识中的噪音之后，它可能至少是两件事：a) 您正在研究的系统的一个特征，因此，您可以为它开发模型/公式的东西你可以“适合”数据。 b) 测量链中某处有限带宽的特性，因此，这里是高频截止。见例如https://en.wikipedia.org/wiki/Ringing_artifacts。不幸的是，对于 a 和 b，都没有包罗万象的通用解决方案。而且您的问题描述（即使有代码和数据）也不足以提出理想的方法。

4. 现在花了大约一个小时处理您的数据并绘制了一些图。我相信（推测）大约 10 秒的极其锐利的特征不能是数据的“物理”属性。它只是过于极端/陡峭。这里发生了一些根本性的事情。我的猜测可能是某些设备刚刚打开（之前关闭）。因此，之前的数据毫无意义，之后有很短的时间来稳定系统。在这种情况下没有真正的替代方案，只能完全丢弃数据，直到系统稳定在 40 秒左右。这也使您的问题变得微不足道。只需删除前40s，然后最大值就很明显了。

那么您可以使用哪些技术解决方案，请不要太沮丧，因为您必须自己考虑并为您的案例组装最佳解决方案。我将您的数据复制到两个 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()

注意如果你有点疯狂（关于数据），请继续阅读下面的内容。

滑动窗口

您提到了“滑动窗口”，它最适合随机噪声（您没有）或周期性波动（您也没有）。滑动窗口只是对连续的 bin 取平均值，平均随机波动。从数学上讲，这是一个卷积。

从技术上讲，您实际上可以像这样解决您的问题（自己尝试更大的 Nwindow 值）：

Nwindow=10
y_slide_10 = np.convolve(y,np.ones((Nwindow,))/Nwindow,mode='same')
Nwindow=20
y_slide_20 = np.convolve(y,mode='same')
Nwindow=30
y_slide_30 = np.convolve(y,mode='same')

plt.xlabel('time')
plt.ylabel('signal')
plt.plot(x,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()

因此，从技术上讲，您可以成功抑制最初的“驼峰”。但不要忘记这是一个手动调整而不是通用的解决方案...

任何滑动窗口解决方案的

另一个警告：这总是会扭曲您的时间安排。由于您根据上升或下降信号对时间间隔进行平均，因此您错综复杂的轨迹会在时间上来回移动（轻微但显着）。在您的特定情况下，这不是问题，因为主要信号区域基本上没有时间依赖性（非常平坦）。

频域

这应该是灵丹妙药，但它也不适用于您的示例。这不能更好地工作的事实是对我的主要暗示，最好丢弃前 40 秒的数据......（即在科学工作中）

您可以使用快速傅立叶变换来检查频域中的数据。

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 秒后的稳定区域，驼峰与时间的（快速）变化相关。您现在可以玩转并更改频谱（--> 过滤器），但我认为频谱太人为，因此在这里不会产生很好的结果。与其他数据一起尝试，您可能会印象深刻！我尝试了两件事，首先切掉高频区域（设置为零），其次，在频域中应用滑动窗口滤波器（将峰值保留在 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:],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()

将其转换回时域：

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

收益

这些解决方案中存在明显的波动，这使得它们对于您的目的来说基本上毫无用处。这导致我进行最后的练习，再次在“频率滑动窗口”时域上应用滑动窗口滤波器

# extra time-domain sliding window
Nwindow=90
y_fft_90 = np.convolve(y_filtered_slide,mode='same')

# final time-domain plot
plt.plot(x[:500],y_fft_90[:500],label="frequency-sliding window,slide")
# plt.xscale('log') # useful
plt.legend()
plt.show()

我对这个结果很满意，但它仍然有很小的振荡，因此不能解决您原来的问题。

结论

多么有趣。浪费了一个小时。也许它对某人有用。甚至对你来说，娜塔莎。请不要生我的气...

,

假设您的数据在 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。

,

如果您正在处理最终单调收敛到某个固定值的确定性数据，则问题非常简单。您的最后一次观察结果应该最接近极限，因此您可以相对于最后一个数据点定义一个可接受的容差阈值，并从后向前扫描您的数据，以找出超出阈值的位置。

一旦将随机噪声添加到图片中，事情就会变得更加糟糕，尤其是在存在序列相关的情况下。这个问题在仿真建模中很常见^{（见下面的 (*)）}，被称为 initial bias 问题。它首先由 Conway in 1963 确定，从那时起一直是一个活跃的研究领域，关于如何处理它没有普遍接受的明确答案。与确定性情况一样，最广泛接受的答案从数据集的右侧开始解决问题，因为这是数据最有可能处于稳定状态的地方。基于这种方法的技术使用数据集的末尾来建立某种统计标准或基线，以测量数据开始看起来显着不同的位置，因为通过向数据集的前面移动而增加了观察。由于存在序列相关性，这会变得非常复杂。

如果时间序列处于稳态，从 covariance stationary 的意义上说，那么数据的简单平均值是对其期望值的无偏估计，但估计均值的标准误差在很大程度上取决于序列相关。正确的标准误差平方不再是 s²/n，而是 (s²/n)*W，其中 W 是自相关值的适当加权和. 1990 年代开发了一种称为 MSER 的方法，它通过尝试确定标准误差在何处最小化来避免尝试正确估计 W 的问题。如果样本量足够大，它将 W 视为事实上的常数，因此，如果您考虑两个标准误差估计值的比率，W 会抵消，并且最小值出现在 s²/n 最小的地方。 MSER 进行如下：

• 从最后开始，计算一半数据集的 s² 以建立基线。
• 现在使用一种有效的技术，例如 Welford's online algorithm，一次更新一次观测值 s² 的估计值，计算 s²/n，其中 n 是到目前为止统计的观察数量。跟踪哪个 n 值产生最小的 s²/n。起泡、冲洗、重复。
• 一旦你从后到前遍历了整个数据集，产生最小 s²/n 的 n 是来自数据集末尾的不可检测的观察数因为受到起始条件的影响。

理由 - 具有足够大的基线（您的数据的一半），只要时间序列保持稳定状态，s²/n 应该相对稳定。由于 n 是单调递增的，因此 s²/n 应继续递减，但受其作为估计值的可变性的限制。但是，一旦您开始获取非稳态观测值，均值和方差的漂移将使 s²/n 的分子膨胀。因此，最小值对应于没有非平稳性迹象的最后一次观察。可以在此 proceedings paper 中找到更多详细信息。一个Ruby implementation is available on BitBucket。

您的数据的变化量如此之小，以至于 MSER 得出结论，它仍在收敛到稳定状态。因此，我建议采用第一段中概述的确定性方法。如果您将来有嘈杂的数据，我绝对建议您尝试 MSER。

---

(*) - 简而言之，模拟模型是一个计算机程序，因此必须将其状态设置为一些初始值。从长远来看，我们通常不知道系统状态会是什么样子，因此我们将其初始化为一组任意但方便的值，然后让系统“预热”。问题是模拟的初始结果不是典型的稳态行为，因此在分析中包含这些数据会使它们产生偏差。解决方案是去除数据的偏差部分，但应该是多少？