Merge pull request #218 from MIT-LCP/read_edf_87

Produces record file from EDF format #87

Author: Lucas-Mc

sample-data/SC4001E0-PSG.edf

@@ -0,0 +1,8 @@
+SC4001E0_PSG 7 1 79500 16:13:00 24/04/1989
+SC4001E0_PSG.dat 16x100 10.6640625/uV 11 0 53 29442 0 EEG Fpz-Cz
+SC4001E0_PSG.dat 16x100 10.4198473282(5)/uV 11 0 -21 -8660 0 EEG Pz-Oz
+SC4001E0_PSG.dat 16x100 2.02923686819/uV 11 0 33 11375 0 EOG horizontal
+SC4001E0_PSG.dat 16 1(1) 11 0 -482 1209 0 Resp oro-nasal
+SC4001E0_PSG.dat 16 500(1)/uV 12 0 1776 19019 0 EMG submental
+SC4001E0_PSG.dat 16 930(-34468)/DegC 12 -58 133 24053 0 Temp rectal
+SC4001E0_PSG.dat 16 1 11 1 920 18040 0 Event marker

sample-data/SC4001E0_PSG.dat

@@ -0,0 +1,6 @@
+n16 5 100 3080000 22:34:47 01/01/2006
+n16.dat 16 44.9629231183/uV 15 0 -161 -9172 0 Fp2-F4
+n16.dat 16 44.9629231183/uV 15 0 1926 29657 0 F4-C4
+n16.dat 16 44.9629231183/uV 15 0 -4623 23837 0 C4-P4
+n16.dat 16 44.9629231183/uV 15 0 2472 25301 0 P4-O2
+n16.dat 16 44.9629231183/uV 15 0 -2354 -17711 0 C4-A1

sample-data/SC4001E0_PSG.hea

@@ -63,6 +63,7 @@
         'pandas>=0.19.1',
         'scipy>=0.19.0',
         'sklearn>=0.0',
+        'mne>=0.20.5'
     ],
 
     # List additional groups of dependencies here (e.g. development

sample-data/n16.dat

@@ -1,9 +1,9 @@
 import os
+import pdb
 import shutil
 import unittest
 
 import numpy as np
-
 import wfdb
 
 
@@ -249,6 +249,99 @@ def test_2e(self):
         sig_target = sig_target.reshape([977, 1])
         assert np.array_equal(sig, sig_target)
 
