plot_items: add sampling_freq and ann_freq parameters.

An annotation file (represented by a wfdb.Annotation object) may have
a "sampling frequency" that differs from the sampling frequency or
frequencies of the signals themselves.  This will be the case, for
example, for annotations generated by programs like sqrs that operate
on an upsampled or downsampled copy of the input signals.

In any event, when plotting annotations, we need to translate the
annotation time into a sample number in order to display it in the
correct location.

Furthermore, in the future, we want to permit plotting multiple
synchronized signals that are sampled at different frequencies.

Therefore, to disambiguate between the many possible "sampling
frequencies" invvolved, add parameters 'sampling_freq' and 'ann_freq'.
Both parameters are optional and default to the value of 'fs'.  Either
may be a list (one entry per channel), allowing the API to accomodate
multi-frequency data.  Currently, if 'sampling_freq' is a list, all of
its elements must be equal.

Author: Benjamin Moody

wfdb/plot/plot.py

@@ -8,12 +8,79 @@
 from wfdb.io.annotation import Annotation
 
 
+def _get_sampling_freq(sampling_freq, n_sig, frame_freq):
+    """
+    Convert application-specified sampling frequency to a list.
+
+    Parameters
+    ----------
+    sampling_freq : number or sequence or None
+        The sampling frequency or frequencies of the signals.  If this is a
+        list, its length must equal `n_sig`.  If unset, defaults to
+        `frame_freq`.
+    n_sig : int
+        Number of channels.
+    frame_freq : number or None
+        Default sampling frequency (record frame frequency).
+
+    Returns
+    -------
+    sampling_freq : list
+        The sampling frequency for each channel (a list of length `n_sig`.)
+
+    """
+    if sampling_freq is None:
+        return [frame_freq] * n_sig
+    elif hasattr(sampling_freq, '__len__'):
+        if len(sampling_freq) != n_sig:
+            raise ValueError('length mismatch: n_sig = {}, '
+                             'len(sampling_freq) = {}'.format(
+                                 n_sig, len(sampling_freq)))
+        return list(sampling_freq)
+    else:
+        return [sampling_freq] * n_sig
+
+
+def _get_ann_freq(ann_freq, n_annot, frame_freq):
+    """
+    Convert application-specified annotation frequency to a list.
+
+    Parameters
+    ----------
+    ann_freq : number or sequence or None
+        The sampling frequency or frequencies of the annotations.  If this
+        is a list, its length must equal `n_annot`.  If unset, defaults to
+        `frame_freq`.
+    n_annot : int
+        Number of channels.
+    frame_freq : number or None
+        Default sampling frequency (record frame frequency).
+
+    Returns
+    -------
+    ann_freq : list
+        The sampling frequency for each channel (a list of length `n_annot`).
+
+    """
+    if ann_freq is None:
+        return [frame_freq] * n_annot
+    elif hasattr(ann_freq, '__len__'):
+        if len(ann_freq) != n_annot:
+            raise ValueError('length mismatch: n_annot = {}, '
+                             'len(ann_freq) = {}'.format(
+                                 n_annot, len(ann_freq)))
+        return list(ann_freq)
+    else:
+        return [ann_freq] * n_annot
+
+
 def plot_items(signal=None, ann_samp=None, ann_sym=None, fs=None,
                time_units='samples', sig_name=None, sig_units=None,
                xlabel=None, ylabel=None, title=None, sig_style=[''],
                ann_style=['r*'], ecg_grids=[], figsize=None,
                sharex=False, sharey=False, return_fig=False, 
-               return_fig_axes=False):
+               return_fig_axes=False, sampling_freq=None,
+               ann_freq=None):
     """
     Subplot individual channels of signals and/or annotations.
 
@@ -87,6 +154,14 @@ def plot_items(signal=None, ann_samp=None, ann_sym=None, fs=None,
         'figsize' argument passed into matplotlib.pyplot's `figure` function.
     return_fig : bool, optional
         Whether the figure is to be returned as an output argument.
+    sampling_freq : number or sequence, optional
+        The sampling frequency or frequencies of the signals.  If this is a
+        list, it must have the same length as the number of channels.  If
+        unspecified, defaults to `fs`.
+    ann_freq : number or sequence, optional
+        The sampling frequency or frequencies of the annotations.  If this
+        is a list, it must have the same length as `ann_samp`.  If
+        unspecified, defaults to `fs`.
 
     Returns
     -------
@@ -113,18 +188,25 @@ def plot_items(signal=None, ann_samp=None, ann_sym=None, fs=None,
     # Figure out number of subplots required
     sig_len, n_sig, n_annot, n_subplots = get_plot_dims(signal, ann_samp)
 
+    # Convert sampling_freq and ann_freq to lists if needed
+    sampling_freq = _get_sampling_freq(sampling_freq, n_sig, fs)
+    ann_freq = _get_ann_freq(ann_freq, n_annot, fs)
+
     # Create figure
     fig, axes = create_figure(n_subplots, sharex, sharey, figsize)
 
     if signal is not None:
-        plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes)
+        plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes,
+                    sampling_freq=sampling_freq)
 
     if ann_samp is not None:
         plot_annotation(ann_samp, n_annot, ann_sym, signal, n_sig, fs,
-                        time_units, ann_style, axes)
+                        time_units, ann_style, axes,
+                        sampling_freq=sampling_freq, ann_freq=ann_freq)
 
     if ecg_grids:
-        plot_ecg_grids(ecg_grids, fs, sig_units, time_units, axes)
+        plot_ecg_grids(ecg_grids, fs, sig_units, time_units, axes,
+                       sampling_freq=sampling_freq)
 
     # Add title and axis labels.
     # First, make sure that xlabel and ylabel inputs are valid
@@ -238,7 +320,8 @@ def create_figure(n_subplots, sharex, sharey, figsize):
     return fig, axes
 
 
-def plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes):
+def plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes,
+                sampling_freq=None):
     """
     Plot signal channels.
 
