如何获取信号的高和低包络? [英] python - How to get high and low envelope of a signal?

问题描述

我的数据非常嘈杂，我正在尝试确定信号的上限和下限.有点像MATLAB中的这个示例:

I have quite a noisy data, and I am trying to work out a high and low envelope to the signal. It is kinda of like this example in MATLAB:

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

在提取峰包络"中. Python中是否有类似的功能可以做到这一点?我的整个项目都是用Python编写的，在最坏的情况下，我可以提取numpy数组并将其放入MATLAB中并使用该示例.但是我更喜欢matplotlib的外观...而且真的是cba在MATLAB和Python之间执行所有这些I/O ...

in "Extracting Peak Envelope". Is there a similar function in Python that can do that? My entire project has been written in Python, worst case scenario I can extract my numpy array and throw it into MATLAB and use that example. But I prefer the look of matplotlib... and really cba doing all of those I/O between MATLAB and Python...

谢谢

推荐答案

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

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

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

As far as I am aware there is no such function in Numpy / Scipy / Python. However, it is not that difficult to create one. The general idea is as follows:

给出值的向量:

1. 找到(s)的峰的位置.我们称它们为(u)
2. 找到s的槽的位置.我们称它们为(l).
3. 使模型适合(u)个值对.我们称它为(u_p)
4. 使模型适合(l)个值对.我们称它为(l_p)
5. 在(s)的范围内求值(u_p)，以获得上包络线的插值. (我们叫他们(q_u))
6. 在(s)的范围内评估(l_p)，以获得较低包络线的插值. (我们叫他们(q_l)).

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

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

要收集(s)的峰"，您需要找到(s)的斜率从正变为负的点，并且要收集(s)的波谷"，需要找到斜率的在点(s)从负变为正.

To collect the "peaks" of (s) you need to locate points where the slope of (s) changes from positive to negative and to collect the "troughs" of (s) you need to locate the points where the slope of (s) changes from negative to positive.

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

A peak example: s = [4,5,4] 5-4 is positive 4-5 is negative

一个低谷例子:s = [5,4,5] 4-5为负5-4为正

A trough example: s = [5,4,5] 4-5 is negative 5-4 is positive

以下是一个示例脚本，可让您开始获得大量内联注释:

Here is an example script to get you started with plenty 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()

这将产生以下输出:

进一步改进的要点:

1. 上面的代码没有过滤可能比某个阈值距离"(Tl)(例如时间)更近的峰或谷.这类似于envelope的第二个参数.通过检查u_x,u_y的连续值之间的差异，可以很容易地添加它.

1. The above code does not filter peaks or troughs that may be occuring closer than some threshold "distance" (Tl) (e.g. time). This is similar to the second parameter of envelope. It is easy to add it though by examining the differences between consecutive values of u_x,u_y.

但是，对上述要点的快速改进是使用移动平均值滤波器之前插值上下包络函数对数据进行低通滤波.您可以通过使用适当的移动平均滤波器对卷积进行卷积来轻松实现此目的.无需在此处进行详细介绍(如果需要，可以做)，以生成可对N个连续样本进行操作的移动平均滤波器，您可以执行以下操作:s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N). (N)越高，数据将显示得越平滑.但是请注意，由于称为此链接.

However, a quick improvement over the point mentioned previously is to lowpass filter your data with a moving average filter BEFORE interpolating an upper and lower envelope functions. You can do this easily by convolving your (s) with a suitable moving average filter. Without going to a great detail here (can do if required), to produce a moving average filter that operates over N consecutive samples, you would do something like this: s_filtered = numpy.convolve(s, numpy.ones((1,N))/float(N). The higher the (N) the smoother your data will appear. Please note however that this will shift your (s) values (N/2) samples to the right (in s_filtered) due to something that is called group delay of the smoothing filter. For more information about the moving average, please see this link.

希望这会有所帮助.

(如果提供了有关原始应用程序的更多信息，请高兴地回应一下.也许可以用更合适的方式对数据进行预处理(?)

(Happy to ammend the response if more information about the original application is provided. Perhaps the data can be pre-processed in a more suitable way (?) )

这篇关于如何获取信号的高和低包络?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！