第一次尝试是利用scipy Hilbert transform来确定幅度包络，但在许多情况下并不符合预期，主要原因如下：
希尔伯特包络，也称为能量-时间曲线（ETC），只适用于窄带波动。产生一个分析信号，然后取其绝对值，是一个线性运算，所以它平等地对待信号的所有频率。如果你给它一个纯正弦波，它确实会返回给你一条直线。但是，当你给它白噪声，你可能会得到噪声回来。
从那时起，由于其他的答案都是使用三次样条插值法，而且往往会变得很麻烦，有点不稳定（寄生振荡），而且对于非常长和嘈杂的数据阵列来说很耗时，我将在这里提供一个简单而有效的版本，看起来工作得很好：

import numpy as np
from matplotlib import pyplot as plt

def hl_envelopes_idx(s, dmin=1, dmax=1, split=False):
    """
    Input :
    s: 1d-array, data signal from which to extract high and low envelopes
    dmin, dmax: int, optional, size of chunks, use this if the size of the input signal is too big
    split: bool, optional, if True, split the signal in half along its mean, might help to generate the envelope in some cases
    Output :
    lmin,lmax : high/low envelope idx of input signal s
    """

    # locals min      
    lmin = (np.diff(np.sign(np.diff(s))) > 0).nonzero()[0] + 1 
    # locals max
    lmax = (np.diff(np.sign(np.diff(s))) < 0).nonzero()[0] + 1 
    

    if split:
        # s_mid is zero if s centered around x-axis or more generally mean of signal
        s_mid = np.mean(s) 
        # pre-sorting of locals min based on relative position with respect to s_mid 
        lmin = lmin[s[lmin]<s_mid]
        # pre-sorting of local max based on relative position with respect to s_mid 
        lmax = lmax[s[lmax]>s_mid]

    # global max of dmax-chunks of locals max 
    lmin = lmin[[i+np.argmin(s[lmin[i:i+dmin]]) for i in range(0,len(lmin),dmin)]]
    # global min of dmin-chunks of locals min 
    lmax = lmax[[i+np.argmax(s[lmax[i:i+dmax]]) for i in range(0,len(lmax),dmax)]]
    
    return lmin,lmax

• • 示例1：* 准周期振动***

t = np.linspace(0,8*np.pi,5000)
s = 0.8*np.cos(t)**3 + 0.5*np.sin(np.exp(1)*t)
lmin, lmax = hl_envelopes_idx(s)

# plot
plt.plot(t,s,label='signal')
plt.plot(t[lmin], s[lmin], 'r', label='low')
plt.plot(t[lmax], s[lmax], 'g', label='high')

• • 示例2：* 噪声衰减信号***

t = np.linspace(0,2*np.pi,5000)
s = 5*np.cos(5*t)*np.exp(-t) + np.random.rand(len(t))

lmin, lmax = hl_envelopes_idx(s,dmin=15,dmax=15)

# plot
plt.plot(t,s,label='signal')
plt.plot(t[lmin], s[lmin], 'r', label='low')
plt.plot(t[lmax], s[lmax], 'g', label='high')

• • 示例3：* 非对称调制啁啾***

18867925样本的复杂得多的信号（此处不包括）：

Python中是否有类似的函数可以做到这一点？
据我所知，Numpy/Scipy/Python中没有这样的函数，但是创建一个并不困难，大致思路如下：
给定一个值向量：
1.求（s）的峰的位置，我们称之为（u）
1.求s的波谷位置，记为（l）。
1.将模型拟合到（u）值对，我们称之为（u_p）
1.对（l）个值对拟合一个模型，我们称之为（l_p）
1.在（s）的定义域上求（u_p），得到上包络的插值（我们称之为（q_u））
1.在（s）的定义域上求（l_p），得到下包络的插值（我们称之为（q_l））。
如您所见，这是三个步骤（查找位置、拟合模型、评估模型）的序列，但应用了两次，一次用于封套的上部，一次用于封套的下部。
要收集（s）的"峰值"，您需要定位（s）的斜率从正变为负的点;要收集（s）的"谷值"，您需要定位（s）的斜率从负变为正的点。
峰值示例：s =[4，5，4] 5 - 4为正4 - 5为负
Flume示例：s =[5，4，5] 4 - 5为负值5 - 4为正值
下面是一个示例脚本，让您开始使用大量的内联注解：

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

这将生成以下输出：