@@ -262,22 +345,39 @@ def plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes):
         will be used for all channels.
     axes : list
         The information needed for each subplot.
+    sampling_freq : number or sequence, optional
+        The sampling frequency or frequencies of the signals.  If this is a
+        list, it must have the same length as the number of channels.  If
+        unspecified, defaults to `fs`.
 
     Returns
     -------
     N/A
 
     """
+    if n_sig == 0:
+        return
+
     # Extend signal style if necessary
     if len(sig_style) == 1:
         sig_style = n_sig * sig_style
 
+    # Convert sampling_freq to a list if needed
+    sampling_freq = _get_sampling_freq(sampling_freq, n_sig, fs)
+
+    if any(f != sampling_freq[0] for f in sampling_freq):
+        raise NotImplementedError(
+            'multiple sampling frequencies are not supported')
+
     # Figure out time indices
     if time_units == 'samples':
         t = np.linspace(0, sig_len-1, sig_len)
     else:
-        downsample_factor = {'seconds':fs, 'minutes':fs * 60,
-                             'hours':fs * 3600}
+        downsample_factor = {
+            'seconds': sampling_freq[0],
+            'minutes': sampling_freq[0] * 60,
+            'hours': sampling_freq[0] * 3600
+        }
         t = np.linspace(0, sig_len-1, sig_len) / downsample_factor[time_units]
 
     # Plot the signals
@@ -289,7 +389,7 @@ def plot_signal(signal, sig_len, n_sig, fs, time_units, sig_style, axes):
 
 
 def plot_annotation(ann_samp, n_annot, ann_sym, signal, n_sig, fs, time_units,
-                    ann_style, axes):
+                    ann_style, axes, sampling_freq=None, ann_freq=None):
     """
     Plot annotations, possibly overlaid on signals.
     ann_samp, n_annot, ann_sym, signal, n_sig, fs, time_units, ann_style, axes