+    def test_2f(self):
+        """
+        EDF format conversion to MIT for uniform sample rates.
+
+        """
+        # Uniform sample rates
+        record_MIT = wfdb.rdrecord('sample-data/n16').__dict__
+        record_EDF = wfdb.rdrecord('sample-data/n16.edf').__dict__
+
+        fields = list(record_MIT.keys())
+        # Original MIT format method of checksum is outdated, sometimes
+        # the same value though
+        fields.remove('checksum')
+        # Original MIT format units are less comprehensive since they
+        # default to mV if unknown.. therefore added more default labels
+        fields.remove('units')
+
+        test_results = []
+        for field in fields:
+            # Signal value will be slightly off due to C to Python type conversion
+            if field == 'p_signal':
+                true_array = np.array(record_MIT[field])
+                pred_array = np.array(record_EDF[field])
+                sig_diff = np.abs((pred_array - true_array) / true_array)
+                sig_diff[sig_diff == -np.inf] = 0
+                sig_diff[sig_diff == np.inf] = 0
+                sig_diff = np.nanmean(sig_diff,0)
+                # 5% tolerance
+                if np.max(sig_diff) <= 5:
+                    test_results.append(True)
+                else:
+                    test_results.append(False)
+            elif field == 'init_value':
+                signal_diff = [abs(record_MIT[field][i] - record_EDF[field][i]) for i in range(len(record_MIT[field]))]
+                if abs(max(min(signal_diff), max(signal_diff), key=abs)) <= 2:
+                    test_results.append(True)
+                else:
+                    test_results.append(False)
+            else:
+                test_results.append(record_MIT[field] == record_MIT[field])
+
+        target_results = len(fields) * [True]
+        assert np.array_equal(test_results, target_results)
+
+    def test_2g(self):
+        """
+        EDF format conversion to MIT for non-uniform sample rates.
+
+        """
+        # Non-uniform sample rates
+        record_MIT = wfdb.rdrecord('sample-data/SC4001E0_PSG').__dict__
+        record_EDF = wfdb.rdrecord('sample-data/SC4001E0-PSG.edf').__dict__
+
+        fields = list(record_MIT.keys())
+        # Original MIT format method of checksum is outdated, sometimes
+        # the same value though
+        fields.remove('checksum')
+        # Original MIT format units are less comprehensive since they
+        # default to mV if unknown.. therefore added more default labels
+        fields.remove('units')
+        # Initial value of signal will be off due to resampling done by
+        # MNE in the EDF reading phase
+        fields.remove('init_value')
+        # Samples per frame will be off due to resampling done by MNE in
+        # the EDF reading phase... I should probably fix this later
+        fields.remove('samps_per_frame')
+
+        test_results = []
+        for field in fields:
+            # Signal value will be slightly off due to C to Python type conversion
+            if field == 'p_signal':
+                true_array = np.array(record_MIT[field])
+                pred_array = np.array(record_EDF[field])
+                sig_diff = np.abs((pred_array - true_array) / true_array)
+                sig_diff[sig_diff == -np.inf] = 0
+                sig_diff[sig_diff == np.inf] = 0
+                sig_diff = np.nanmean(sig_diff,0)
+                # 5% tolerance
+                if np.max(sig_diff) <= 5:
+                    test_results.append(True)
+                else:
+                    test_results.append(False)
+            elif field == 'init_value':
+                signal_diff = [abs(record_MIT[field][i] - record_EDF[field][i]) for i in range(len(record_MIT[field]))]
+                if abs(max(min(signal_diff), max(signal_diff), key=abs)) <= 2:
+                    test_results.append(True)
+                else:
+                    test_results.append(False)
+            else:
+                test_results.append(record_MIT[field] == record_MIT[field])
+
+        target_results = len(fields) * [True]
+        assert np.array_equal(test_results, target_results)
 
     # --------------------- 3. Multi-dat records --------------------- #
 
@@ -258,12 +351,12 @@ def test_3a(self):
         Target file created with:
             rdsamp -r sample-data/s0010_re | cut -f 2- > record-3a
         """
-        record= wfdb.rdrecord('sample-data/s0010_re', physical=False)
+        record = wfdb.rdrecord('sample-data/s0010_re', physical=False)
         sig = record.d_signal
         sig_target = np.genfromtxt('tests/target-output/record-3a')
 
         # Compare data streaming from Physionet
-        record_pn= wfdb.rdrecord('s0010_re', physical=False,
+        record_pn = wfdb.rdrecord('s0010_re', physical=False,
                                  pn_dir='ptbdb/patient001')
 
         # Test file writing
@@ -514,7 +607,7 @@ def test_multi_variable_b(self):
         from several segments.
 
         Target file created with:
-        rdsamp -r sample-data/multi-segment/s00001/s00001-2896-10-10-00-31 -f s14428364 -t s14428375 -P | cut -f 2- > record-multi-variable-b
+            rdsamp -r sample-data/multi-segment/s00001/s00001-2896-10-10-00-31 -f s14428364 -t s14428375 -P | cut -f 2- > record-multi-variable-b
         """
         record = wfdb.rdrecord('sample-data/multi-segment/s00001/s00001-2896-10-10-00-31',
                                sampfrom=14428364, sampto=14428375)

sample-data/n16.edf

