RMS RMS signal smoothing

Question

I have an electromyographic data signal that I suppose (explicit recommendations from scientific articles) to smooth out the use of RMS.

I have the following working code creating the desired result, but it is slower than I think is possible.

#!/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 

I saw some deque and itertools regarding optimizing hinges of movable windows as well as convolve from numpy, but I could not figure out how to accomplish what I want to use.

Also, I don’t want to avoid border problems anymore, because in the end I have large arrays and relatively small sliding windows.

Thanks for reading

Answer 1

You can use convolution to perform the operation you are referring to. I have done this several times to process EEG signals.

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

By breaking it, the np.power(a, 2) creates a new array with the same size as a , but where each value is quadratic. np.ones(window_size)/float(window_size) creates an array or the length of window_size , where each element is 1/window_size . Thus, convolution effectively creates a new array, where each element i is equal to

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

which is the RMS value of the array elements in the moving window. That should work just fine.

Note that np.power(a, 2) creates a new array of the same size. If a has a size really big, I mean big enough so that it cannot fit twice in memory, you may need a strategy in which each element is modified in place. In addition, the argument 'valid' indicates a reset of border effects, resulting in a smaller array created by np.convolve() . You can save all this by specifying 'same' instead (see documentation ).

Answer 2

Since this is not a linear transformation, I do not believe that np.convolve () can be used.

Here is a function that should do what you want. Note that the first element of the returned array is the rms of the first full window; those. for array a in the example, the returned array is the rms value of the subframes [1,2],[2,3],[3,4],[4,5] and does not include the partial windows [1] and [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]) 

Answer 3

I found that my machine is struggling with convolution, so I offer the following solution:

Fast RMS moving window calculation

Suppose we have analog voltage samples a0 ... a99 (one hundred samples), and we need to go through RMS 10 samples.

Initially, the window will be scanned from elements a0 through 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 ... 

Simplification: we have a = [a0,... a99] To create an RMS of 10 samples, we can take sqrt of addition 10 a^2 and multiply by 1/10.

In other words, if we have

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

To get rms^2 just add a group of 10 people.

Let there be 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 get an rms array

The 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 calculate an array of 128 megapixels in less than 1 minute.