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

Можете ли вы предложить функцию модуля из numpy / scipy, которая может находить локальные максимумы / минимумы в массиве 1D numpy? Очевидно, самый простой подход - взглянуть на ближайших соседей, но я хотел бы иметь приемлемое решение, которое является частью дистрибутива numpy.

Если вы ищете все записи в массиве 1d, aменьшие, чем их соседи, вы можете попробовать

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

Вы также можете сгладить свой массив перед этим шагом, используя numpy.convolve().

Я не думаю, что для этого есть специальная функция.

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

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

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

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

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

Поскольку сглаживание, в простейшем смысле, является фильтром нижних частот, сглаживание часто лучше (ну, проще всего) выполняется с использованием ядра свертки, и «формирование» этого ядра может обеспечить удивительное количество возможностей сохранения / улучшения функций. , Процесс поиска оптимального ядра можно автоматизировать с помощью различных средств, но лучшим может быть простой перебор (достаточно быстро для поиска небольших ядер). Хорошее ядро ​​(как и предполагалось) сильно исказит исходные данные, но НЕ повлияет на расположение интересующих пиков / впадин.

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

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

Определение «правильной» (оптимальной) степени сглаживания (усиления ядра свертки) можно даже автоматизировать: сравните стандартное отклонение данных первой производной со стандартным отклонением сглаженных данных. То, как соотношение двух стандартных отклонений изменяется с изменением степени сглаживания кулачка, можно использовать для прогнозирования эффективных значений сглаживания. Все, что нужно, - это несколько ручных прогонов данных (которые действительно репрезентативны).

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

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

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

Начиная с версии 1.1 SciPy, вы также можете использовать 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()

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

С уважением

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

Пока этот вопрос действительно старый. Я считаю, что в 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). Поэтому сначала находим разницу. Затем мы находим знак, а затем мы находим изменения знака, снова взяв разницу. (Это что-то вроде первой и второй производных в исчислении, только у нас есть дискретные данные и нет непрерывной функции.)

Результат в моем примере не содержит экстремумов (первое и последнее значения в списке). Кроме того, как и в исчислении, если вторая производная отрицательна, у вас есть max, а если она положительна, у вас есть min.

Таким образом, у нас получился следующий матч:

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

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

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 

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

Другое решение, использующее по существу оператор расширения:

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вас могут заменить rank_filter(x, -1, size=3)с grey_dilationи rank_filter(x, 0, size=3)с grey_erosion. Это не требует локальной сортировки, поэтому выполняется немного быстрее.

Другой:

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