Numpy Root-Mean-Squared (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 предложения по оптимизации движущихся оконных петель, а также convolve от NumPy, но я не мог понять, как добиться того, что я хочу, используя их.

Кроме того, я больше не хочу избегать проблем с границами, потому что у меня большие массивы и относительно маленькие раздвижные окна.

Спасибо за прочтение

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

которое является среднеквадратичным значением элементов массива в движущемся окне. Так должно быть очень хорошо.

Обратите внимание, что np.power(a, 2) создает новый массив того же размера. Если a действительно большой, я имею в виду достаточно большой, чтобы он не мог поместиться дважды в памяти, вам может понадобиться стратегия, в которой каждый элемент изменяется на месте. Так же 'valid' Аргумент указывает, что нужно отбрасывать эффекты границы, в результате чего получается меньший массив np.convolve(), Вы можете сохранить все это, указав 'same' вместо этого (см. документацию).

Answer 2

Я обнаружил, что моя машина борется со сверткой, поэтому предлагаю следующее решение:

Быстрое вычисление движущегося окна RMS

Предположим, у нас есть выборки аналогового напряжения a0 ... a99 (сто выборок), и нам нужно пройти через них скользящее среднеквадратичное значение 10 выборок.

Окно будет сначала сканировать от элементов a0 до a9 (десять отсчетов), чтобы получить среднеквадратичное значение 0.

    # 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 выборок, мы можем взять sqrt из добавления 10 a^2и умножить на 1/10.

Другими словами, если у нас есть

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

Получить rms^2 просто добавьте группу из 10 п.

Имеем акуммулятор аку:

    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)

Я могу вычислить массив из 128 мегасэмплов менее чем за 1 минуту.

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])