Merge pull request #331 from MIT-LCP/multisegment-multifrequency

Read multi-frequency multi-segment records

Author: tompollard

sample-data/multi-segment/041s/041s.hea

@@ -0,0 +1,3 @@
+041s/2 7 125 2000 8:26:04 26/10/1994
+041s01 1000
+041s02 1000

sample-data/multi-segment/041s/041s01.dat

@@ -0,0 +1,9 @@
+041s01 7 125 1000  8:26:04 26/10/1994
+041s01.dat 212x4 2000 12 0 168 -2716 0 III 
+041s01.dat 212x4 2000 12 0 2 -25019 0 I 
+041s01.dat 212x4 2000 12 0 155 -12467 0 V 
+041s01.dat 212 20(-1600)/mmHg 12 0 -242 -18875 0 ABP 
+041s01.dat 212 80(-1600)/mmHg 12 0 706 -5338 0 PAP 
+041s01.dat 212 2000 12 0 -841 30145 0 PLETH 
+041s01.dat 212 2000 12 0 401 3712 0 RESP 
+#Produced by xform from record mimicdb/041/04100001, beginning at s74000

sample-data/multi-segment/041s/041s01.hea

@@ -0,0 +1,9 @@
+041s02 7 125 1000  8:26:12 26/10/1994
+041s02.dat 212x4 2000 12 0 -103 -862 0 III 
+041s02.dat 212x4 2000 12 0 -64 14967 0 I 
+041s02.dat 212x4 2000 12 0 89 13162 0 V 
+041s02.dat 212 20(-1600)/mmHg 12 0 -715 -21117 0 ABP 
+041s02.dat 212 80(-1600)/mmHg 12 0 -583 -31770 0 PAP 
+041s02.dat 212 2000 12 0 -840 -31041 0 PLETH 
+041s02.dat 212 2000 12 0 -861 -31272 0 RESP 
+#Produced by xform from record mimicdb/041/04100002, beginning at 0:0

sample-data/multi-segment/041s/041s02.dat

@@ -746,6 +746,43 @@ def test_multi_fixed_c(self):
         np.testing.assert_equal(sig_round, sig_target)
         assert record.__eq__(record_named)
 
