Search code examples
pythonnumpyenvelope

Why is the envelope curve erroneous at the beginning?


I implemented a function that gives the envelope curve of discrete values. I think there might be an error as when I tested it over a date I will provide you with at the bottom of the post, I get dissimilarities between the real data points ant the envelope curve as sown in this figure example of the code applied to a given data

from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt

def enveloppe(s):
    u_x = [0,]        
    u_y = [s[0],]
    q_u = np.zeros(s.shape)
    for k in xrange(1,len(s)-1):
        if (np.sign(s[k]-s[k-1])==1) and (np.sign(s[k]-s[k+1])==1):
            u_x.append(k)
            u_y.append(s[k])
    u_x.append(len(s)-1)
    u_y.append(s[-1]) 
    u_p = interp1d(u_x,u_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)
    return q_u   

fig, ax = plt.subplots()
ax.plot(S, '-o', label = 'magnitude')
ax.plot(envelope(S), '-o', label = 'enveloppe magnitude')
ax.legend()

Data S : array([  9.12348621e-11,   6.69568658e-10,   6.55973768e-09,
         1.26822485e-06,   4.50553316e-09,   5.06526113e-07,
         2.96728433e-09,   2.36088205e-07,   1.90802318e-09,
         1.15867354e-07,   1.18504790e-09,   5.72888034e-08,
         6.98672478e-10,   2.75361324e-08,   3.82391643e-10,
         1.25393143e-08,   1.96697343e-10,   5.96979943e-09,
         1.27009013e-10,   4.46365555e-09,   1.31769958e-10,
         4.42024233e-09,   1.42514400e-10,   4.17757107e-09,
         1.41640360e-10,   3.65170558e-09,   1.29784598e-10,
         2.99790514e-09,   1.11732461e-10])

Solution

  • I would make two modifications to your enveloppe function to get a more monotone output

    enter image description here

    The idea is to avoid implicit addition of the left and right ends to the list of peaks that is used to construct the envelope

    def enveloppe(s):
        u_x = [] # do not add 0
        u_y = []
        q_u = np.zeros(s.shape)
        for k in range(1,len(s)-1):
            if (np.sign(s[k]-s[k-1])==1) and (np.sign(s[k]-s[k+1])==1):
                u_x.append(k)
                u_y.append(s[k])
        print(u_x)
        u_p = interp1d(u_x,u_y, kind = 'cubic',
                  bounds_error = False, 
                  fill_value="extrapolate") # use fill_value="extrapolate"
        for k in range(0,len(s)):
            q_u[k] = u_p(k)
        return q_u