将模拟信号数字化

Question

我有一个CSV值数组，表示数字输出。它是使用模拟示波器收集的，因此不是完美的数字信号。我正在尝试过滤数据，以获得完美的数字信号，以计算周期（可能会变化）。

我还想定义从此过滤中获得的最大误差。

类似这样：

想法

对数据应用阈值。以下是伪代码：

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

有两个问题困扰着我。

1. 这似乎是数字信号处理中的常见问题，但我没有找到预定义的标准函数。这种过滤的方法可行吗？
2. 我该如何获得最大误差？

- TheMeaningfulEngineer

---

完全离题，但模拟示波器在工作台上比数字示波器要好得多。数字示波器具有最小粒度（输入和显示输出都是如此），并且往往会忽略异常值。 - L0j1k

你考虑过将 Matlab 与 Python 结合使用吗？可以参考这个链接：Python vs Matlab 和 Matplotlib。 - Inbar Rose

@L0j1k 嗯，实际上它是数字的，我们只是使用模拟探头，这样我们就可以自己进行转换 :) - TheMeaningfulEngineer

@InbarRose 抱歉，Matlab不在考虑范围内。不知道你为什么推荐了matplotlib的链接？问题中的图是由它生成的，并在Libre Office Draw中进行了一些编辑。 - TheMeaningfulEngineer

Answer 1

这里有一段可能有帮助的代码。

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) 返回的数组。

- Warren Weckesser

Answer 2

这不是针对你的问题的答案，只是一个可能有所帮助的建议。我在这里写下它，因为无法在评论中添加图片。

我认为您应该在任何处理之前规范化数据。

将数据规范化到0...1范围后，应用过滤器。

- Kamil

---

谢谢你的建议。不过，在这个具体的例子中，我将信号转换为布尔值，所以我看不出归一化的好处。 - TheMeaningfulEngineer

Answer 3

如果你只对周期感兴趣，你可以绘制傅里叶变换，你会看到一个峰值，代表信号的频率（也就是周期）。在傅里叶域中峰值越宽，你的周期测量误差就越大。

import numpy as np

data = np.asarray(my_data)

np.fft.fft(data)

- danodonovan

---

很酷的想法，但我们希望获得一个时间段的直方图，以便我们可以知道它在多大程度上适合我们的需求。这是在树莓派上实现不同软件PWM的领域中进行实验的。 - TheMeaningfulEngineer

Answer 4

你的过滤器很好，基本上与Schmitt触发器相同，但可能会遇到的主要问题是速度。使用Numpy的好处在于它可以像C一样快，而你必须对每个元素进行一次迭代。

你可以使用SciPy的中值过滤器实现类似的效果。以下内容应该可以实现类似的结果（并且不依赖于任何量级）：

filtered = scipy.signal.medfilt(raw)
filtered = numpy.where(filtered > numpy.mean(filtered), 1, 0)

你可以使用medfilt(raw, n_samples)来调整中值滤波器的强度，n_samples默认为3。

至于错误，这将是非常主观的。一种方法是在不进行过滤的情况下离散化信号，然后比较差异。例如：

discrete = numpy.where(raw > numpy.mean(raw), 1, 0)
errors = np.count_nonzero(filtered != discrete)
error_rate = errors / len(discrete)

- jleahy