在1D numpy数组中使用Numpy查找局部最大值/最小值

Question

你能否建议numpy/scipy的模块函数可以在1D numpy数组中找到局部最大值/最小值？显然,最简单的方法是看看最近的邻居,但我希望有一个公认的解决方案,这是numpy发行版的一部分.

Answer 1

在SciPy> = 0.11

import numpy as np
from scipy.signal import argrelextrema

x = np.random.random(12)

# for local maxima
argrelextrema(x, np.greater)

# for local minima
argrelextrema(x, np.less)

产生

>>> x
array([ 0.56660112,  0.76309473,  0.69597908,  0.38260156,  0.24346445,
    0.56021785,  0.24109326,  0.41884061,  0.35461957,  0.54398472,
    0.59572658,  0.92377974])
>>> argrelextrema(x, np.greater)
(array([1, 5, 7]),)
>>> argrelextrema(x, np.less)
(array([4, 6, 8]),)

注意,这些是x的索引,是本地最大/最小值.要获取值,请尝试:

>>> x[argrelextrema(x, np.greater)[0]]

scipy.signal还提供argrelmax和argrelmin分别查找最大值和最小值.

Answer 2

如果您要查找1d数组a中小于其邻居的所有条目,您可以尝试

numpy.r_[True, a[1:] < a[:-1]] & numpy.r_[a[:-1] < a[1:], True]

您还可以在使用此步骤之前平滑阵列numpy.convolve().

我认为没有专门的功能.

Answer 3

对于没有太多噪音的曲线,我建议使用以下小代码片段:

from numpy import *

# example data with some peaks:
x = linspace(0,4,1e3)
data = .2*sin(10*x)+ exp(-abs(2-x)**2)

# that's the line, you need:
a = diff(sign(diff(data))).nonzero()[0] + 1 # local min+max
b = (diff(sign(diff(data))) > 0).nonzero()[0] + 1 # local min
c = (diff(sign(diff(data))) < 0).nonzero()[0] + 1 # local max


# graphical output...
from pylab import *
plot(x,data)
plot(x[b], data[b], "o", label="min")
plot(x[c], data[c], "o", label="max")
legend()
show()

这+1很重要,因为diff减少了原始索引号.

Answer 4

另一种方法(更多的话,更少的代码)可能会有所帮助:

局部最大值和最小值的位置也是一阶导数的过零点的位置.找到零交叉通常比直接找到局部最大值和最小值要容易得多.

不幸的是,一阶导数倾向于"放大"噪声,因此当原始数据中存在显着噪声时,只有在原始数据已经应用了一定程度的平滑之后才能最好地使用一阶导数.

由于平滑在最简单的意义上是低通滤波器,因此通过使用卷积内核进行平滑通常是最好的(很好,最容易),并且"整形"内核可以提供惊人数量的特征保留/增强功能.找到最佳内核的过程可以使用各种方法自动完成,但最好的方法可能是简单的强力(很快找到小内核).一个好的内核将(按预期)大量扭曲原始数据,但它不会影响感兴趣的峰/谷的位置.

幸运的是,通常可以通过简单的SWAG("有根据的猜测")创建合适的内核.平滑内核的宽度应该比原始数据中最宽的预期"有趣"峰值宽一些,并且其形状将类似于该峰值(单级小波).对于保持平均值的内核(任何好的平滑滤波器应该是什么),内核元素的总和应该精确地等于1.00,并且内核应该关于其中心对称(意味着它将具有奇数个元素.

给定最佳平滑内核(或针对不同数据内容优化的少量内核),平滑程度成为卷积内核("增益")的缩放因子.

确定"正确"(最佳)平滑程度(卷积核增益)甚至可以自动化:将一阶导数数据的标准偏差与平滑数据的标准偏差进行比较.两个标准偏差的比率如何随平滑度的变化而变化,可用于预测有效平滑值.一些手动数据运行(真正具有代表性)应该是所需要的.

上面公布的所有先前解决方案都计算出一阶导数,但它们并未将其视为统计测量,上述解决方案也不会尝试执行特征保留/增强平滑(以帮助微妙的峰值"跳过"噪声).

最后,坏消息:当噪声也具有看起来像真正的峰值(重叠带宽)的特征时,找到"真正的"峰值变成了皇家的痛苦.下一个更复杂的解决方案通常是使用更长的卷积核("更宽的核孔径"),其考虑相邻"实际"峰值之间的关系(例如峰值出现的最小或最大速率),或使用多个卷积通过使用具有不同宽度的内核(但只有当它更快时:基本的数学事实是,顺序执行的线性卷积总是可以卷积成单个卷积).但是,通常更容易找到一系列有用的内核(宽度不同)并将它们卷积在一起,而不是直接在一个步骤中找到最终的内核.

希望这提供足够的信息让谷歌(也许是一个好的统计文本)填补空白.我真的希望我有时间提供一个有用的例子,或一个链接到一个.如果有人在网上遇到一个,请在这里发布!

Answer 5

从SciPy 1.1版开始，您还可以使用find_peaks。以下是从文档本身获取的两个示例。

使用该height参数，可以选择高于某个阈值的所有最大值（在此示例中，是所有非负最大值；这在必须处理嘈杂的基线时非常有用；如果要查找最小值，只需将输入乘以通过-1）：

import matplotlib.pyplot as plt
from scipy.misc import electrocardiogram
from scipy.signal import find_peaks
import numpy as np

x = electrocardiogram()[2000:4000]
peaks, _ = find_peaks(x, height=0)
plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.plot(np.zeros_like(x), "--", color="gray")
plt.show()

另一个非常有用的参数是distance，它定义了两个峰之间的最小距离：

peaks, _ = find_peaks(x, distance=150)
# difference between peaks is >= 150
print(np.diff(peaks))
# prints [186 180 177 171 177 169 167 164 158 162 172]

plt.plot(x)
plt.plot(peaks, x[peaks], "x")
plt.show()

Answer 6

为什么不使用Scipy内置函数signal.find_peaks_cwt来完成这项工作？

from scipy import signal
import numpy as np

#generate junk data (numpy 1D arr)
xs = np.arange(0, np.pi, 0.05)
data = np.sin(xs)

# maxima : use builtin function to find (max) peaks
max_peakind = signal.find_peaks_cwt(data, np.arange(1,10))

# inverse  (in order to find minima)
inv_data = 1/data
# minima : use builtin function fo find (min) peaks (use inversed data)
min_peakind = signal.find_peaks_cwt(inv_data, np.arange(1,10))

#show results
print "maxima",  data[max_peakind]
print "minima",  data[min_peakind]

结果:

maxima [ 0.9995736]
minima [ 0.09146464]

问候

Answer 7

我相信在 numpy（单行）中有一种更简单的方法。

import numpy as np

list = [1,3,9,5,2,5,6,9,7]

np.diff(np.sign(np.diff(list))) #the one liner

#output
array([ 0, -2,  0,  2,  0,  0, -2])

为了找到局部最大值或最小值，我们本质上想要找到列表中的值之间的差异 (3-1, 9-3...) 何时从正变为负 (max) 或从负变为正 (min)。因此，我们首先找出不同之处。然后我们找到符号，然后我们通过再次取差来找到符号的变化。（有点像微积分中的一阶和二阶导数，只有我们有离散数据并且没有连续函数。）

我的示例中的输出不包含极值（列表中的第一个和最后一个值）。此外，就像微积分一样，如果二阶导数为负，则为最大值，如果为正，则为最小值。

因此，我们有以下匹配：

[1,  3,  9,  5,  2,  5,  6,  9,  7]
    [0, -2,  0,  2,  0,  0, -2]
        Max     Min         Max

Answer 8

更新: 我对渐变感到不满意,所以我发现使用它更可靠numpy.diff.请告诉我它是否符合您的要求.

关于噪声问题,数学问题是如果我们想要查看噪声我们可以使用前面提到的卷积之类的噪声来定位最大值/最小值.

import numpy as np
from matplotlib import pyplot

a=np.array([10.3,2,0.9,4,5,6,7,34,2,5,25,3,-26,-20,-29],dtype=np.float)

gradients=np.diff(a)
print gradients


maxima_num=0
minima_num=0
max_locations=[]
min_locations=[]
count=0
for i in gradients[:-1]:
        count+=1

    if ((cmp(i,0)>0) & (cmp(gradients[count],0)<0) & (i != gradients[count])):
        maxima_num+=1
        max_locations.append(count)     

    if ((cmp(i,0)<0) & (cmp(gradients[count],0)>0) & (i != gradients[count])):
        minima_num+=1
        min_locations.append(count)


turning_points = {'maxima_number':maxima_num,'minima_number':minima_num,'maxima_locations':max_locations,'minima_locations':min_locations}  

print turning_points

pyplot.plot(a)
pyplot.show()