+    def test_multi_fixed_d(self):
+        """
+        Multi-segment, fixed layout, multi-frequency, selected channels
+
+        Target file created with:
+            rdsamp -r sample-data/multi-segment/041s/ -s 3 2 1 -H |
+            cut -f 2- | sed s/-32768/-2048/ |
+            gzip -9 -n > tests/target-output/record-multi-fixed-d.gz
+        """
+        record = wfdb.rdrecord(
+            "sample-data/multi-segment/041s/041s",
+            channels=[3, 2, 1],
+            physical=False,
+            smooth_frames=False,
+        )
+
+        # Convert expanded to uniform array (high-resolution)
+        sig = np.zeros((record.sig_len * 4, record.n_sig), dtype=int)
+        for i, s in enumerate(record.e_d_signal):
+            sig[:, i] = np.repeat(s, len(sig[:, i]) // len(s))
+
+        sig_target = np.genfromtxt(
+            "tests/target-output/record-multi-fixed-d.gz"
+        )
+
+        record_named = wfdb.rdrecord(
+            "sample-data/multi-segment/041s/041s",
+            channel_names=["ABP", "V", "I"],
+            physical=False,
+            smooth_frames=False,
+        )
+
+        # Sample values should match the output of rdsamp -H
+        np.testing.assert_array_equal(sig, sig_target)
+        # channel_names=[...] should give the same result as channels=[...]
+        self.assertEqual(record, record_named)
+
     def test_multi_variable_a(self):
         """
         Multi-segment, variable layout, selected duration, samples read
@@ -788,7 +825,7 @@ def test_multi_variable_b(self):
 
     def test_multi_variable_c(self):
         """
-        Multi-segment, variable layout, entire signal, physical
+        Multi-segment, variable layout, entire signal, physical, expanded
 
         The reference signal creation cannot be made with rdsamp
         directly because the WFDB c package (10.5.24) applies the single
@@ -811,9 +848,14 @@ def test_multi_variable_c(self):
 
         """
         record = wfdb.rdrecord(
-            "sample-data/multi-segment/s25047/s25047-2704-05-04-10-44"
+            "sample-data/multi-segment/s25047/s25047-2704-05-04-10-44",
+            smooth_frames=False,
         )
-        sig_round = np.round(record.p_signal, decimals=8)
+
+        # convert expanded to uniform array and round to 8 digits
+        sig_round = np.zeros((record.sig_len, record.n_sig))
+        for i in range(record.n_sig):
+            sig_round[:, i] = np.round(record.e_p_signal[i], decimals=8)
 
         sig_target_a = np.full((25740, 3), np.nan)
         sig_target_b = np.concatenate(

sample-data/multi-segment/041s/041s02.hea

@@ -295,7 +295,7 @@ def set_p_features(self, do_dac=False, expanded=False):
                 self.e_p_signal = self.dac(expanded)
 
             # Use e_p_signal to set fields
-            self.check_field("e_p_signal", channels="all")
+            self.check_field("e_p_signal", "all")
             self.sig_len = int(
                 len(self.e_p_signal[0]) / self.samps_per_frame[0]
             )
@@ -361,7 +361,7 @@ def set_d_features(self, do_adc=False, single_fmt=True, expanded=False):
         if expanded:
             # adc is performed.
             if do_adc:
-                self.check_field("e_p_signal", channels="all")
+                self.check_field("e_p_signal", "all")
 
                 # If there is no fmt set it, adc_gain, and baseline
                 if self.fmt is None:
@@ -393,7 +393,7 @@ def set_d_features(self, do_adc=False, single_fmt=True, expanded=False):
                 self.d_signal = self.adc(expanded)
 
             # Use e_d_signal to set fields
-            self.check_field("e_d_signal", channels="all")
+            self.check_field("e_d_signal", "all")
             self.sig_len = int(
                 len(self.e_d_signal[0]) / self.samps_per_frame[0]
             )

tests/target-output/record-multi-fixed-d.gz

@@ -510,13 +510,6 @@ def check_read_inputs(
                 "return_res must be one of the following when physical is True: 64, 32, 16"
             )
 
-        # Cannot expand multiple samples/frame for multi-segment records
-        if isinstance(self, MultiRecord):
-            if not smooth_frames:
-                raise ValueError(
-                    "This package version cannot expand all samples when reading multi-segment records. Must enable frame smoothing."
-                )
-
     def _adjust_datetime(self, sampfrom):
         """
         Adjust date and time fields to reflect user input if possible.
@@ -1269,7 +1262,7 @@ def _arrange_fields(
         self.n_seg = len(self.segments)
         self._adjust_datetime(sampfrom=sampfrom)
 
-    def multi_to_single(self, physical, return_res=64):
+    def multi_to_single(self, physical, return_res=64, expanded=False):
         """
         Create a Record object from the MultiRecord object. All signal
         segments will be combined into the new object's `p_signal` or
@@ -1283,6 +1276,11 @@ def multi_to_single(self, physical, return_res=64):
         return_res : int, optional
             The return resolution of the `p_signal` field. Options are:
             64, 32, and 16.
+        expanded : bool, optional
+            If false, combine the sample data from `p_signal` or `d_signal`
+            into a single two-dimensional array. If true, combine the
+            sample data from `e_p_signal` or `e_d_signal` into a list of
+            one-dimensional arrays.
 
         Returns
         -------
@@ -1300,7 +1298,14 @@ def multi_to_single(self, physical, return_res=64):
         # Figure out single segment fields to set for the new Record
         if self.layout == "fixed":
             # Get the fields from the first segment
-            for attr in ["fmt", "adc_gain", "baseline", "units", "sig_name"]:
+            for attr in [
+                "fmt",
+                "adc_gain",
+                "baseline",
+                "units",
+                "sig_name",
+                "samps_per_frame",
+            ]:
                 fields[attr] = getattr(self.segments[0], attr)
         else:
             # For variable layout records, inspect the segments for the
@@ -1311,9 +1316,14 @@ def multi_to_single(self, physical, return_res=64):
             # must have the same fmt, gain, baseline, and units for all
             # segments.
 
+            # For either physical or digital conversion, all signals
+            # of the same name must have the same samps_per_frame,
+            # which must match the value in the layout header.
+
             # The layout header should be updated at this point to
-            # reflect channels. We can depend on it for sig_name, but
-            # not for fmt, adc_gain, units, and baseline.
+            # reflect channels. We can depend on it for sig_name and
+            # samps_per_frame, but not for fmt, adc_gain, units, and
+            # baseline.
 
             # These signal names will be the key
             signal_names = self.segments[0].sig_name
@@ -1325,6 +1335,7 @@ def multi_to_single(self, physical, return_res=64):
                 "adc_gain": n_sig * [None],
                 "baseline": n_sig * [None],
                 "units": n_sig * [None],
+                "samps_per_frame": self.segments[0].samps_per_frame,
             }
 
             # For physical signals, mismatched fields will not be copied
@@ -1346,7 +1357,16 @@ def multi_to_single(self, physical, return_res=64):
                             reference_fields[field][ch] = item_ch
                         # mismatch case
                         elif reference_fields[field][ch] != item_ch:
-                            if physical:
+                            if field == "samps_per_frame":
+                                expected = reference_fields[field][ch]
+                                raise ValueError(
+                                    f"Incorrect samples per frame "
+                                    f"({item_ch} != {expected}) "
+                                    f"for signal {signal_names[ch]} "
+                                    f"in segment {seg.record_name} "
+                                    f"of {self.record_name}"
+                                )
+                            elif physical:
                                 mismatched_fields.append(field)
                             else:
                                 raise Exception(
@@ -1361,18 +1381,31 @@ def multi_to_single(self, physical, return_res=64):
 
         # Figure out signal attribute to set, and its dtype.
         if physical:
-            sig_attr = "p_signal"
+            if expanded:
+                sig_attr = "e_p_signal"
+            else:
+                sig_attr = "p_signal"
             # Figure out the largest required dtype
             dtype = _signal._np_dtype(return_res, discrete=False)
             nan_vals = np.array([self.n_sig * [np.nan]], dtype=dtype)
         else:
-            sig_attr = "d_signal"
+            if expanded:
+                sig_attr = "e_d_signal"
+            else:
+                sig_attr = "d_signal"
             # Figure out the largest required dtype
             dtype = _signal._np_dtype(return_res, discrete=True)
             nan_vals = np.array([_signal._digi_nan(fields["fmt"])], dtype=dtype)
 
+        samps_per_frame = fields["samps_per_frame"]
+
         # Initialize the full signal array
-        combined_signal = np.repeat(nan_vals, self.sig_len, axis=0)
+        if expanded:
+            combined_signal = []
+            for nan_val, spf in zip(nan_vals[0], samps_per_frame):
+                combined_signal.append(np.repeat(nan_val, spf * self.sig_len))
+        else:
+            combined_signal = np.repeat(nan_vals, self.sig_len, axis=0)
 
         # Start and end samples in the overall array to place the
         # segment samples into
@@ -1383,9 +1416,16 @@ def multi_to_single(self, physical, return_res=64):
             # Copy over the signals directly. Recall there are no
             # empty segments in fixed layout records.
             for i in range(self.n_seg):
-                combined_signal[start_samps[i] : end_samps[i], :] = getattr(
-                    self.segments[i], sig_attr
-                )
+                signals = getattr(self.segments[i], sig_attr)
+                if expanded:
+                    for ch in range(self.n_sig):
+                        start = start_samps[i] * samps_per_frame[ch]
+                        end = end_samps[i] * samps_per_frame[ch]
+                        combined_signal[ch][start:end] = signals[ch]
+                else:
+                    start = start_samps[i]
+                    end = end_samps[i]
+                    combined_signal[start:end, :] = signals
         else:
             # Copy over the signals into the matching channels
             for i in range(1, self.n_seg):
@@ -1396,12 +1436,20 @@ def multi_to_single(self, physical, return_res=64):
                     segment_channels = _get_wanted_channels(
                         fields["sig_name"], seg.sig_name, pad=True
                     )
+                    signals = getattr(seg, sig_attr)
                     for ch in range(self.n_sig):
                         # Copy over relevant signal
                         if segment_channels[ch] is not None:
-                            combined_signal[
-                                start_samps[i] : end_samps[i], ch
-                            ] = getattr(seg, sig_attr)[:, segment_channels[ch]]
+                            if expanded:
+                                signal = signals[segment_channels[ch]]
+                                start = start_samps[i] * samps_per_frame[ch]
+                                end = end_samps[i] * samps_per_frame[ch]
+                                combined_signal[ch][start:end] = signal
+                            else:
+                                signal = signals[:, segment_channels[ch]]
+                                start = start_samps[i]
+                                end = end_samps[i]
+                                combined_signal[start:end, ch] = signal
 
         # Create the single segment Record object and set attributes
         record = Record()
@@ -1411,9 +1459,9 @@ def multi_to_single(self, physical, return_res=64):
 
         # Use the signal to set record features
         if physical:
-            record.set_p_features()
+            record.set_p_features(expanded=expanded)
         else:
-            record.set_d_features()
+            record.set_d_features(expanded=expanded)
 
         return record
 
@@ -4168,6 +4216,7 @@ def rdrecord(
                     channels=seg_channels[i],
                     physical=physical,
                     pn_dir=pn_dir,
+                    smooth_frames=smooth_frames,
                     return_res=return_res,
                 )
 
@@ -4184,7 +4233,9 @@ def rdrecord(
         # Convert object into a single segment Record object
         if m2s:
             record = record.multi_to_single(
-                physical=physical, return_res=return_res
+                physical=physical,
+                expanded=(not smooth_frames),
+                return_res=return_res,
             )
 
     # Perform dtype conversion if necessary