paulvangentcom / heartrate_analysis_python

Python Heart Rate Analysis Package, for both PPG and ECG signals

HeartPy - Python Heart Rate Analysis Toolkit

Like HeartPy? Don't forget to leave a star!

Structural update

HeartPy V1.2 has landed! The structure of the package has been reworked to be in separate modules now in preparation of the next big update, which will feature many analysis expansions and the first steps towards a GUI for HeartPy. HeartPy has been growing steadily and had reached the point where it became cluttered and unwieldy to keep in a single file. The API remains unchanged.

An 'Examples' folder has been added to the repo which will be expanded soon. Now there's two notebooks explaining how to analyse ppg signals from smartwatches and smart rings.

Colorblind support has been added, see this notebook in the examples folder

Installation

python setup.py install

Alternatively, we're also on PIP:

python -m pip install heartpy

That's it! Note that Github always has the newest version.

Documentation

The official documentation is online! You can find the official documentation here

Python 2.7

The module compiles and and runs fine on Python 2.7, but the some unit tests fail.

Tutorial notebooks are now available in Examples/

These show how to handle various analysis tasks with HeartPy, from smartwatch data, smart ring data, regular PPG, and regular (and very noisy) ECG. The notebooks sometimes don't render through the github engine, so either open them locally, or use an online viewer like nbviewer.

We recommend you follow the notebooks in order:

• 1. Analysing a PPG signal, a notebook for starting out with HeartPy using built-in examples.
• 2. Analysing an ECG signal, a notebook for working with HeartPy and typical ECG data.
• 3. Analysing smartwatch data, a notebook on analysing low resolution PPG data from a smartwatch.
• 4. Analysing smart ring data, a notebook on analysing smart ring PPG data.
• 5. Analysing noisy ECG data, an advanced notebook on working with very noisy ECG data, using data from the MIT-BIH noise stress test dataset.
• 6. Colorblind mode - How To and Styles

More information

HeartPy, the Python Heart Rate Analysis Toolkit is a module for heart rate analysis in Python. It started as pure-python implementation to analyse physiological data taken in naturalistic driving and cycling experiments.

The module takes a discrete heart rate signal and outputs time-domain and frequency-domain measures often found in scientific literature:

Time domain:

• beats per minute, BPM
• interbeat interval, IBI
• standard deviation if intervals between adjacent beats, SDNN
• standard deviation of successive differences between adjacent R-R intervals, SDSD
• root mean square of successive differences between adjacend R-R intervals, RMSSD
• proportion of differences between R-R intervals greater than 20ms, 50ms, pNN20, pNN50
• median absolute deviation, MAD
• Poincare analysis (SD1, SD2, S, SD1/SD2)
• Poincare plotting

