improve breathing rate estimation

Breathing rate is now estimated after upsampling array of peak-peak intervals to 1000Hz, then extracting the dominant frequency over 0.04Hz.

This should significantly reduce np.nan occurrences with breathing rate estimation, and increase precision.

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={}):

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
 

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
+    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
+