我有相当嘈杂的数据，我正在尝试计算信号的高低包络。这有点像 MATLAB 中的这个例子：
http://uk.mathworks.com/help/signal/examples/signal-smoothing.html
在“提取峰包络”中。 Python中是否有类似的功能可以做到这一点？我的整个项目都是用 Python 编写的，最坏的情况是我可以提取我的 numpy 数组并将其放入 MATLAB 并使用该示例。但我更喜欢 matplotlib 的外观......而且真的是 cba 在 MATLAB 和 Python 之间完成所有这些 I/O......
谢谢，
【问题讨论】：
标签： python matlab numpy matplotlib signal-processing
Python中有没有类似的函数可以做到这一点？
据我所知，Numpy / Scipy / Python 中没有这样的功能。但是，创建一个并不难。大致思路如下：
给定一个值向量：
如您所见，它是三个步骤（查找位置、拟合模型、评估模型）的顺序，但应用了两次，一次用于信封的上部，一次用于下部。
要收集 (s) 的“峰值”，您需要定位 (s) 的斜率从正变为负的点，并收集 (s) 的“谷”，您需要定位斜率所在的点(s) 由负变为正。
峰值示例：s = [4,5,4] 5-4 为正 4-5 为负
一个低谷示例：s = [5,4,5] 4-5 是负数 5-4 是正数
这是一个示例脚本，可帮助您开始使用大量内联 cmets：
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()
这会产生这个输出：
进一步改进的要点：
上述代码不会过滤可能出现在比某个阈值“距离”（Tl）（例如时间）更近的峰或谷。这类似于envelope的第二个参数。通过检查u_x,u_y 的连续值之间的差异，很容易添加它。
但是，对前面提到的一点的快速改进是使用移动平均滤波器在插入上包络函数和下包络函数来对数据进行低通滤波。您可以通过将您的 (s) 与合适的移动平均滤波器进行卷积来轻松做到这一点。无需在这里详细介绍（如果需要，可以这样做），要生成对 N 个连续样本进行操作的移动平均滤波器，您可以执行以下操作：s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N)。 (N) 越高，您的数据就越平滑。但是请注意，由于平滑滤波器的称为group delay 的东西，这会将您的 (s) 值 (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])
用于较低的信封。
【讨论】：
第一次尝试是使用 scipy Hilbert transform 来确定幅度包络，但这在许多情况下并没有按预期工作，主要是因为引用 digital signal processing answer：
希尔伯特包络，也称为能量时间曲线 (ETC)，仅适用于 用于窄带波动。产生一个分析信号，其中 你后来取绝对值，是一个线性运算，所以它对待 信号的所有频率均等。如果你给它一个纯正弦 波，它确实会以一条直线返回给你。当你给它 但是，您可能会得到白噪声。
从那时起，由于其他答案使用三次样条插值并且确实会变得很麻烦，有点不稳定（虚假振荡）并且对于非常长且嘈杂的数据阵列非常耗时，我将在这里提供一个简单且高效的 numpy似乎工作得很好的版本：
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)
high_idx, low_idx = hl_envelopes_idx(s)
# plot
plt.plot(t,s,label='signal')
plt.plot(t[high_idx], s[high_idx], 'r', label='low')
plt.plot(t[low_idx], s[low_idx], '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))
high_idx, low_idx = hl_envelopes_idx(s,dmin=15,dmax=15)
# plot
plt.plot(t,s,label='signal')
plt.plot(t[high_idx], s[high_idx], 'r', label='low')
plt.plot(t[low_idx], s[low_idx], 'g', label='high')
示例 3：非对称调制啁啾
18867925 样本的复杂得多的信号（此处不包括在内）：
【讨论】：
您可能想查看希尔伯特变换，这可能是 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()
【讨论】：
或者你使用熊猫。这里我只需要两行代码：
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))
输出：
【讨论】：