Frequency domain (ranges per Shaffer and Ginsberg: https://doi.org/10.3389/fpubh.2017.00258)

• very low frequency component (0.0033–0.04 Hz), VLF
• low frequency component (0.04–0.15 Hz), LF
• high frequency component (0.15–0.4 Hz), HF
• lf/hf ratio, LF/HF

When using the package in your research, please cite:

van Gent, P., Farah, H., van Nes, N., & van Arem, B. (2019). Analysing Noisy Driver Physiology Real-Time Using Off-the-Shelf Sensors: Heart Rate Analysis Software from the Taking the Fast Lane Project. Journal of Open Research Software, 7(1), 32. DOI: http://doi.org/10.5334/jors.241

van Gent, P., Farah, H., van Nes, N., & van Arem, B. (2019). HeartPy: A novel heart rate algorithm for the analysis of noisy signals. Transportation Research Part F: Traffic Psychology and Behaviour, 66, 368–378. https://doi.org/10.1016/j.trf.2019.09.015

Documentation

You can find the official documentation here

The module is also to some extent described in my tutorial series:

• Analyzing a Discrete Heart Rate Signal Using Python - Part 1
• Analyzing a Discrete Heart Rate Signal Using Python - Part 2
• Analyzing a Discrete Heart Rate Signal Using Python - Part 3
• Analyzing a Discrete Heart Rate Signal Using Python - Part 4: in development

License

The module is licensed under the MIT License

Validation

Initial results of the validation have been reported in [1, 2].

[1]van Gent, P., Farah, H., van Nes, N., & van Arem, B. (2018). Heart Rate Analysis for Human Factors: Development and Validation of an Open Source Toolkit for Noisy Naturalistic Heart Rate Data. In Proceedings of the 6th HUMANIST Conference (pp. 173–178).

[2] van Gent, P., Farah, H., van Nes, N., & van Arem, B. (2019). HeartPy: A novel heart rate algorithm for the analysis of noisy signals. Transportation Research Part F: Traffic Psychology and Behaviour, 66, 368–378. https://doi.org/10.1016/j.trf.2019.09.015

To-do

The module is still in active development. See the changelog for past changes. The to-do for the coming months is:

to do before V1.3

• Same but for PPG - morphology too variable, method unstable
• Add 'strictness parameter' to affect how HeartPy evaluates peaks for acceptance/rejection
• Add method to handle NaN data automatically
• clean_rr method now removes incorrect values, update to allow for replacement by median of surrounding data points
  • add method that can fill in missing R-peaks, settable to search for either local optimum OR mean imputation.
• Report validation performance on repo (published paper + key-points document once published)
• Change backend structure in anticipation of GUI development
• Develop GUI for HeartPy

heartrate_analysis_python's Issues

BPM of every R peaks and time of this bpm?

Hi!
I'm a university student. In my dissertation, an ECG displaying and processing program must be created. Something similar as in this video:
https://www.youtube.com/watch?v=9jKccn9JaGc

I found the HeartPy yesterday and thought the HeartPy could help with my work.
There was already a similar question: #7
But I’m not sure the solution there is right for my problem. The solution to that problem has been the proccess_segmentwise method, which is known to split the ecg into segments and the width of the segment is given in seconds.

The BPM value is calculated from two RR peaks. So it's okay to split it into segments based on seconds using the process_segmentwise () method, but this only gives the bpm of the segments. I want to get a bpm value after each RR peak. As in this picture:

After all, I want to split it segments by RR peaks and get the bpm values of these RR segments and at what point in time I get these bpm values, which in this case is equal to at what point in time a peak R occurs.

I hope understable what I would like to do. I don't know if this is feasible.

Another problem is that the segment_plotter () method doesn't work for me, but I'm sure I'm setting something wrong.
I tried with this example:

import heartpy as hp

data, timer = hp.load_exampledata(2)
sample_rate = hp.get_samplerate_datetime(timer, timeformat='%Y-%m-%d %H:%M:%S.%f')
working_data, measures = hp.process_segmentwise(data, sample_rate, segment_width = 40, segment_overlap = 0.25)

hp.segment_plotter(working_data,measures)

And I get this error:

Traceback (most recent call last):
  File "F:/.../.../.../01.py", line 7, in <module>
    hp.segment_plotter(working_data,measures)
  File "F:\...\...\...\heartpy\visualizeutils.py", line 198, in segment_plotter
    p = plotter(wd_segment, m_segment, show=False)
  File "F:\...\...\...\heartpy\visualizeutils.py", line 80, in plotter
    fs = working_data['sample_rate']
KeyError: 'sample_rate'

only integer arrays error

When i give the input HR data to the library, i got an error:
peaksy = hrdata[np.where((hrdata > rol_mean))[0]]

TypeError: only integer arrays with one element can be converted to an index

My code is:

import bitalino
import numpy as np
import time
import matplotlib.pyplot as plt
import peakdetect
import heartbeat as hb

Mac OS

macAddress = "/dev/tty.BITalino-XX-XX-DevB"

Windows

#macAddress = "20:16:07:18:13:64"
macAddress= "20:17:09:18:46:90"

#device = bitalino.BITalino(macAddress)
time.sleep(1)

srate = 1000
nframes = 100
threshold = 5

data=np.loadtxt('baselineHr.txt')

ECG=data[:, -1]

ECGmean=np.mean(ECG)

peaks = peakdetect.peakdetect(np.array(EDA), lookahead=2,delta=threshold)
indexes = []
for posOrNegPeaks in peaks:
for peak in posOrNegPeaks:
indexes.append(peak[0])

peaksECG = peakdetect.peakdetect(np.array(ECG),delta=ECGmean/20)
indexesECG = []
for posOrNegPeaks in peaks:
for peak in posOrNegPeaks:
indexesECG.append(peak[0])

ECG_data_arr=[]

for i in range(0,len(ECG)):
ECG_data_arr.append(np.int(ECG[i]))

sp=np.fft.fft(ECG_data_arr)

freq = np.fft.fftfreq(len(ECG_data_arr), d=1/srate)

data = ECG_data_arr
ECG_data.xlsx

measures=hb.process(data,1000)

allow control over what happens to final bin from make_windows and process_segmentwise

I have a bit of 'spare' data at the end of each recording so I don't mind a larger bin at the last index hence the period set to [0,1] position, but it would be nice to have some control over this through a parameter

Originally posted by @KrisITB in #7 (comment)

return indexes of useg segments from process_segmentwise

it works for me but I think it will help others who seek 'plug and play' functionality if there would be some return for indexes to map continuous values to the indexes in the original dataset

Originally posted by @KrisITB in #7 (comment)

spectral.py user warning

When using heartpy v1.2.4 with python v3.7.4 on Windows 64bit to analyse heart rate from PPG data with the hp.process_segmentwise method I get the following warning from scipy\signal:

..\spectral.py:1969: UserWarning: nperseg = 100000 is greater than input length = 38779, using nperseg = 38779 .format(nperseg, input_length))

