Numpy均方根（RMS）信号平滑处理

Question

我有一个肌电数据信号，根据科学论文的明确建议，我应该使用RMS平滑处理。

我有以下工作代码，可以产生所需的输出，但速度比我认为可能的要慢得多。

#!/usr/bin/python
import numpy
def rms(interval, halfwindow):
    """ performs the moving-window smoothing of a signal using RMS """
    n = len(interval)
    rms_signal = numpy.zeros(n)
    for i in range(n):
        small_index = max(0, i - halfwindow)  # intended to avoid boundary effect
        big_index = min(n, i + halfwindow)    # intended to avoid boundary effect
        window_samples = interval[small_index:big_index]

        # here is the RMS of the window, being attributed to rms_signal 'i'th sample:
        rms_signal[i] = sqrt(sum([s**2 for s in window_samples])/len(window_samples))

    return rms_signal

我看到一些有关优化移动窗口循环的建议，包括使用deque和itertools，以及来自numpy的convolve。但是我不知道如何使用它们来完成我想要的操作。

此外，我不再关心边界问题，因为我最终会得到大数组和相对较小的滑动窗口。

谢谢阅读。

- heltonbiker

---

你能提供这篇论文的链接吗？我从未听说过通过计算移动窗口上点的RMS来平滑信号。一般来说，这不会看起来像原始信号的平滑版本。 - user545424

建议使用这种平滑方法，因为它与信号功率（能量）相关，这可以用来推断肌肉的努力。链接：http://www.isek-online.org/standards_emg.html “提供幅度信息的另一种可接受方法是“均方根”或RMS。就像移动平均值一样，这个数量被定义为一个特定的时间间隔（移动窗口）T，必须指出。” 根据Noraxon手册（由我的公司拥有的闭源）的说法，这是平滑的首选，时间窗口大约在50到100毫秒之间。 - heltonbiker

移动窗口的均方根是音频电平表背后的理念。 - endolith

Answer 1

使用卷积来执行您所提到的操作是可能的。我之前也用它处理脑电图信号。

import numpy as np
def window_rms(a, window_size):
  a2 = np.power(a,2)
  window = np.ones(window_size)/float(window_size)
  return np.sqrt(np.convolve(a2, window, 'valid'))

通过分解，np.power(a, 2) 部分生成一个与 a 相同维度的新数组，但每个值都被平方。 np.ones(window_size)/float(window_size) 生成一个长度为 window_size 的数组，其中每个元素都是 1/window_size。因此，卷积有效地生成一个新数组，其中每个元素 i 等于

(a[i]^2 + a[i+1]^2 + … + a[i+window_size]^2)/window_size

这是在移动窗口内计算数组元素RMS值。这种方法的性能非常好。

需要注意的是，np.power(a, 2)会产生一个与输入数组同维度的新数组。如果a非常大，即使只是一次在内存中存储不下，你可能需要采用逐个修改每个元素的策略。另外，'valid'参数指定舍弃边缘效应，从而由np.convolve()生成一个较小的数组。你可以通过选择'same'来保留所有内容（详见文档）。

- matehat

---

太棒了，非常聪明！EMG信号比EMG短得多，4通道录制不到50mb，因此内存使用不会受到限制。我计划使用“same”来在堆栈中显示原始/平滑/滤波信号。 - heltonbiker

仅供记录， RMS 现在快了500倍（使用 cProfile 测量，对于 16 个样本的 window_size，时间从0.002秒变为1.000秒），而且即使是更大的窗口（100+）也没有明显变慢。我说过它很棒吗？ - heltonbiker

你不应该在代码中指定除“valid”以外的任何内容，因为np.convolve()的第二个参数包含完整窗口长度的倒数。 - user545424

只是稍微扩展一下，如果需要一些更奇特的行为，可以将“window”作为核心，其总和为“1.0”，例如归一化高斯核。实际上，函数中的三行代码执行了一些DSP文本通常称为“去线性化”，“解调”和“重新线性化”的操作，这些操作可以使用不同的幂次（除了二次方），核函数（除了单位平方或高斯核），统计算子（除了加权平均）和窗口大小来完成。 - heltonbiker

np.convolve(a2, window, 'valid') 应该替换为 np.convolve(a2, window_size, 'valid') 吗？看起来是在 np.convolve 中的 window 打错了。 - Matias Andina

Answer 2

我发现我的机器在卷积方面遇到了困难，因此我提出以下解决方案：

快速计算移动 RMS 窗口

Suppose we have analog voltage samples a0 ... a99 (one hundred samples) and we need to take moving RMS of 10 samples through them.

The window will scan initially from elements a0 to a9 (ten samples) to get rms0.

    # rms = [rms0, rms1, ... rms99-9] (total of 91 elements in list):
    (rms0)^2 = (1/10) (a0^2 + ...         + a9^2)            # --- (note 1)
    (rms1)^2 = (1/10) (...    a1^2 + ...  + a9^2 + a10^2)    # window moved a step, a0 falls out, a10 comes in
    (rms2)^2 = (1/10) (              a2^2 + ... + a10^2 + a11^2)     # window moved another step, a1 falls out, a11 comes in
    ...

Simplifying it: We havea = [a0, ... a99] To create moving RMS of 10 samples, we can take sqrt of the addition of 10 a^2's and multiplied by 1/10.

In other words, if we have

    p = (1/10) * a^2 = 1/10 * [a0^2, ... a99^2]

To get rms^2 simply add a group of 10 p's.

Let's have an acummulator acu:

    acu = p0 + ... p8     # (as in note 1 above)

Then we can have

    rms0^2 =  p0 + ...  p8 + p9 
           = acu + p9
    rms1^2 = acu + p9 + p10 - p0
    rms2^2 = acu + p9 + p10 + p11 - p0 - p1
    ...

we can create:

    V0 = [acu,   0,   0, ...  0]
    V1 = [ p9, p10, p11, .... p99]          -- len=91
    V2 = [  0, -p0, -p1, ... -p89]          -- len=91

    V3 = V0 + V1 + V2

if we run itertools.accumulate(V3) we will get rms array

Code:

    import numpy as np
    from   itertools import accumulate

    a2 = np.power(in_ch, 2) / tm_w                  # create array of p, in_ch is samples, tm_w is window length
    v1 = np.array(a2[tm_w - 1 : ])                  # v1 = [p9, p10, ...]
    v2 = np.append([0], a2[0 : len(a2) - tm_w])     # v2 = [0,   p0, ...]
    acu = list(accumulate(a2[0 : tm_w - 1]))        # get initial accumulation (acu) of the window - 1
    v1[0] = v1[0] + acu[-1]                         # rms element #1 will be at end of window and contains the accumulation
    rmspw2 = list(accumulate(v1 - v2))

    rms = np.power(rmspw2, 0.5)

I can compute an array of 128 Mega samples in less than 1 minute.

- Ketut Wiadnyana

Answer 3

由于这不是线性变换，我认为不可能使用np.convolve()。

这里有一个函数可以实现你想要的功能。请注意，返回数组的第一个元素是第一个完整窗口的均方根；即对于示例中的数组a，返回数组是子窗口[1,2]，[2,3]，[3,4]，[4,5]的均方根，不包括部分窗口[1]和[5]。

>>> def window_rms(a, window_size=2):
>>>     return np.sqrt(sum([a[window_size-i-1:len(a)-i]**2 for i in range(window_size-1)])/window_size)
>>> a = np.array([1,2,3,4,5])
>>> window_rms(a)
array([ 1.41421356,  2.44948974,  3.46410162,  4.47213595])

- user545424

---

可以使用卷积，详见我的回答。 - matehat