进一步改进的要点：
1.上述代码不 * 过滤 * 峰或谷，这些峰或谷可能出现在比某个阈值"距离"（Tl）（例如时间）更近的地方。这类似于envelope的第二个参数。不过，通过检查u_x,u_y的连续值之间的差异，很容易添加它。
1.但是，对前面提到的一个快速改进是使用移动平均滤波器对数据进行低通滤波在插值上下包络函数之前。（s）与一个合适的移动平均滤波器。这里不去太多的细节（如果需要，可以这样做），要产生一个在N个连续样本上运行的移动平均滤波器，您可以这样做：s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N)。（N）值越大，数据显示越平滑。但请注意，由于平滑滤波器的group delay，这会将值向右移动（N/2）个样本（s_filtered）。有关移动平均值的详细信息，请参阅this link。
希望这个有用。
（如果提供了有关原始应用程序的更多信息，我很乐意修改回复。也许可以用更合适的方式对数据进行预处理（？））

基于@A_A的答案，用nim/max测试替换符号检查，使其更健壮。

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as pt
%matplotlib inline

t = np.multiply(list(range(1000)), .1)
s = 10*np.sin(t)*t**.5

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 range(2,len(s)-1):
    if s[k] >= max(s[:k-1]):
        u_x.append(t[k])
        u_y.append(s[k])

for k in range(2,len(s)-1):
    if s[k] <= min(s[:k-1]):
        l_x.append(t[k])
        l_y.append(s[k])

u_p = scipy.interpolate.interp1d(u_x, u_y, kind = 'cubic', bounds_error = False, fill_value=0.0)
l_p = scipy.interpolate.interp1d(l_x, l_y, kind = 'cubic', bounds_error = False, fill_value=0.0)

q_u = np.zeros(s.shape)
q_l = np.zeros(s.shape)
for k in range(0,len(s)):
    q_u[k] = u_p(t[k])
    q_l[k] = l_p(t[k])

pt.plot(t,s)
pt.plot(t, q_u, 'r')
pt.plot(t, q_l, 'g')

如果您希望函数递增，请尝试：

for k in range(1,len(s)-2):
    if s[k] <= min(s[k+1:]):
        l_x.append(t[k])
        l_y.append(s[k])

用于下包络。

或者你用Pandas。这里我只需要两行代码：

import pandas as pd
import numpy as np

x=np.linspace(0,5*np.pi,1000)
y=np.sin(x)+0.4*np.cos(x/4)*np.sin(x*20)

df=pd.DataFrame(data={"y":y},index=x)

windowsize = 20
df["y_upperEnv"]=df["y"].rolling(window=windowsize).max().shift(int(-windowsize/2))
df["y_lowerEnv"]=df["y"].rolling(window=windowsize).min().shift(int(-windowsize/2))

df.plot(figsize=(20,10))

• • 输出：**

您可能想看看希尔伯特变换，它可能是MATLAB中包络函数背后的实际代码。scipy的信号子模块内置希尔伯特变换，文档中有一个很好的示例，其中提取了振荡信号的包络：https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html

我发现使用scipy函数的组合比其他替代方法性能更好

def envelope(sig, distance):
    # split signal into negative and positive parts
    u_x = np.where(sig > 0)[0]
    l_x = np.where(sig < 0)[0]
    u_y = sig.copy()
    u_y[l_x] = 0
    l_y = -sig.copy()
    l_y[u_x] = 0
    
    # find upper and lower peaks
    u_peaks, _ = scipy.signal.find_peaks(u_y, distance=distance)
    l_peaks, _ = scipy.signal.find_peaks(l_y, distance=distance)
    
    # use peaks and peak values to make envelope
    u_x = u_peaks
    u_y = sig[u_peaks]
    l_x = l_peaks
    l_y = sig[l_peaks]
    
    # add start and end of signal to allow proper indexing
    end = len(sig)
    u_x = np.concatenate((u_x, [0, end]))
    u_y = np.concatenate((u_y, [0, 0]))
    l_x = np.concatenate((l_x, [0, end]))
    l_y = np.concatenate((l_y, [0, 0]))
    
    # create envelope functions
    u = scipy.interpolate.interp1d(u_x, u_y)
    l = scipy.interpolate.interp1d(l_x, l_y)
    return u, l

def test():
    x = np.arange(200)
    sig = np.sin(x)
    u, l = envelope(sig, 1)
    
    plt.figure(figsize=(25,5))
    plt.plot(x, u(x))
    plt.plot(x, l(x))
    plt.plot(x, sig*0.9)
    plt.show()
    
test()