python - 如何获得信号的高低信封？

Question

我有相当嘈杂的数据,我正试图找出信号的高低信封.它有点像MATLAB中的这个例子:

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

在"提取峰值信封"中.Python中是否有类似的功能可以做到这一点？我的整个项目都是用Python编写的,最糟糕的情况我可以提取我的numpy数组并将其抛入MATLAB并使用该示例.但是我更喜欢matplotlib的外观......而且真正的cba在MATLAB和Python之间完成所有这些I/O ......

谢谢,

Answer 1

Python中是否有类似的功能可以做到这一点？

据我所知,Numpy/Scipy/Python中没有这样的功能.但是,创建一个并不困难.总体思路如下:

给定值的向量:

1. 找到(s)的峰的位置.我们叫他们(你)
2. 找到s的低谷的位置.我们称之为(l).
3. 使模型适合(u)值对.我们称之为(u_p)
4. 使模型适合(l)值对.我们称之为(l_p)
5. 评估(s)域上的(u_p)以获得上包络的内插值.(我们称之为(q_u))
6. 评估(s)域上的(l_p)以获得下包络的内插值.(我们称之为(q_l)).

如您所见,它是三个步骤(查找位置,拟合模型,评估模型)的序列,但应用了两次,一次用于信封的上半部分,一次用于下部.

要收集(s)的"峰值",您需要找到(s)的斜率从正变为负的点,并收集(​​s)的"波谷",您需要找到斜率为(s)的点. )从负面变为正面.

峰值示例:s = [4,5,4] 5-4为正4-5为负

一个低谷的例子: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. 上述代码不过滤可能比某个阈值"距离"(T1)(例如时间)更近的峰值或波谷.这类似于第二个参数envelope.通过检查连续值之间的差异,可以很容易地添加它u_x,u_y.

2. 但是,与前面提到的相比,快速改进的方法是在插入上下包络函数之前,使用移动平均滤波器对数据进行低通滤波.您可以通过将您的(s)与合适的移动平均滤波器进行卷积来轻松完成此操作.如果没有详细说明(如果需要可以做的话),要产生一个在N个连续样本上运行的移动平均滤波器,你可以这样做:s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N).数据(N)越高,您的数据就会越平滑.但请注意,s_filtered由于平滑滤波器的群延迟,会将您的值(N/2)样本移到右侧(in ).有关移动平均线的更多信息,请参阅此链接.

希望这可以帮助.

(如果提供了有关原始申请的更多信息,请尽快提出答复.也许数据可以更合适的方式进行预处理(？))

Answer 2

第一次尝试是利用 scipy Hilbert 变换来确定幅度包络，但这在许多情况下没有按预期工作，主要是因为引用了这个数字信号处理答案：

希尔伯特包络，也称为能量-时间曲线 (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

例一：准周期振动

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样本信号（这里不包括）：