neuropsychology / NeuroKit

NeuroKit2: The Python Toolbox for Neurophysiological Signal Processing
https://neuropsychology.github.io/NeuroKit
MIT License
1.6k stars 426 forks source link

Different HRV features using Neurokit2 and Heartpy #448

Closed LaxmanSinghTomar closed 2 years ago

LaxmanSinghTomar commented 3 years ago

I'm trying to obtain HRV features from both Neurokit2 and HeartPy. Although, both are making use of same peaks list, the results are different.

Neurokit2 Screenshot_2021-04-28 Google Colaboratory(3)

HeartPy Screenshot_2021-04-28 Google Colaboratory(4)

I tried looking up the source code and the implementation is exactly the same too. Am I missing something here? Any explanation for this behavior or possible solution will be really helpful!

zen-juen commented 3 years ago

Hi @LaxmanSinghTomar, could you include your code here so that we can reproduce the results and investigate? 😃

LaxmanSinghTomar commented 3 years ago

Sure. Here they are:

Neurokit2:

import neurokit2 as nk
import heartpy as hp
from heartpy.datautils import rolling_mean, _sliding_window

def peakslist(el):
    data = np.array(el)
    new_data = hp.enhance_peaks(data, iterations=2)
    filtered = hp.filtering.filter_signal(new_data, cutoff = 0.75, sample_rate = 100.0, order = 3, filtertype='highpass')
    rol_mean = rolling_mean(filtered, windowsize = 0.75, sample_rate = 100.0)
    wd = hp.peakdetection.fit_peaks(filtered, rol_mean, sample_rate = 100.0)
    return wd['peaklist']

info = peakslist(list(chain(*total.iloc[0, 1:2].values)))
print(info)
hrv_indices = nk.hrv(info, sampling_rate = 100, show=False)
hrv_indices

It generates the first output in my previous comment. For the sake of reference, this is the info or list of peaks:

[  28  197  270  342  412  481  548  615  681  747  805  865  928  993
 1059 1121 1185 1223 1252 1253 1325 1348 1393 1413 1474 1493 1494 1547
 1588 1638 1669 1704 1770]

HeartPy:

data = np.array(list(chain(*total.iloc[0, 1:2].values)))
new_data = hp.enhance_peaks(data, iterations=2)
filtered = hp.filtering.filter_signal(new_data, cutoff = 0.75, sample_rate = 100.0, order = 3, filtertype='highpass')
wd, m = hp.process(filtered, sample_rate = 100.0)
print(wd['peaklist'])
wd, m = hp.analysis.calc_fd_measures(method = 'fft', measures = m, working_data = wd)
m

It generates the second output in my previous comment. And again, this is the list of peaks:

[  28  197  270  342  412  481  548  615  681  747  805  865  928  993
 1059 1121 1185 1223 1252 1253 1325 1348 1393 1413 1474 1493 1494 1547
 1588 1638 1669 1704 1770]

Hoping this clarifies the issue a little more!

Tam-Pham commented 3 years ago

hey @LaxmanSinghTomar

I would like to clarify a little bit about the code you shared. Could you share the exact input for us (i.e. list(chain(*total.iloc[0, 1:2].values))) so that we can fully replicate your output?

I tried looking up the source code and the implementation is exactly the same too.

if the implementation is the same, it's usually the case that there might be an extra step done somewhere in the analysis pipeline or that the input into the two functions is not the same. I'm not very familiar with the implementation pipeline in heartpy but to have such drastic differences in indices, especially in the time-domain, seems strange.

LaxmanSinghTomar commented 3 years ago

Sure, that line specifies that from a pandas dataframe read the first row and second column which is PPG signal. List and Chain have been used to only get the PPG signal values in a list. I'm attaching the CSV file which contains PPG values in IRLED column. total is nothing but the variable name assigned to this dataframe.

606863ad786baa05288d6d29.csv

Tam-Pham commented 3 years ago

hey @LaxmanSinghTomar

Would appreciate it if you can give us the exact reproducible piece of code so that replicate your exact output and sort out what's wrong. Specifically, how's your chain function defined?

import numpy as np
import pandas as pd

import neurokit2 as nk
import heartpy as hp
from heartpy.datautils import rolling_mean, _sliding_window

def peakslist(el):
    data = np.array(el)
    new_data = hp.enhance_peaks(data, iterations=2)
    filtered = hp.filtering.filter_signal(new_data, cutoff = 0.75, sample_rate = 100.0, order = 3, filtertype='highpass')
    rol_mean = rolling_mean(filtered, windowsize = 0.75, sample_rate = 100.0)
    wd = hp.peakdetection.fit_peaks(filtered, rol_mean, sample_rate = 100.0)
    return wd['peaklist']

total = pd.read_csv('total.csv')

peakslist(list(chain(*total.iloc[0, 1:2].values)))
print(info)
hrv_indices = nk.hrv(info, sampling_rate = 100, show=False)
hrv_indices

data = np.array(list(chain(*total.iloc[0, 1:2].values)))
new_data = hp.enhance_peaks(data, iterations=2)
filtered = hp.filtering.filter_signal(new_data, cutoff = 0.75, sample_rate = 100.0, order = 3, filtertype='highpass')
wd, m = hp.process(filtered, sample_rate = 100.0)
print(wd['peaklist'])
wd, m = hp.analysis.calc_fd_measures(method = 'fft', measures = m, working_data = wd)
m
DominiqueMakowski commented 3 years ago

bump @LaxmanSinghTomar

stale[bot] commented 2 years ago

This issue has been automatically marked as inactive because it has not had recent activity. It will eventually be closed if no further activity occurs.

stale[bot] commented 2 years ago

Is this still relevant? If so, what is blocking it? Is there anything you can do to help move it forward?

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs.