13

我有一个肌电图数据信号，我应该（科学论文的明确建议）使用 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和建议，也来自 numpy，但我无法弄清楚如何使用它们来完成我想要的。itertoolsconvolve

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

谢谢阅读

4

3 回答 3

17

可以使用卷积来执行您所指的操作。我也做了几次处理脑电图信号。

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'（请参阅文档）。

于 2011-11-24T16:49:15.893 回答
1

我发现我的机器在 convolve 上苦苦挣扎，所以我提出了以下解决方案：

快速计算移动 RMS 窗口

假设我们有模拟电压样本 a0 ... a99（一百个样本），我们需要通过它们获取 10 个样本的移动 RMS。

该窗口最初将从元素 a0 扫描到 a9（十个样本）以获得 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
    ...

简化它：我们必须a = [a0, ... a99] 创建 10 个样本的移动 RMS，我们可以取 sqrt 的加法 10a^2并乘以 1/10。

换句话说，如果我们有

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

rms^2只需添加一组 10 个 p 。

让我们有一个累加器 acu：

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

然后我们可以有

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

我们可以创建：

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

    V3 = V0 + V1 + V2

如果我们运行 itertools.accumulate(V3) ，我们将得到 rms 数组

代码：

    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)

我可以在不到 1 分钟的时间内计算出包含 128 个 Mega 样本的数组。

于 2019-07-16T05:34:03.317 回答
0

由于这不是线性变换，我不相信可以使用 np.convolve()。

这是一个应该做你想做的功能。请注意，返回数组的第一个元素是第一个完整窗口的 rms；即对于a示例中的数组，返回数组是子窗口的 rms，[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])
于 2011-11-23T21:11:41.143 回答