我具有表示数字输出CSV值的数组。 已使用模拟示波器,所以它不是一个完美的数字信号采集。 我试图滤出数据有一个完美的数字信号用于计算周期(其可能会发生变化)。 我也想定义的最大的错误,我从这个过滤得到。
事情是这样的:
理念
套用treshold OD数据。 下面是一个伪代码:
for data_point_raw in data_array:
if data_point_raw < 0.8: data_point_perfect = LOW
if data_point_raw > 2 : data_point_perfect = HIGH
else:
#area between thresholds
if previous_data_point_perfect == Low : data_point_perfect = LOW
if previous_data_point_perfect == HIGH: data_point_perfect = HIGH
有两个问题困扰着我。
- 这似乎是在数字信号处理中的常见问题,但是我还没有为它找到一个预定义的标准功能。 这是执行过滤的方式OK?
- 我将如何得到最大的错误?
下面是一些代码,这可能有助于。
from __future__ import division
import numpy as np
def find_transition_times(t, y, threshold):
"""
Given the input signal `y` with samples at times `t`,
find the times where `y` increases through the value `threshold`.
`t` and `y` must be 1-D numpy arrays.
Linear interpolation is used to estimate the time `t` between
samples at which the transitions occur.
"""
# Find where y crosses the threshold (increasing).
lower = y < threshold
higher = y >= threshold
transition_indices = np.where(lower[:-1] & higher[1:])[0]
# Linearly interpolate the time values where the transition occurs.
t0 = t[transition_indices]
t1 = t[transition_indices + 1]
y0 = y[transition_indices]
y1 = y[transition_indices + 1]
slope = (y1 - y0) / (t1 - t0)
transition_times = t0 + (threshold - y0) / slope
return transition_times
def periods(t, y, threshold):
"""
Given the input signal `y` with samples at times `t`,
find the time periods between the times at which the
signal `y` increases through the value `threshold`.
`t` and `y` must be 1-D numpy arrays.
"""
transition_times = find_transition_times(t, y, threshold)
deltas = np.diff(transition_times)
return deltas
if __name__ == "__main__":
import matplotlib.pyplot as plt
# Time samples
t = np.linspace(0, 50, 501)
# Use a noisy time to generate a noisy y.
tn = t + 0.05 * np.random.rand(t.size)
y = 0.6 * ( 1 + np.sin(tn) + (1./3) * np.sin(3*tn) + (1./5) * np.sin(5*tn) +
(1./7) * np.sin(7*tn) + (1./9) * np.sin(9*tn))
threshold = 0.5
deltas = periods(t, y, threshold)
print("Measured periods at threshold %g:" % threshold)
print(deltas)
print("Min: %.5g" % deltas.min())
print("Max: %.5g" % deltas.max())
print("Mean: %.5g" % deltas.mean())
print("Std dev: %.5g" % deltas.std())
trans_times = find_transition_times(t, y, threshold)
plt.plot(t, y)
plt.plot(trans_times, threshold * np.ones_like(trans_times), 'ro-')
plt.show()
输出:
Measured periods at threshold 0.5:
[ 6.29283207 6.29118893 6.27425846 6.29580066 6.28310224 6.30335003]
Min: 6.2743
Max: 6.3034
Mean: 6.2901
Std dev: 0.0092793
可以使用numpy.histogram
和/或matplotlib.pyplot.hist
进一步分析由返回的数组periods(t, y, threshold)
。
这不是你的问题的答案,只是与建议,可以帮助。 我在这里写它,因为我不能把图像中的评论。
我想你应该处理之前不知何故数据标准化。
正常化的0范围之后...... 1你应该申请您的过滤器。
如果你真的只在月经感兴趣,你可以绘制傅立叶变换,你就必须在信号的频率出现的峰值(等你的时间)。 在傅里叶域中的宽峰,在周期测量中的误差较大
import numpy as np
data = np.asarray(my_data)
np.fft.fft(data)
你的过滤是好的,这是基本相同的施密特触发器,但你可能有它的主要问题是速度。 使用NumPy的的好处在于,它可以为C一样快,而你有超过每个元素重复一次。
您可以使用从SciPy的中值滤波器实现类似的东西。 所述应以下实现了类似的结果(而不是依赖于任何量值):
filtered = scipy.signal.medfilt(raw)
filtered = numpy.where(filtered > numpy.mean(filtered), 1, 0)
可以调中值与滤波的强度medfilt(raw, n_samples)
, n_samples
缺省值为3。
至于错误,那将是非常主观的。 一种方法是不过滤,discretise信号,然后比较不同之处。 例如:
discrete = numpy.where(raw > numpy.mean(raw), 1, 0)
errors = np.count_nonzero(filtered != discrete)
error_rate = errors / len(discrete)