diff --git a/heartpy/analysis.py b/heartpy/analysis.py index 41e2cd3..fbf0913 100644 --- a/heartpy/analysis.py +++ b/heartpy/analysis.py @@ -538,7 +538,8 @@ def calc_fd_measures(method='welch', square_spectrum=True, measures={}, working_ return working_data, measures -def calc_breathing(rrlist, hrdata, sample_rate, measures={}, working_data={}): +def calc_breathing(rrlist, hrdata, sample_rate, method='fft', + measures={}, working_data={}): '''estimates breathing rate Function that estimates breathing rate from heart rate signal. @@ -582,28 +583,42 @@ def calc_breathing(rrlist, hrdata, sample_rate, measures={}, working_data={}): Breathing is then computed with the function - >>> m = calc_breathing(wd['RR_list_cor'], data, sample_rate = 100.0, measures = m, working_data = wd) + >>> m, wd = calc_breathing(wd['RR_list_cor'], data, sample_rate = 100.0, measures = m, working_data = wd) >>> round(m['breathingrate'], 3) - 0.161 + 0.128 There we have it, .16Hz, or about one breathing cycle in 6.25 seconds. ''' + #resample RR-list to 1000Hz x = np.linspace(0, len(rrlist), len(rrlist)) - x_new = np.linspace(0, len(rrlist), len(rrlist)*10) + x_new = np.linspace(0, len(rrlist), np.sum(rrlist)) interp = UnivariateSpline(x, rrlist, k=3) breathing = interp(x_new) - breathing_rolling_mean = rolling_mean(breathing, 0.75, sample_rate) - peaks, working_data = hp.peakdetection.detect_peaks(breathing, breathing_rolling_mean, 1, sample_rate, - update_dict=False) - - if len(peaks) > 1: - signaltime = len(hrdata) / sample_rate - measures['breathingrate'] = len(peaks) / signaltime + + if method.lower() == 'fft': + datalen = len(breathing) + frq_ = np.fft.fftfreq(datalen, d=((1/1000.0))) + frq_ = frq_[range(int(datalen/2))] + Y = np.fft.fft(breathing)/datalen + Y = Y[range(int(datalen/2))] + psd_ = np.power(np.abs(Y), 2) + elif method.lower() == 'welch': + frq_, psd_ = welch(breathing, fs=1000, nperseg=len(breathing)) else: - measures['breathingrate'] = np.nan # pragma: no cover + raise ValueError('Breathing rate extraction method not understood! Must be \'welch\' or \'fft\'!') - return measures + #take out lowest peak + frq = frq_[frq_ >= 0.04] + psd = psd_[frq_ >= 0.04] + + #find max + measures['breathingrate'] = frq[np.argmax(psd)] + working_data['breathing_signal'] = breathing + working_data['breathing_psd'] = psd + working_data['breathing_frq'] = frq + + return measures, working_data def calc_poincare(rr_list, rr_mask=[], measures={}, working_data={}): diff --git a/heartpy/heartpy.py b/heartpy/heartpy.py index c0a2ba9..3b47979 100644 --- a/heartpy/heartpy.py +++ b/heartpy/heartpy.py @@ -20,7 +20,7 @@ remove_baseline_wander, smooth_signal from .peakdetection import make_windows, append_dict, fit_peaks, check_peaks, \ check_binary_quality, interpolate_peaks -from .visualizeutils import plotter, segment_plotter, plot_poincare +from .visualizeutils import plotter, segment_plotter, plot_poincare, plot_breathing from .analysis import calc_rr, calc_rr_segment, clean_rr_intervals, calc_ts_measures, \ calc_fd_measures, calc_breathing, calc_poincare @@ -36,6 +36,7 @@ 'hampel_filter', 'load_exampledata', 'plotter', + 'plot_breathing', 'plot_poincare', 'process', 'process_rr', @@ -54,7 +55,7 @@ def process(hrdata, sample_rate, windowsize=0.75, report_time=False, calc_freq=False, freq_method='welch', freq_square=True, interp_clipping=False, clipping_scale=False, interp_threshold=1020, hampel_correct=False, bpmmin=40, bpmmax=180, reject_segmentwise=False, - high_precision=False, high_precision_fs=1000.0, + high_precision=False, high_precision_fs=1000.0, breathing_method='fft', clean_rr=False, clean_rr_method='quotient-filter', measures={}, working_data={}): '''processes passed heart rate data. @@ -134,6 +135,10 @@ def process(hrdata, sample_rate, windowsize=0.75, report_time=False, the sample rate to which to upsample for more accurate peak position estimation default : 1000 Hz + breathing_method : str + method to use for estimating breathing rate, should be 'welch' or 'fft' + default : fft + clean_rr : bool if true, the RR_list is further cleaned with an outlier rejection pass default : false @@ -273,7 +278,9 @@ def process(hrdata, sample_rate, windowsize=0.75, report_time=False, working_data = working_data) try: - measures = calc_breathing(working_data['RR_list_cor'], hrdata, sample_rate, measures=measures) + measures, working_data = calc_breathing(working_data['RR_list_cor'], hrdata, sample_rate, + method = breathing_method, measures=measures, + working_data=working_data) except: measures['breathingrate'] = np.nan diff --git a/heartpy/visualizeutils.py b/heartpy/visualizeutils.py index 016dc13..d841711 100644 --- a/heartpy/visualizeutils.py +++ b/heartpy/visualizeutils.py @@ -12,7 +12,8 @@ __all__ = ['plotter', 'segment_plotter', - 'plot_poincare'] + 'plot_poincare', + 'plot_breathing'] def plotter(working_data, measures, show=True, title='Heart Rate Signal Peak Detection', moving_average=False): # pragma: no cover @@ -201,6 +202,10 @@ def plot_poincare(working_data, measures, show = True, dictionary object used by heartpy to store computed measures. Will be created if not passed to function + show : bool + whether to show the plot right away, or return a matplotlib object for + further manipulation + title : str the title used in the plot @@ -310,4 +315,55 @@ def rotate_vec(x, y, angle): x_rot = (x * cs) - (y * sn) y_rot = (x * sn) + (y * cs) - return x_rot, y_rot \ No newline at end of file + return x_rot, y_rot + + +def plot_breathing(working_data, measures, show=True): # pragma: no cover + '''plots extracted breathing signal and spectrogram + + Function that plots the breathing signal extracted from RR-intervals alongside + its computed spectrogram representation. + + Parameters + ---------- + working_data : dict + dictionary object that contains all heartpy's working data (temp) objects. + will be created if not passed to function + + measures : dict + dictionary object used by heartpy to store computed measures. Will be created + if not passed to function + + show : bool + whether to show the plot right away, or return a matplotlib object for + further manipulation + + Returns + ------- + out : matplotlib plot object + only returned if show == False. + + Examples + -------- + This function has no examples. See documentation of heartpy for more info. + ''' + + plt.subplot(211) + plt.plot(working_data['breathing_signal'], label='breathing signal') + plt.xlabel('ms') + plt.title('breathing signal extracted from RR-intervals') + + plt.subplot(212) + plt.plot(working_data['breathing_frq'], working_data['breathing_psd'], label='spectrogram') + plt.xlim(0, 2) + plt.xlabel('Hz') + plt.title('spectrogram extracted from breathing rate signal') + + plt.legend() + plt.tight_layout() + + if show: + plt.show() + else: + return plt +