Поиск локальных максимумов / минимумов с помощью Numpy в 1D массиве numpy

В 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 для нахождения максимумов и минимумов соответственно.

для кривых с не слишком много шума, я рекомендую следующий небольшой фрагмент кода:

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

The +1 это важно, потому что diff уменьшает исходный номер индекса.

другой подход (больше слов, меньше кода), который может помочь:

Расположение локальных максимумов и минимумов также в местах пересечения нуля первой производной. Как правило, намного легче найти нулевые пересечения, чем непосредственно найти локальные максимумы и минимумы.

к сожалению, первая производная имеет тенденцию "усиливать" шум, поэтому, когда значительный шум присутствует в исходных данных, первая производная лучше всего используется только после оригинала данные имели некоторую степень сглаживания применяется.

к счастью, довольно часто подходящее ядро может быть создано с помощью простого SWAG ("образованное предположение"). Ширина ядра сглаживания должна быть немного шире, чем самый широкий ожидаемый "интересный" пик в исходных данных, и его форма будет напоминать этот пик (одномерный вейвлет). Для средне-сохраняющих ядер (каким должен быть любой хороший сглаживающий фильтр) сумма элементов ядра должна быть точно равна 1,00, и ядро должно быть симметричным относительно своего центра (это означает, что оно будет иметь нечетное число элементов.

при наличии оптимального ядра сглаживания (или небольшого числа ядер, оптимизированных для различного содержимого данных) степень сглаживания становится масштабирующим фактором для ("усиления") ядра свертки.

определение "правильной" (оптимальной) степени сглаживания (усиления ядра свертки) может быть даже автоматизированный: сравните стандартное отклонение данных первой производной со стандартным отклонением сглаженных данных. Как соотношение двух стандартных отклонений изменяется с изменением степени сглаживания кулачка, можно использовать для прогнозирования эффективных значений сглаживания. Несколько ручных запусков данных (которые действительно представительны) должны быть всем, что нужно.

все предыдущие решения, опубликованные выше, вычисляют первую производную, но они не рассматривают ее как статистическую меру, а также вышеуказанные решения пытаются выполнить функцию сохранения / усиления сглаживания (чтобы помочь тонким пикам "прыгнуть выше" шума).

наконец, плохая новость: поиск "реальных" пиков становится королевской болью, когда шум также имеет функции, которые выглядят как настоящие пики (перекрывающиеся полосы пропускания). Следующим более сложным решением, как правило, является использование более длинного ядра свертки ("более широкая апертура ядра"), которая учитывает взаимосвязь между соседними "реальными" пиками (такими как минимум или максимум скорости возникновения пиков), или использовать несколько проходов свертки с использованием ядер, имеющих разную ширину (но только если это быстрее: это фундаментальная математическая истина, что линейные свертки, выполняемые последовательно, всегда могут быть свернуты вместе в одну свертку). Но часто гораздо легче сначала найти последовательность полезных ядер (различной ширины) и свернуть их вместе, чем непосредственно найти окончательное ядро за один шаг.

надеюсь, это предоставляет достаточно информации, чтобы позволить Google (и, возможно, хороший текст статистики) заполнить пробелы. Я действительно хотел бы иметь время, чтобы предоставить работающий пример или ссылку на него. Если кто-то сталкивается с одним онлайн, пожалуйста, разместите его здесь!

обновление: Я не был доволен градиентом, поэтому я нашел его более надежным в использовании 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()

почему бы не использовать встроенную функцию Scipy сигнал.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))

#generate an inverse numpy 1D arr (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])

чтобы найти локальный max или min, мы по существу хотим найти, когда разница между значениями в списке (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

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 содержат индексы минимумов и максимумов соответственно. Для огромного набора данных он даст много максимумов / минимумов, поэтому в этом случае сначала сгладьте кривую, а затем примените этот алгоритм.