使用 Numpy 在一维 Numpy 数组中查找局部极大值/极小值

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

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

也可以在此步骤之前使用 numpy.convolve()创建数组 光滑。

我不认为有专门的功能。

更新: 我不满意梯度，所以我发现它更可靠的使用 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()

对于没有太多噪音的曲线，我推荐以下小代码片段:

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减少了原始索引数。

In 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 的局部 max/min 索引:

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

scipy.signal还分别提供 argrelmax和 argrelmin来寻找极值。

另一种方法(增加词汇量，减少代码量)可能会有所帮助:

本地极值的位置也是一阶导数过零点的位置。与直接找到当地极值相比，找到零交叉点通常要容易得多。

不幸的是，一阶导数往往会“放大”噪声，因此当原始数据中存在显著噪声时，最好在原始数据经过一定程度的平滑处理后才使用一阶导数。

由于平滑，在最简单的意义上，是一个低通滤波器，平滑往往是最好(嗯，最容易)使用卷积内核，并“塑造”内核可以提供惊人数量的特征保持/增强能力。寻找最佳内核的过程可以通过多种方法自动化，但最好的方法可能是简单的蛮力(足够快地寻找小内核)。一个好的内核将(如预期的那样)大规模地扭曲原始数据，但是它不会影响感兴趣的峰/谷的位置。

幸运的是，通常可以通过一个简单的 SWAG (“有根据的猜测”)来创建合适的内核。平滑核的宽度应该比原始数据中预期的最宽的“有趣”峰稍宽一些，它的形状将类似于那个峰(一个单一尺度的小波)。对于均值保持的内核(任何好的平滑滤波器都应该是这样的) ，内核元素的和应该精确地等于1.00，内核应该围绕它的中心对称(这意味着它将有奇数个元素。

给定一个最优的平滑内核(或少量针对不同数据内容优化的内核) ，平滑程度成为卷积内核的缩放因子(“增益”)。

确定“正确”(最佳)的平滑程度(卷积内核增益)甚至可以自动化: 比较一阶导数数据的标准差和平滑数据的标准差。如何利用两个标准差的比值随光顺程度的变化来预测有效的光顺值。一些手动数据运行(真正具有代表性的)应该就足够了。

上面提到的所有先前的解决方案都计算了一阶导数，但是他们并没有把它当作一个统计度量，上面的解决方案也没有尝试执行特征保持/增强平滑(以帮助微妙的峰值“跳过”噪音)。

最后，坏消息是: 当噪音也具有看起来像真正的峰值(重叠带宽)的特征时，找到“真正的”峰值就成了皇家的痛苦。下一个更复杂的解决方案通常是使用一个更长的卷积内核(一个“更宽的内核孔径”) ，考虑到相邻的“实”峰之间的关系(例如峰出现的最小或最大速率) ，或者使用具有不同宽度的内核的多次卷积传递(但只有在它更快的情况下: 这是一个基本的数学真理，按顺序进行的线性卷积总是可以卷积在一起成为一个单一的卷积)。但是，首先找到一系列有用的内核(宽度不同)并将它们卷积在一起，往往比在一个步骤中直接找到最终的内核要容易得多。

希望这提供了足够的信息，让谷歌(也许是一个很好的统计文本)填补空白。我真的希望我有时间提供一个工作的例子，或一个链接到一个。如果有人在网上看到，请把它贴在这里！

import numpy as np
x=np.array([6,3,5,2,1,4,9,7,8])
y=np.array([2,1,3,5,3,9,8,10,7])
sortId=np.argsort(x)
x=x[sortId]
y=y[sortId]
minm = np.array([])
maxm = np.array([])
i = 0
while i < length-1:
if i < length - 1:
while i < length-1 and y[i+1] >= y[i]:
i+=1


if i != 0 and i < length-1:
maxm = np.append(maxm,i)


i+=1


if i < length - 1:
while i < length-1 and y[i+1] <= y[i]:
i+=1


if i < length-1:
minm = np.append(minm,i)
i+=1




print minm
print maxm

minm和 maxm分别包含极值指数。对于一个庞大的数据集，它会给出大量的最大值/最小值，因此在这种情况下，首先平滑曲线，然后应用该算法。

为什么不使用 Scipy 内置函数 Find _ Peak _ 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]

问候

这些解决方案都不适合我，因为我也想在重复值的中心找到峰值。例如，在

ar = np.array([0,1,2,2,2,1,3,3,3,2,5,0])

答案应该是

array([ 3,  7, 10], dtype=int64)

我用了一个循环，我知道它不是很干净，但它完成了任务。

def findLocalMaxima(ar):
# find local maxima of array, including centers of repeating elements
maxInd = np.zeros_like(ar)
peakVar = -np.inf
i = -1
while i < len(ar)-1:
#for i in range(len(ar)):
i += 1
if peakVar < ar[i]:
peakVar = ar[i]
for j in range(i,len(ar)):
if peakVar < ar[j]:
break
elif peakVar == ar[j]:
continue
elif peakVar > ar[j]:
peakInd = i + np.floor(abs(i-j)/2)
maxInd[peakInd.astype(int)] = 1
i = j
break
peakVar = ar[i]
maxInd = np.where(maxInd)[0]
return maxInd

我相信 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

从 SciPy 版本1.1开始，您也可以使用 找到山峰。

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

还有一个:

def local_maxima_mask(vec):
"""
Get a mask of all points in vec which are local maxima
:param vec: A real-valued vector
:return: A boolean mask of the same size where True elements correspond to maxima.
"""
mask = np.zeros(vec.shape, dtype=np.bool)
greater_than_the_last = np.diff(vec)>0  # N-1
mask[1:] = greater_than_the_last
mask[:-1] &= ~greater_than_the_last
return mask

另一种基本上使用膨胀算子的解决方案:

import numpy as np
from scipy.ndimage import rank_filter


def find_local_maxima(x):
x_dilate = rank_filter(x, -1, size=3)
return x_dilate == x

至于最低限度:

def find_local_minima(x):
x_erode = rank_filter(x, -0, size=3)
return x_erode == x

另外，从 scipy.ndimage你可以用 grey_dilation代替 rank_filter(x, -1, size=3)，用 grey_erosion代替 rank_filter(x, 0, size=3)。这不需要本地排序，因此速度稍快。

还有一个答案。

这个程序不需要额外的包(numpy 除外),

points = [ 0, 0, 1, 2, 3, 3, 2, 2, 3, 1, 1 ]
minimums   ^  ^              ^  ^     ^  ^

将返回所有本地最小值的列表

result = [ 0, 1, 6, 7, 9, 10 ]

它可以很容易地扩展到也寻找最大值。

def find_valleys(points: np.ndarray, edges=True) -> list:
"""
Find the indices of all points that are local minimums.


:param np.ndarray points: a 1D array of numeric data
:param bool edges: allows the first and last indices to be returned, defaults to True
:return list: a list of integers, indices into the array
"""
dif = np.diff(points)
p = -1 if edges else 1
s = 0
result = []
for i,d in enumerate(dif):
if d < 0: s = i + 1
if p < 0 and d > 0:   # found a valley
result.extend(range(s,i + 1))
if d: p = d
if p < 0 and edges:
result.extend(range(s,i + 2))
return result