NeuroKit icon indicating copy to clipboard operation
NeuroKit copied to clipboard

Different HRV features using Neurokit2 and Heartpy

Open LaxmanSinghTomar opened this issue 3 years ago • 7 comments

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!

LaxmanSinghTomar avatar Apr 28 '21 13:04 LaxmanSinghTomar

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

zen-juen avatar Apr 30 '21 05:04 zen-juen

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!

LaxmanSinghTomar avatar Apr 30 '21 06:04 LaxmanSinghTomar

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.

Tam-Pham avatar May 04 '21 04:05 Tam-Pham

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

LaxmanSinghTomar avatar May 05 '21 04:05 LaxmanSinghTomar

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

Tam-Pham avatar May 18 '21 01:05 Tam-Pham

bump @LaxmanSinghTomar

DominiqueMakowski avatar Oct 19 '21 00:10 DominiqueMakowski

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] avatar Apr 18 '22 18:04 stale[bot]

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.

stale[bot] avatar Sep 08 '22 17:09 stale[bot]