Как получить 2-ю производную сигнала PPG с использованием solve_ivp или аналогичной?

Question

Я хотел бы получить 2-ю производную от необработанного сигнала PPG с использованием solve_ivp, как мне ввести значения сигнала в функцию?

У меня есть набор данных (панды) необработанного сигнала со значениями "y" дляамплитуда для сигнала.Значения «x» или значения времени вычисляются из частоты дискретизации, которая составляет 1000 Гц, что может сделать значения «x» в качестве выборок или преобразовать их в единицы времени.Я подумал, что мне нужно использовать функцию solve_ivp, но может быть и другой путь.Я пробовал np.gradient, но результаты не такие, как ожидалось.Прежде чем пытаться получить 2-ю производную, я использовал полосу пропускания фильтра Баттерворта 0,25-3 Гц, чтобы получить чистый начальный сигнал и снова фильтровать после 1-й производной.

Вот ссылка для необработанного сигнала: https://www.dropbox.com/s/4mhz3tphdzcp4e2/2_1.txt?dl=0

import neurokit as nk
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter, lfilter_zi
import scipy as sc
from scipy import signal


plt.rcParams['figure.figsize'] = [20, 5]

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='bandpass')
    return b, a

def butter_bandpass_filter_zi(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    zi = lfilter_zi(b, a)
    y,zo = lfilter(b, a, data, zi=zi*data[0])
    return y

filename = "0_subject/2_1.txt"

data = pd.read_csv(filename, header=None, sep="\t")
data = data.transpose()
data.plot()
plt.show()

fs = 1000.0
lowcut = 0.15
highcut = 3.0
f0 = 600

T = (len(data)-1)
nsamples = len(data)
# t = np.linspace(0, T, nsamples, endpoint=False)
t = np.linspace(0, T, nsamples, endpoint=False) / fs

y = np.array(data)[:, 0]
plt.figure(2)
plt.clf()
plt.plot(t, y, label='Noisy signal')

x = butter_bandpass_filter_zi(y, lowcut, highcut, fs, order=1)
# print(x)
plt.figure(3)
plt.plot(t, x, label='Filtered signal (%g Hz)' % f0)
plt.xlabel('time (seconds)')
# plt.hlines([-a, a], 0, T, linestyles = '--')
plt.grid(True)
plt.axis('tight')
plt.legend(loc='upper left')

plt.show()

first_derivative = np.gradient(x)
first_derivative_filter = butter_bandpass_filter_zi(first_derivative, lowcut, highcut, fs, order=1)
plt.figure(4)
plt.plot(t, first_derivative_filter)
plt.show()

second_derivative = np.gradient(first_derivative_filter)
second_derivative_filter = butter_bandpass_filter_zi(first_derivative_filter, lowcut, highcut, fs, order=1)
plt.figure(5)
plt.plot(t, second_derivative_filter)
plt.show()

И это то, что я пробовал с solve_ivp, который также дает неправильный результат:

def second_der_func(t, y): return y
sol = solve_ivp(second_der_func, [t[0], len(t)-1], [x[0]])
print(sol.y[0])

plt.plot(sol.t, sol.y[0])
plt.show()

Я должен получить кривую 2-й производной, как это показано на изображении: