Как найти несколько максимумов числовых временных рядов - PullRequest
Новая
       16

Как найти несколько максимумов числовых временных рядов

0 голосов
/ 19 января 2019

Я численно (на питоне) решаю дифференциальное уравнение первого порядка, используя метод Эйлера.Я строю решение от времени t = 0 до некоторого произвольного времени, продолжающегося с шагом времени, скажем, размером 0,05.Один из примеров полученного решения приведен на рисунке ниже enter image description here

Я хотел бы найти максимумы, которые мы видим на этом рисунке (а также время их появления), и сохранитьих в словаре.Если бы во всем временном диапазоне был только один максимум, я мог бы использовать этот код 

y=[xinitial,viinitial,ainitial]            
t=0
maximum=-20000
maxai={}
h=0.05
ai=-2.1
for i in range(0,3701):
   dydt=computederivs(y)
   y = euler(y,dydt,h)
   t+=h
   if y[0]>maximum:
      maximum = y[0]
 maxai[ai]=maximum

Поскольку у меня есть несколько локальных максимумов, я должен их обнаруживать по мере движения во времени t, каким-то образом проверяя, когда функция завершается.вниз после восхождения на несколько шагов.Мне также нужно хранить максимумы в списке, который является значением ключа словаря.Я полагаю, что это достаточно распространенная задача, что должны быть хорошо известные способы сделать это более или менее просто?

Ответы [ 2 ]

0 голосов
/ 19 января 2019

Я бы предложил использовать find_peaks из scipy.signal модуля. Эта функция принимает одномерный массив и находит все локальные максимумы путем простого сравнения соседних значений. При желании подмножество этих пиков может быть выбрано путем указания условий для свойств пика.

Вот фрагмент кода для начала работы: 

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

Fs = 8000
f = 5
sample = 8000
x = np.arange(sample)
y = np.sin(2 * np.pi * f * x / Fs)
peaks = find_peaks(y)

plt.scatter(peaks[0], np.ones(f), c='red')
plt.plot(x, y)
plt.xlabel('sample(n)')
plt.ylabel('voltage(V)')
plt.show()

enter image description here

Ваши максимумы: 

print(peaks[0])
[ 400 2000 3600 5200 6800]
0 голосов
/ 19 января 2019

Событие, которое вы должны искать, - это когда первая производная v меняет свой знак.Чтобы получить максимум, проверьте, является ли вторая производная a отрицательной. 

 yp, y = y, euler(y,dydt,h)
 if py[1]*y[1] < 0 and y[2] < 0:
      # do what you need to do to list the maximum

Вы можете использовать линейную интерполяцию, чтобы получить более точную ординату для максимальной позиции.

ЕслиВы хотите получить серьезные результаты, используйте метод более высокого порядка для интеграции ODE.

...