Interestingly I get this warning on about 80% of the returns, not all of them.
My data is roughly half an hour long, my sample rate is 204.8Hz and my function call is as follows:
working_data, measures = hp.process_segmentwise(matlabdata, sample_rate=fs, segment_width = 40, segment_overlap = 0.25, calc_freq=True, reject_segmentwise=True, report_time=True)

temp containers nn20, nn50 incorrectly appended to output measures dict

nn20 and nn50 returns lists of values e.g,:

Originally posted by @KrisITB in #7 (comment)

BadSignalWarning on data with 0.5s intervals

Hello!

I'm trying to use your tool for the first time.

I took the data from here: http://ecg.mit.edu/time-series/ ("series 2")
Direct link for downloading data: http://ecg.mit.edu/time-series/hr.7275

From the description:

Each series contains 1800 evenly-spaced measurements of instantaneous heart rate from a single subject. The two subjects were engaged in comparable activities for the duration of each series. The measurements (in units of beats per minute) occur at 0.5 second intervals, so that the length of each series is exactly 15 minutes

I'm trying to use this code to load the data:

data = pd.read_csv('hr.7257', header=None, index_col=None, sep='\n')
data.columns = ['data', ]
data = df.data

fs = hp.get_samplerate_mstimer([(x/2) * 1000 for x in data.index.tolist()])
print(fs)  # according to http://ecg.mit.edu/time-series/, is sampled at 0.5[s] intervals
working_data, measures = hp.process(data.values, fs)

The sample rate is 2.0 Hz, but I get the BadSignalWarning. What I am doing wrong? Thanks for help!

Best regards,
Oskar

scale_sections function bug

Dear HeartPy team,

I am a maker. I am designing a PPG wearable device for real-time HRV analysis, data collecting and visualization.

When developing the PC application, I integrated your library. When working on preprocessing for rolling data, I find that "scale_sections" function issue may analysis failed randomly.

Let me share what I've found:

output = np.empty(len(data))
// length of output array is the same as length of data array.

so if data_end is not equal to length of data. Then some data at the end of output numpy array will not be assigned. They will be random value. If one of them is nan, then the process fuction will crash.

To fix this:
just change return np.array(output) -> return np.array(output[0:data_start])

Let me know if I described the issue clearly.
Feel free to let me know your thoughts.

Best,
Lauren

error

hi, I was doing Analysing a Noisy ECG Signal but when i run the code below I'm getting the error
"AttributeError: module 'HeartPy' has no attribute 'enhance_ecg_peaks'"....can u help me with it

filtered =hp.enhance_ecg_peaks(hp.scale_data(ecg), sample_rate,
aggregation='median', iterations=5)

#show filtered signal
plt.figure(figsize=(12,4))
plt.plot(filtered)
plt.show()

#zoom in on signal section and overlay filtered segment
plt.figure(figsize=(12,4))
plt.title('original signal zoom in')
plt.plot(hp.scale_data(ecg[15000:17000]), label='original data')
plt.title('processed signal zoom in')
plt.plot(hp.scale_data(filtered[15000:17000]), alpha=0.5, label='processed data')
plt.legend()
plt.show()

Unique data

Hi Paul,

Nice articles and source code. Is the heart data unique for each user?

Please explain which data points are unique for each user.

Regards
Shafi

Breathing rate return "nanHz"

hi Paul
i try to estimate respiration rate using data from link below. but it return nanHz breathing rate. any suggestion to overcome this.

http://peterhcharlton.github.io/RRest/synthetic_dataset.html

---------------CODE start--------------
import heartpy as hp
data = hp.get_data('rrest-syn001_data.csv', delim=',', column_name='ppg')
fs = 500
working_data, measures = hp.process(data, fs, report_time=True)
print('breathing rate is:%sHz' %measures['breathingrate'])

Threshold algorythm

