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ますが、各値は 2 乗されます。各要素が でnp.ones(window_size)/float(window_size)ある配列または長さを生成します。したがって、畳み込みは、各要素が等しい新しい配列を効果的に生成しますwindow_size1/window_sizei

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

これは、移動ウィンドウ内の配列要素の RMS 値です。このように非常にうまく機能するはずです。

ただし、同じ次元の新しい配列がnp.power(a, 2)生成されることに注意してください。が本当に大きい場合、つまりメモリに 2 回収まらないほど大きい場合、各要素をその場で変更する戦略が必要になる場合があります。また、引数は境界効果を破棄するように指定しているため、 によって生成される配列が小さくなります。代わりに指定することで、すべてを保持できます(ドキュメントを参照)。a'valid'np.convolve()'same'

于 2011-11-24T16:49:15.893 に答える
1

マシンが畳み込みに苦労していることがわかったので、次の解決策を提案します。

移動 RMS ウィンドウの迅速な計算

アナログ電圧サンプル a0 ～ a99 (100 サンプル) があり、それらを介して 10 サンプルの移動 RMS を取得する必要があるとします。

ウィンドウは最初に要素 a0 から a9 (10 サンプル) をスキャンして 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 を作成するには、10a^2の足し算の sqrt を取り、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 メガ サンプルの配列を計算できます。

于 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 に答える