Python - How to get high and low envelope of a signal?

Question

I have pretty noisy data and am trying to develop a high and low envelope signal. This is similar to this example in MATLAB:

http://uk.mathworks.com/help/signal/examples/signal-smoothing.html

in the "Removing the Peak Envelope" section. Is there a similar function in Python that can do this? My whole project is written in Python, in the worst case I can extract my numpy array and throw it in MATLAB and use this example. But I prefer the look of matplotlib ... and really cba does all these I / O between MATLAB and Python ...

Thanks,

Answer 1

Is there a similar function in Python that can do this?

As far as I know, there is no such function in Numpy / Scipy / Python. However, creating it is not so difficult. The general idea is this:

For the vector of values ​​(s):

• Find the location of the peaks (s). Let us call them (u)
• Find the location of the s depressions. Let us call them (l).
• Set the model to value pairs (u). Let me call him (u_p)
• Set the model to value pairs (1). Let me call him (l_p)
• Rate (u_p) over domain (s) to get interpolated upper shell values. (Let me call them (q_u))
• Rate (l_p) in domain (s) to get the interpolated values ​​of the lower envelope. (Let me call them (q_l)).

As you can see, this is a sequence of three steps (find location, fit the model, evaluate the model), but it is applied twice, once for the upper part of the envelope and the other for the lower.

To collect the “peaks” (s), you need to find the points where the slope changes from positive to negative, and to collect the “peaks” (s) you need to find the points where the slope (s) changes from negative to positive.

Peak example: s = [4,5,4] 5-4 positive 4-5 negative

Gutter example: s = [5,4,5] 4-5 negative 5-4 positive

Here is an example script to get you started with a lot of inline comments:

from numpy import array, sign, zeros from scipy.interpolate import interp1d from matplotlib.pyplot import plot,show,hold,grid s = array([1,4,3,5,3,2,4,3,4,5,4,3,2,5,6,7,8,7,8]) #This is your noisy vector of values. q_u = zeros(s.shape) q_l = zeros(s.shape) #Prepend the first value of (s) to the interpolating values. This forces the model to use the same starting point for both the upper and lower envelope models. u_x = [0,] u_y = [s[0],] l_x = [0,] l_y = [s[0],] #Detect peaks and troughs and mark their location in u_x,u_y,l_x,l_y respectively. for k in xrange(1,len(s)-1): if (sign(s[k]-s[k-1])==1) and (sign(s[k]-s[k+1])==1): u_x.append(k) u_y.append(s[k]) if (sign(s[k]-s[k-1])==-1) and ((sign(s[k]-s[k+1]))==-1): l_x.append(k) l_y.append(s[k]) #Append the last value of (s) to the interpolating values. This forces the model to use the same ending point for both the upper and lower envelope models. u_x.append(len(s)-1) u_y.append(s[-1]) l_x.append(len(s)-1) l_y.append(s[-1]) #Fit suitable models to the data. Here I am using cubic splines, similarly to the MATLAB example given in the question. u_p = interp1d(u_x,u_y, kind = 'cubic',bounds_error = False, fill_value=0.0) l_p = interp1d(l_x,l_y,kind = 'cubic',bounds_error = False, fill_value=0.0) #Evaluate each model over the domain of (s) for k in xrange(0,len(s)): q_u[k] = u_p(k) q_l[k] = l_p(k) #Plot everything plot(s);hold(True);plot(q_u,'r');plot(q_l,'g');grid(True);show() 

This produces this conclusion:

Points for further improvement:

• The code above does not filter peaks or troughs that may occur closer than some threshold “distance” (Tl) (for example, time). This is similar to the second envelope parameter. It is easy to add if you analyze the differences between successive values u_x,u_y .

• However, a quick improvement over the aforementioned point is a low-pass filter for your data with a moving average filter BEFORE interpolating the upper and lower envelope functions. You can do this easily by folding yours with a suitable moving average filter. Without going into details (you can do it if necessary) to create a moving average filter that works on N consecutive patterns, you would do something like this: s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N) . The higher (N) the more smooth the information will be displayed. Note that this will shift your (N / 2) values ​​to the right (in s_filtered ) due to what is called group delay of the smoothing filter. For more information on the moving average, see this link .

Hope this helps.

(If you want more information about the original application, you can compensate for the answer. Perhaps the data can be pre-processed in a more appropriate way (?))