@@ -320,6 +420,14 @@ def plot_annotation(ann_samp, n_annot, ann_sym, signal, n_sig, fs, time_units,
         will be used for all channels.
     axes : list
         The information needed for each subplot.
+    sampling_freq : number or sequence, optional
+        The sampling frequency or frequencies of the signals.  If this is a
+        list, it must have the same length as the number of channels.  If
+        unspecified, defaults to `fs`.
+    ann_freq : number or sequence, optional
+        The sampling frequency or frequencies of the annotations.  If this
+        is a list, it must have the same length as `ann_samp`.  If
+        unspecified, defaults to `fs`.
 
     Returns
     -------
@@ -330,25 +438,45 @@ def plot_annotation(ann_samp, n_annot, ann_sym, signal, n_sig, fs, time_units,
     if len(ann_style) == 1:
         ann_style = n_annot * ann_style
 
-    # Figure out downsample factor for time indices
-    if time_units == 'samples':
-        downsample_factor = 1
-    else:
-        downsample_factor = {'seconds':float(fs), 'minutes':float(fs)*60,
-                             'hours':float(fs)*3600}[time_units]
+    # Convert sampling_freq and ann_freq to lists if needed
+    sampling_freq = _get_sampling_freq(sampling_freq, n_sig, fs)
+    ann_freq = _get_ann_freq(ann_freq, n_annot, fs)
 
     # Plot the annotations
     for ch in range(n_annot):
+        afreq = ann_freq[ch]
+        if ch < n_sig:
+            sfreq = sampling_freq[ch]
+        else:
+            sfreq = afreq
+
+        # Figure out downsample factor for time indices
+        if time_units == 'samples':
+            if afreq is None and sfreq is None:
+                downsample_factor = 1
+            else:
+                downsample_factor = afreq / sfreq
+        else:
+            downsample_factor = {
+                'seconds': float(afreq),
+                'minutes': float(afreq) * 60,
+                'hours': float(afreq) * 3600
+            }[time_units]
+
         if ann_samp[ch] is not None and len(ann_samp[ch]):
             # Figure out the y values to plot on a channel basis
 
             # 1 dimensional signals
             try:
                 if n_sig > ch:
+                    if sfreq == afreq:
+                        index = ann_samp[ch]
+                    else:
+                        index = (sfreq / afreq * ann_samp[ch]).astype('int')
                     if signal.ndim == 1:
-                        y = signal[ann_samp[ch]]
+                        y = signal[index]
                     else:
-                        y = signal[ann_samp[ch], ch]
+                        y = signal[index, ch]
                 else:
                     y = np.zeros(len(ann_samp[ch]))
             except IndexError:
@@ -364,7 +492,7 @@ def plot_annotation(ann_samp, n_annot, ann_sym, signal, n_sig, fs, time_units,
                                           y[i]))
 
 
-def plot_ecg_grids(ecg_grids, fs, units, time_units, axes):
+def plot_ecg_grids(ecg_grids, fs, units, time_units, axes, sampling_freq=None):
     """
     Add ECG grids to the axes.
     
@@ -381,6 +509,10 @@ def plot_ecg_grids(ecg_grids, fs, units, time_units, axes):
         and 'hours'.
     axes : list
         The information needed for each subplot.
+    sampling_freq : number or sequence, optional
+        The sampling frequency or frequencies of the signals.  If this is a
+        list, it must have the same length as the number of channels.  If
+        unspecified, defaults to `fs`.
 
     Returns
     -------
@@ -390,15 +522,18 @@ def plot_ecg_grids(ecg_grids, fs, units, time_units, axes):
     if ecg_grids == 'all':
         ecg_grids = range(0, len(axes))
 
+    # Convert sampling_freq to a list if needed
+    sampling_freq = _get_sampling_freq(sampling_freq, len(axes), fs)
+
     for ch in ecg_grids:
         # Get the initial plot limits
         auto_xlims = axes[ch].get_xlim()
         auto_ylims= axes[ch].get_ylim()
 
         (major_ticks_x, minor_ticks_x, major_ticks_y,
             minor_ticks_y) = calc_ecg_grids(auto_ylims[0], auto_ylims[1],
-                                            units[ch], fs, auto_xlims[1],
-                                            time_units)
+                                            units[ch], sampling_freq[ch],
+                                            auto_xlims[1], time_units)
 
         min_x, max_x = np.min(minor_ticks_x), np.max(minor_ticks_x)
         min_y, max_y = np.min(minor_ticks_y), np.max(minor_ticks_y)
@@ -439,7 +574,7 @@ def calc_ecg_grids(minsig, maxsig, sig_units, fs, maxt, time_units):
     sig_units : list
         The units used for plotting each signal.
     fs : float
-        The sampling frequency of the record.
+        The sampling frequency of the signal.
     maxt : float
         The max time of the signal.
     time_units : str