If the first stage of signal have the distortions, the second (good) stage fully marked as rejected
bad_samp.pdf
good_samp.pdf

Integrate PPG methods with NeuroKit

Hi @renatosc,

I'd heard of it but have not used it so far. The repo looks interesting indeed! Right now I'm working on adding an ECG-specific peak finder to HeartPy based on an often used template matching method.

I'll have a look later today at NeuroKit. The integration between our two packages should be relatively simple. Might indeed be nice to use this to add PPG capability to NeuroKit :)

Cheers,
Paul

Originally posted by @renatosc in #7 (comment)

square_spectrum keyword argument not used

As flagged in #42, requesting a squared frequency spectrum will return the regular spectrum.

frequency estimations

Are the frequency estimations computed from the timeseries signal, or by an interpolated signal constructed from the peak R-R annotations?

plot_poincare() method error

I'm getting an error when I try to call the method plot_poincare(): "alpha must be a float or None".
This error also occurs in the examples like 2_regular_ECG and 5_noisy_ECG of the Github page.

Issue #44 seems to have fixed this but I still get the error.
My HeartPy Version is 1.2.6

process_segment_wise not working since update

Hi Paul,

I seem to be having some serious problems with the process_segment_wise since the latest update...

dataset:
01-05-19_00-25-52.txt

I have a dataset of which I'm using the first 500000 readings. My sampling speed is 1000Hz so that should correspond to the following segments when running:

working_data_segment, measures_segment = hp.process_segmentwise(data, sample_rate=1000.0, segment_width = 300, segment_overlap = 0.9)

