Numpy Root-Mean-Squared (RMS) сглаживания сигнала - PullRequest
Новая
       15

Numpy Root-Mean-Squared (RMS) сглаживания сигнала

8 голосов
/ 23 ноября 2011

У меня есть сигнал электромиографических данных, которые я должен (явная рекомендация научных трудов) сгладить с помощью 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, но я не мог понять, как добиться того, что я хочу, используя их.

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

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

Ответы [ 3 ]

12 голосов
/ 24 ноября 2011

В можно использовать свертку для выполнения операции, на которую вы ссылаетесь. Я делал это несколько раз и для обработки сигналов ЭЭГ. 

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' (см. документацию ).

1 голос
/ 24 ноября 2011

Поскольку это не линейное преобразование, я не верю, что можно использовать 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])
0 голосов
/ 16 июля 2019

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

Быстрое вычисление окна RMS с быстрым вычислением

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

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

Другими словами, если мы имеем

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

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

Давайте получим acummulator 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)

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

...