@@ -1,5 +1,5 @@
 from .io.record import (Record, MultiRecord, rdheader, rdrecord, rdsamp,
-                        wrsamp, dl_database)
+                        wrsamp, dl_database, edf2mit)
 from .io.annotation import (Annotation, rdann, wrann, show_ann_labels,
                             show_ann_classes)
 from .io.download import get_dbs, get_record_list, dl_files, set_db_index_url

sample-data/n16.hea

@@ -726,6 +726,7 @@ def convert_dtype(self, physical, return_res, smooth_frames):
                         self.e_d_signal[ch] = self.e_d_signal[ch].astype(return_dtype, copy=False)
         return
 
+
     def calc_checksum(self, expanded=False):
         """
         Calculate the checksum(s) of the input signal.
@@ -859,7 +860,7 @@ def smooth_frames(self, sigtype='physical'):
 
 def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
                 samps_per_frame, skew, sampfrom, sampto, channels,
-                smooth_frames, ignore_skew, return_res=64):
+                smooth_frames, ignore_skew, no_file=False, sig_data=None, return_res=64):
     """
     Read the digital samples from a single segment record's associated
     dat file(s).
@@ -897,6 +898,12 @@ def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
         Specifies whether to apply the skew to align the signals in the
         output variable (False), or to ignore the skew field and load in
         all values contained in the dat files unaligned (True).
+    no_file : bool, optional
+        Used when using this function with just an array of signal data
+        and no associated file to read the data from.
+    sig_data : ndarray, optional
+        The signal data that would normally be imported using the associated
+        .dat and .hea files. Should only be used when no_file is set to True.
     return_res : int, optional
         The numpy array dtype of the returned signals. Options are: 64,
         32, 16, and 8, where the value represents the numpy int or float
@@ -918,6 +925,10 @@ def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
     specifications of the segment.
 
     """
+    # Check for valid inputs
+    if no_file and sig_data is None:
+        raise Exception('signal_dat empty: No signal data provided')
+
     # Avoid changing outer variables
     byte_offset = byte_offset[:]
     samps_per_frame = samps_per_frame[:]
@@ -984,10 +995,17 @@ def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
 
         # Read each wanted dat file and store signals
         for fn in w_file_name:
-            signals[:, out_dat_channel[fn]] = _rd_dat_signals(fn, dir_name, pn_dir,
-                w_fmt[fn], len(datchannel[fn]), sig_len, w_byte_offset[fn],
-                w_samps_per_frame[fn], w_skew[fn], sampfrom, sampto,
-                smooth_frames)[:, r_w_channel[fn]]
+            if no_file:
+                signals[:, out_dat_channel[fn]] = _rd_dat_signals(fn, dir_name,
+                    pn_dir, w_fmt[fn], len(datchannel[fn]), sig_len,
+                    w_byte_offset[fn], w_samps_per_frame[fn], w_skew[fn],
+                    sampfrom, sampto, smooth_frames, no_file=True,
+                    sig_data=sig_data)[:, r_w_channel[fn]]
+            else:
+                signals[:, out_dat_channel[fn]] = _rd_dat_signals(fn, dir_name,
+                    pn_dir, w_fmt[fn], len(datchannel[fn]), sig_len,
+                    w_byte_offset[fn], w_samps_per_frame[fn], w_skew[fn],
+                    sampfrom, sampto, smooth_frames)[:, r_w_channel[fn]]
 
     # Return each sample in signals with multiple samples/frame, without smoothing.
     # Return a list of numpy arrays for each signal.
@@ -996,10 +1014,16 @@ def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
 
         for fn in w_file_name:
             # Get the list of all signals contained in the dat file
-            datsignals = _rd_dat_signals(fn, dir_name, pn_dir, w_fmt[fn],
-                len(datchannel[fn]), sig_len, w_byte_offset[fn],
-                w_samps_per_frame[fn], w_skew[fn], sampfrom, sampto,
-                smooth_frames)
+            if no_file:
+                datsignals = _rd_dat_signals(fn, dir_name, pn_dir, w_fmt[fn],
+                    len(datchannel[fn]), sig_len, w_byte_offset[fn],
+                    w_samps_per_frame[fn], w_skew[fn], sampfrom, sampto,
+                    smooth_frames, no_file=True, sig_data=sig_data)
+            else:
+                datsignals = _rd_dat_signals(fn, dir_name, pn_dir, w_fmt[fn],
+                    len(datchannel[fn]), sig_len, w_byte_offset[fn],
+                    w_samps_per_frame[fn], w_skew[fn], sampfrom, sampto,
+                    smooth_frames)
 
             # Copy over the wanted signals
             for cn in range(len(out_dat_channel[fn])):
@@ -1010,7 +1034,7 @@ def _rd_segment(file_name, dir_name, pn_dir, fmt, n_sig, sig_len, byte_offset,
 
 def _rd_dat_signals(file_name, dir_name, pn_dir, fmt, n_sig, sig_len,
                    byte_offset, samps_per_frame, skew, sampfrom, sampto,
-                   smooth_frames):
+                   smooth_frames, no_file=False, sig_data=None):
     """
     Read all signals from a WFDB dat file.
 
