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

Question

我想在带有噪声的时间序列数据中找到达到某个值的时刻。如果数据中没有峰值，我可以在 MATLAB 中执行以下操作。

代码从这里

% 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

但是当有噪音时，我不确定必须做什么。

In the time-series data plotted in the above image I want to find the time instant at which the y-axis value 5 is reached. The data actually stabilizes to 5 at t>=100 s. But due to the presence of noise in the data, we see a peak that reaches 5 somewhere around 20 s . I would like to know how to detect e.g 100 seconds as the right time and not 20 s . The code posted above will only give 20 s as the answer. I saw a post here that explains using a sliding window to find when the data equilibrates. However, I am not sure how to implement the same. Suggestions will be really helpful.

The sample data plotted in the above image can be found here

Suggestions on how to implement in Python or MATLAB code will be really helpful.

EDIT: I don't want to capture when the peak (/noise/overshoot) occurs. I want to find the time when equilibrium is reached. For example, around 20 s the curve rises and dips below 5. After ~100 s the curve equilibrates to a steady-state value 5 and never dips or peaks.

Answer 1

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

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

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

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

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

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

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

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

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

频域

这应该是灵丹妙药，但它也不适用于您的示例。这不能更好地工作的事实是对我的主要暗示，最好丢弃前 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:], 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()

将其转换回时域：

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

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

结论

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