Expected results:
[(0, 300000),(30000, 330000),(60000, 360000),(90000, 390000),(12000, 420000),(150000, 450000),(180000, 480000),(210000, 500000)

However - it produces the following segments:

[(0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (0, 300000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (29999, 330000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (59999, 360000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (89999, 390000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999, 420000), (119999,420000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (149999, 450000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (179999, 480000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000), (209999, 500000)]

It worked well before the latest update, not sure what's going wrong now. I tested it on another dataset and got similar results

Plot of data2 with exact same code is a different (and peak detection is worse) than in documentation

Hi Paul,

I tried using your documentation (https://buildmedia.readthedocs.org/media/pdf/python-heart-rate-analysis-toolkit/latest/python-heart-rate-analysis-toolkit.pdf) to make the plot using data2 (pp. 11-12). However, when I plot it (using the same exact code as in the documentation), I get the following plot which is clearly different from the one in the documentation:

Picture from the documentation:

In my version, the rolling average is clearly visible (and seems rather high), and, thus, it misses some peaks that your version of the plot includes.

Could this be due to something in the toolkit? Or maybe I am doing something wrong?

I downloaded and ran it today.

Thank you in advance!

Best,

SigneKB

Only working fine on the 14500:20500,

It's only working on the range given, for other ranges, it's not working.

KrisITB

If running make_windows() on the data attached with the default values that you have above (windowsize = 120 and sample rate =100) make_windows() return 11 pairs of int values but the last one is 72000 while my sample has 77953 records

This is correct. For now the method computes only on full windows and discards the final (partial) window. I'll flag this as a new issue and will add something and push it tomorrow.

Originally posted by @paulvangentcom in #7 (comment)

C/C++ code status

I read somewhere (maybe in the tutorials?) that you would be porting this code to c/c++. What is the status of that?

Sample Rate for Frequency Domain Calculation

Dear Paul:

I'm not sure whether this is an issue or not, but I found that in the function calc_fd_measures(), there no sample rate as a input variable. And for all three methods (fft, periodogram and welch), the sample rate is fixed as 1000. Could you please tell me the reason for this?

Thanks so much for your help!

attribute missing

hi, I have downloaded the latest version=1.2.4 of heartpy package and still I'm getting this error:

filtered = hp.enhance_ecg_peaks(hp.scale_data(ecg), sample_rate, iterations=4, aggregation='median', notch_filter=True)

AttributeError: module 'HeartPy' has no attribute 'enhance_ecg_peaks

I appreciate your help!!!!😄

Visualisation of rejected RR values

Hey!

It would be super useful to be able to visualise the rejected RR-interval values (RR_list_cor) on the graph - to make sure that the intervals are being correctly rejected and for other troubleshooting issues.

Is there a function for this in the code? I couldn't find one..

Thanks! :)

Plot visualisation when using hp.process_segmentwise

Hi Paul,

When using the hp.process_segmentwise() function I have not been able to figure out how to plot separate graphs for each segment (other than separating the data beforehand).
Is there a way to do this?

Cheers!

fix issue where poincare plotting fails

alpha passed as string rather than float.

Access to continous bpm?

Hi!

First: congratulations! your ppg to hr conversion seems to be handling my data better than the application that came with the sensor %) and many thanks for the tutorial on your website!

Quick question: is there a way to access BPM values as they change over time (ideally averaged over a specified window)? From what I understand measures['bpm'] only returns the average, or am I missing something?

Am I correct to understand that I need to use indexes from working_data['peaklist_cor'] to calculate this myself? If so please consider this a feature request ;) I'm sure I don't have to stress how valuable such feature is in research practice and you already have an awesome tool here so that little feature would certainly make this the hottest stuff in town IMHO... at least for those of us cursed with working with ppg ;)

P.s.: any tips on how large the window should be?

Best wishes!

Add info for output measures to API

It would be helpful for better documentation regarding function outputs. For instance, calc_fd_measures generates the ouput 'measures' dict but there is no documentation as to what units LF and HF are given in, or whether LF and HF are peak or power measures.
If this documentation is already there, then it is not easy to find.

segment_plotter fails due to export type error

When using heartpy v1.2.4 with python v3.7.4 on Windows 64bit to visualise segmented heartrate estimation data called using the following statements:
working_data, measures = hp.process_segmentwise(matlabdata, sample_rate=fs, segment_width = 40, segment_overlap = 0.25, calc_freq=True, reject_segmentwise=True, report_time=True) segpath = 'E:\GIT\PhD\Experiment1\TSAnalysis\data\shimmer\segments' hp.segment_plotter(working_data, measures, title='Heart Rate Signal Peak Detection', path=segpath, start=0, end=None, step=1)

The following error occurs:
File "......\Python\Python37\lib\site-packages\matplotlib\backend_bases.py", line 1956, in _get_output_canvas .format(fmt, ", ".join(sorted(self.get_supported_filetypes())))) ValueError: Format 'jpg' is not supported (supported formats: eps, pdf, pgf, png, ps, raw, rgba, svg, svgz)
The error is caused by line 169 in visualizeutils.py:
p.savefig('%s%i.jpg' %(path, filenum))
As the error states, matplotlib does not support exporting the jpeg format with version of 3.1.1 the tools.

The error can be fixed by using one of the supported file types, such as png:
p.savefig('%s%i.png' % (path, filenum))
Resulting in the expected output.

Data quality question

Hi,

I'm experimenting with using your very nice toolbox for a project where I want to remotely extract heartrate variability. In my current setup I'm limited to a 30Hz/fps sampling rate - and I think this might be a cause for the issue I'm experiencing.

From my dataset (rawHR) I get quite a few false negatives I wouldn't expect.

For example from samples 0 to 1000 (0 tot 33 seconds)

where on other spots it is detected fine

I'm using this datafile (as a .csv but I can't upload that here) which is a 4minute snippet from a video where the PPG is extracted from the green channel:
hr sample.txt

and preprocessing it using your butterworth 'filtersignal(data,5,30,3)' and 'scale_data' functions on it.

What would be your advice for better results? Upsample, more filtering or is the signal quality too poor.

Thank's in advance!

ValueError: min() arg is an empty sequence

Hey. I am running this code but it is throwing error.
import pandas as pd
import heartpy as hp
df = pd.read_csv('E:/Autistic_data/drive03'+'.csv')
working_data, measures = hp.process(df.iloc[:,5], 100, calc_freq=True)

ValueError Traceback (most recent call last)
in
1 df = pd.read_csv('E:/Autistic_data/drive4'+'.csv')
----> 2 working_data, measures = hp.process(df.iloc[:,5], 100, calc_freq=True)

~\AppData\Roaming\Python\Python37\site-packages\heartpy\heartpy.py in process(hrdata, sample_rate, windowsize, report_time, calc_freq, freq_method, interp_clipping, clipping_scale, interp_threshold, hampel_correct, bpmmin, bpmmax, reject_segmentwise, measures, working_data)
647 working_data['hr'] = hrdata
648 rol_mean = rolmean(hrdata, windowsize, sample_rate)
--> 649 working_data = fit_peaks(hrdata, rol_mean, sample_rate, working_data=working_data)
650 working_data = calc_rr(working_data['peaklist'], sample_rate, working_data=working_data)
651 working_data = check_peaks(working_data['RR_list'], working_data['peaklist'], working_data['ybeat'],

~\AppData\Roaming\Python\Python37\site-packages\heartpy\heartpy.py in fit_peaks(hrdata, rol_mean, sample_rate, bpmmin, bpmmax, working_data)
377
378
--> 379 working_data['best'] = min(valid_ma, key=lambda t: t[0])[1]
380 working_data = detect_peaks(hrdata, rol_mean, min(valid_ma, key=lambda t: t[0])[1],
381 sample_rate, update_dict=True, working_data=working_data)

ValueError: min() arg is an empty sequence

Below is the link for dataset I have used.

https://drive.google.com/open?id=1e08C3kyeq0kkuFdhQ9J7a0vUsZRGY0D8

Crash on noicy signal

1.txt
In a noicy signal in fit_peaks function if calculated _bpm always upper than bpmmax (180 bpm), then valid_ma is empty list and I have error "min() arg is an empty sequence"

Accept .bdf file?

Can this program accept .bdf file as inputs?

fix issue where multiple analyses without extracting values fail to overwrite each other

working_data and measures dicts seem to not be replaced properly

AttributeError: module 'HeartPy' has no attribute 'remove_baseline_wander'

I am trying to run noise removal example code. I have installed heatypy using

python -m pip install heartpy
successfully. Still I am getting this error message

AttributeError                            Traceback (most recent call last)
<ipython-input-3-41d5b2afa0be> in <module>
     23     return filtered
     24 
---> 25 filtered = filter_and_visualise(ecg, sample_rate)

<ipython-input-3-41d5b2afa0be> in filter_and_visualise(data, sample_rate)
      5     '''
      6 
----> 7     filtered = hp.remove_baseline_wander(data, sample_rate)
      8 
      9     plt.figure(figsize=(12,3))

AttributeError: module 'HeartPy' has no attribute 'remove_baseline_wander'

Multiple channels?

Hi,

Thank you for this awesome library.

I was wondering if there's a way to analyze multiple channels. I've seen some papers where linear combinations, averages of different channels are taken to reduce noises or motion artifacts. For instance when recording raw ppg from smartwatches, or arm bands such as Polar OH, we get values from 3 channels. Is there anything in the library that handles such multi-channel inputs?

I'd appreciate if you could please suggest what's currently the go-to approach for handling raw ppg signals from different channels? Visualizing the signals and only using the less noisy one or averaging all or something else altogether. Any insight into this would really be helpful.

Thank you!

Frequency Analysis Result

First, I sincerely appreciate for the properly working BVP/ECG analysis codes you shared.
I have a question about the frequency domain result.
When I try to analyze BVP data with very long length (>3hrs),
The result shows EXTREAMLY large value for hf and lf power (line 44).
But the lf/hf makes sense, and the value is reasonable if I divide it by 1000^2 (maybe for ms^2).
Would you let me know what is the unit of raw HF, LF and is it correct to divide it by 1000^2 for get the trivial ms^2 result?

some output measures are different in v1.2.6

Hi! I was running through the smart watch data and following your Analysing Smartwatch Data notebook but some of the output values differ a lot such as the sdsd and the rmssd, hoping there could be some analysis as to why or what these numbers would mean.

From the notebook:
bpm: 78.981709
ibi: 759.669563
sdnn: 79.989043
sdsd: 99.586554
rmssd: 99.586688
pnn20: 0.385797
pnn50: 0.205374
hr_mad: 41.049817
sd1: 36.948123
sd2: 73.723698
s: 8557.548377
sd1/sd2: 0.501170
breathingrate: 0.239659

My output:
bpm: 75.543450
ibi: 794.244901
sdnn: 92.544456
sdsd: 33.387928
rmssd: 49.947643
pnn20: 0.644444
pnn50: 0.244444
hr_mad: 35.056933
sd1: 34.389733
sd2: 85.202755
s: 9205.180571
sd1/sd2: 0.403622
breathingrate: 0.200000

Customize the plotter window length

If you add to the plotter function a parameter length = 6 and the corresponding plt.figure(figsize=(length ,6)) to the function, one can dynamically adjust the length of the resulting pyplot (for example: plotter(df,20) )plot in inches.
Something that might be handy for quickly viewing the result of a larger set.

Accessing IBI intervals

Hi Paul!
Great work with this toolbox :) I've been using it with my Pulse Sensor and seem to be getting some very good results on the large part.
However, I have had a couple of problems while attempting to collect some results:

1 - List of IBIs from detected peaks
I have run the peak detection functions 'hp.enhance_peaks' and 'hp.process' functions on my data. I would now like to access the list of IBI values calculated from all the consecutively accepted peaks (provided no peaks were rejected in between the accepted peaks). However when accessing measures['RR_list_cor'] and measures['RR_list'] all the values produced seem to be rounded to the nearest 20mS (values such as 860, 880, 900 etc..). Would it be possible to the nearest mS instead?

2- Varying signal amplitude
My sensor signal amplitude varies quite a lot over time. To correct this, I used the enhanced_peaks functions which improved the results. However, I am still having some issues with incorrectly rejected peaks and undetected peaks. Images included below. Are there any other functions I could be using to rectify this? (I have included a screenshot of both a well-performing sample and one where the problems occur below)

Thanks for the help! I look forward to using this toolbox and hopefully contributing in the future :)

Francois

AVR simple logger with scaling + python peak detection

Hi Paul,

I've implemented the AVR simple logger with scaling with my Arduino UNO this morning.
I then used this python toolbox for peak detection (the Arduino version is limited to 250Hz whereas I am sampling at 1000Hz for 1ms RR-interval precision). The scaler-peak detection combination seems to work perfectly until it encounters some noise. I've tried to illustrate this with the pictures below:

This represents the first 100k data entries (sampled at 1kHz). Peak detection working perfectly!

At approx ~110k samples, the sensor encounters a little bit of noise (probably due to movement) leading to two problems:
--> The systolic and diastolic are both detected - doubling the heart rate
-->As a result the RR intervals are incorrect

My understanding is that this happens because, following the noise, the moving average moves below the diastolic peaks.

Graph showing the first 200k samples, with the noise @110k.

Here's my data file:

sensor_scaled_data.txt

cant install python -m pip install heartpy(python -m pip install heartpy)

I can install in base environment but not in virtual envs
i have tried to install heartpy on windows gpu in virtual environment with pytorch 1.6 but its showing SSL error although i can install it in base environment
in virtual env i have to specify the python version instead of - m I think

Bug on setup.py

Hey Paul. In your last commit, you updated the version to 1.1.5 but forgot to close the double quotes:

setup.py line 8 right now is:

version="1.1.5,

dfitpack error from process_segmentwise

When using heartpy v1.2.4 with python v3.7.4 on Windows 64bit to analyse heart rate from PPG data with the hp.process_segmentwise method I get the following error from the scipy interpolate lib:
dfitpack.error: (m>k) failed for hidden m: fpcurf0:m=3
Following a debug trace, the error occurs from the call of calculate_fd_measures in analysis.py line 426:
interpolated_func = UnivariateSpline(rr_x, rr_list, k=3)
which in turn calls where the error originates in scipy:
data = dfitpack.fpcurf0(x, y, k, w=w, xb=bbox[0], xe=bbox[1], s=s)
In the case that the error might be caused by the same issues as the following thread:
https://stackoverflow.com/questions/32230362/python-interpolate-univariatespline-package-error-mk-failed-for-hidden-m

I have checked the the lengxth of x is 3, the length of y is 3 and k = 3, w bbox[0] bbox[1] and s = None

The value of x is

The value of rr_list is:

Reducing the order of the spline to 2 returns successfully, but 3 fails in this case. According to the linked stackoverflow, k in this case needs to be at least one bigger than the length of x, so it would appear that the mechanism which creates RR_list_cor in the working_data dict is causing an error later in the chain.

interpolating clipping while using process_segmentwise

Hi Paul,

Is there a way of interpolating clipping while using
process_segmentwise(data)

I know it's possible when using proces sby setting
process(data,interp_clipping=False)

Cheers,
Francois

The function get_samplerate_datetime needs to be modified. It should accept datetime input as well (or only imo).

As referenced in #29 by @J535D165

Function to be expanded to accept both datetime objects as well as strings.

not able to use Notch filter

hi Paul

im using heartpy v1.2.6. i try to apply notch filter using signal_filter and remove_baseline_wander. the code as below. i attached sample data (1.txt). when i apply lowpass filter the code run without issue but when i try to apply notch filter using remove_baseline_wander function it pop up error as below

when i try to apply notch filter using signal_filter function the error as below

`import numpy as np
import heartpy as hp
import matplotlib.pyplot as plt

input = '1.txt'
dataPPG = numpy.loadtxt(input, delimiter=",", usecols=0)
plt.plot(dataPPG, label='Raw_data')
plt.legend()
plt.show()

########################baseline wander###########################
bw = hp.remove_baseline_wander(dataPPG, sample_rate= 100, cutoff=60.0)
plt.plot(bw, label='Raw_data')
plt.show()`

########################## signal_filter##########################
nf = hp.filter_signal(dataPPG, cutoff=60,sample_rate=100,filtertype='notch')
plt.plot(nf, label='notch filter')
plt.show()

#########################signal_filter_lowpass####################

lp = hp.filter_signal(dataPPG, cutoff=0.05,sample_rate=100,filtertype='lowpass')
plt.plot(nf, label='Lowpass filter')
plt.show()

Thanks in advance.

1.txt

Sampling Frequency of 44100 Hz ValueError: min() arg is an empty sequence

Hi

I am running into an issue with using this module to calculate the BPM of a .wav file sampled at 44100Hz (Python 3.6)

This is the error which I get: 
 Traceback (most recent call last):
  File "/Users/JohnDoe/Desktop/TensorFlow_Programs/Heart-Sounds-Deep-Learning/heart-rate-analysis-module/UsingHeartbeatModule.py", line 15, in <module>
    measures = hb.process(data, 44100)
  File "/Users/JohnDoe/Desktop/TensorFlow_Programs/Heart-Sounds-Deep-Learning/heart-rate-analysis-module/heartbeat.py", line 205, in process
    fit_peaks(hrdata, rol_mean, fs)
  File "/Users/JohnDoe/Desktop/TensorFlow_Programs/Heart-Sounds-Deep-Learning/heart-rate-analysis-module/heartbeat.py", line 124, in fit_peaks
    working_data['best'] = min(valid_ma, key = lambda t: t[0])[1]
ValueError: min() arg is an empty sequence

I am using your heart-rate module to calculate the BPM (shown below):

import heartbeat as hb
from scipy.io import wavfile as wav
import numpy
import pandas as pd


RAW_CSV_NAME = "/Users/JohnDoe/Desktop/TensorFlow_Programs/Heart-Sounds-Deep-Learning/pls.csv"
fs, data = wav.read('/Users/JohnDoe/Desktop/output_test.wav')
df = pd.DataFrame(data)
df.columns=['hart']
df.to_csv(RAW_CSV_NAME, index=False, header="hart" )

print("Sampling frequency " + str(fs)) --> prints 44100
hrdata = hb.get_data(RAW_CSV_NAME, column_name = 'hart')
measures = hb.process(data, 44100)

print(measures['bpm']) #returns BPM value
print(measures['lf/hf']) # returns LF:HF ratio

#Alternatively, use dictionary stored in module:
print(hb.measures['bpm']) #returns BPM value
print(hb.measures['lf/hf']) # returns LF:HF ratio

Here is the code for recording with PyAudio

import pyaudio
import wave

CHUNK = 1024
FORMAT = pyaudio.paInt16
CHANNELS = 1
RATE = 44100
RECORD_SECONDS = 30
WAVE_OUTPUT_FILENAME = "/Users/sreeharirammohan/Desktop/output_test.wav"

p = pyaudio.PyAudio()

stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)

print("* recording")

frames = []

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    frames.append(data)

print("* done recording")

stream.stop_stream()
stream.close()
p.terminate()

wf = wave.open(WAVE_OUTPUT_FILENAME, 'wb')
wf.setnchannels(CHANNELS)
wf.setsampwidth(p.get_sample_size(FORMAT))
wf.setframerate(RATE)
wf.writeframes(b''.join(frames))
wf.close()

I was wondering if you could help me fix this Sampling Rate Issue. Any Guidance would be greatly appreciated!

Hi Sreehari, how can i import heartbeat_audio_experimental in python

Hi Sreehari,

Apologies for the delay, I broke my thumb earlier so I've not had much coding time available.

I've pushed an experimental function to the repo (in the folder 'experimental'). Please be aware that the peak detection is not super accurate yet and will be thrown off if there's a lot of noise.

Since the variability measures are noise-sensitive they should not be trusted. The BPM value as you requested, however, is much more noise-resistant and should be a pretty good estimate of the real BPM in the audio file (given there is not too much noise).

Usage is similar to the main module:

import heartbeat_audio_experimental as hb_audio

measures = hb_audio.process('heartbeat.wav')
print(measures['bpm'])
hb_audio.plotter()

This gives:

Let me know if this helps your project.

• Paul van Gent

Originally posted by @paulvangentcom in #2 (comment)

Welch and periodogram methods in calc_fd_measures

Hi,
In calc_fd_measures, if using welch or periodogram method, it is needed to add dx = freq[1]-freq[0] parameters for np.trapz function for calculating HF and LF. It works for fft since you divide fft by datalen.