@@ -1042,6 +1066,12 @@ def _rd_dat_signals(file_name, dir_name, pn_dir, fmt, n_sig, sig_len,
         The final sample number to be read from the signals.
     smooth_frames : bool
         Whether to smooth channels with multiple samples/frame.
+    no_file : bool, optional
+        Used when using this function with just an array of signal data
+        and no associated file to read the data from.
+    sig_data : ndarray, optional
+        The signal data that would normally be imported using the associated
+        .dat and .hea files. Should only be used when no_file is set to True.
 
     Returns
     -------
@@ -1058,6 +1088,10 @@ def _rd_dat_signals(file_name, dir_name, pn_dir, fmt, n_sig, sig_len,
     specifications of the segment.
 
     """
+    # Check for valid inputs
+    if no_file and sig_data is None:
+        raise Exception('signal_dat empty: No signal data provided')
+
     # Total number of samples per frame
     tsamps_per_frame = sum(samps_per_frame)
     # The signal length to read (per channel)
@@ -1086,26 +1120,27 @@ def _rd_dat_signals(file_name, dir_name, pn_dir, fmt, n_sig, sig_len,
     # already load samples.
 
     # Read values from dat file. Append bytes/samples if needed.
+    if no_file:
+        data_to_read = sig_data
+    else:
+        data_to_read = _rd_dat_file(file_name, dir_name, pn_dir, fmt,
+                                    start_byte, n_read_samples)
+
     if extra_flat_samples:
         if fmt in UNALIGNED_FMTS:
             # Extra number of bytes to append onto the bytes read from
             # the dat file.
             n_extra_bytes = total_process_bytes - total_read_bytes
 
-            sig_data = np.concatenate((_rd_dat_file(file_name, dir_name,
-                                                     pn_dir, fmt, start_byte,
-                                                     n_read_samples),
+            sig_data = np.concatenate((data_to_read,
                                         np.zeros(n_extra_bytes,
-                                                 dtype=np.dtype(DATA_LOAD_TYPES[fmt]))))
+                                                dtype=np.dtype(DATA_LOAD_TYPES[fmt]))))
         else:
-            sig_data = np.concatenate((_rd_dat_file(file_name, dir_name,
-                                                     pn_dir, fmt, start_byte,
-                                                     n_read_samples),
+            sig_data = np.concatenate((data_to_read,
                                         np.zeros(extra_flat_samples,
-                                                 dtype=np.dtype(DATA_LOAD_TYPES[fmt]))))
+                                                dtype=np.dtype(DATA_LOAD_TYPES[fmt]))))
     else:
-        sig_data = _rd_dat_file(file_name, dir_name, pn_dir, fmt, start_byte,
-                                 n_read_samples)
+        sig_data = data_to_read
 
     # Finish processing the read data into proper samples if not already