Merge pull request #1383 from manimoh/extract_events_from_SpikeGLX

Extract events from spike glx

Author: alejoe91

neo/rawio/spikeglxrawio.py

@@ -177,6 +177,22 @@ def _parse_header(self):
 
         # No events
         event_channels = []
+        # This is true only in case of 'nidq' stream
+        for stream_name in stream_names:
+            if "nidq" in stream_name:
+                info = self.signals_info_dict[0, stream_name]
+                if len(info["digital_channels"]) > 0:
+                    # add event channels
+                    for local_chan in info["digital_channels"]:
+                        chan_name = local_chan
+                        chan_id = f"{stream_name}#{chan_name}"
+                        event_channels.append((chan_name, chan_id, "event"))
+                    # add events_memmap
+                    data = np.memmap(info["bin_file"], dtype="int16", mode="r", offset=0, order="C")
+                    data = data.reshape(-1, info["num_chan"])
+                    # The digital word is usually the last channel, after all the individual analog channels
+                    extracted_word = data[:, len(info["analog_channels"])]
+                    self._events_memmap = np.unpackbits(extracted_word.astype(np.uint8)[:, None], axis=1)
         event_channels = np.array(event_channels, dtype=_event_channel_dtype)
 
         # No spikes
@@ -272,6 +288,44 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea
 
         return raw_signals
 
+    def _event_count(self, event_channel_idx, block_index=None, seg_index=None):
+        timestamps, _, _ = self._get_event_timestamps(block_index, seg_index, event_channel_idx, None, None)
+        return timestamps.size
+
+    def _get_event_timestamps(self, block_index, seg_index, event_channel_index, t_start=None, t_stop=None):
+        timestamps, durations, labels = [], None, []
+        info = self.signals_info_dict[0, "nidq"]  # There are no events that are not in the nidq stream
+        dig_ch = info["digital_channels"]
+        if len(dig_ch) > 0:
+            event_data = self._events_memmap
+            channel = dig_ch[event_channel_index]
+            ch_idx = 7 - int(channel[2:])  # They are in the reverse order
+            this_stream = event_data[:, ch_idx]
+            this_rising = np.where(np.diff(this_stream) == 1)[0] + 1
+            this_falling = (
+                np.where(np.diff(this_stream) == 255)[0] + 1
+            )  # because the data is in unsigned 8 bit, -1 = 255!
+            if len(this_rising) > 0:
+                timestamps.extend(this_rising)
+                labels.extend([f"{channel} ON"] * len(this_rising))
+            if len(this_falling) > 0:
+                timestamps.extend(this_falling)
+                labels.extend([f"{channel} OFF"] * len(this_falling))
+        timestamps = np.asarray(timestamps)
+        if len(labels) == 0:
+            labels = np.asarray(labels, dtype="U1")
+        else:
+            labels = np.asarray(labels)
+        return timestamps, durations, labels
+
+    def _rescale_event_timestamp(self, event_timestamps, dtype, event_channel_index):
+        info = self.signals_info_dict[0, "nidq"]  # There are no events that are not in the nidq stream
+        event_times = event_timestamps.astype(dtype) / float(info["sampling_rate"])
+        return event_times
+
+    def _rescale_epoch_duration(self, raw_duration, dtype, event_channel_index):
+        return None
+
 
 def scan_files(dirname):
     """
@@ -507,4 +561,14 @@ def extract_stream_info(meta_file, meta):
     info["channel_offsets"] = np.zeros(info["num_chan"])
     info["has_sync_trace"] = has_sync_trace
 
+    if "nidq" in device:
+        info["digital_channels"] = []
+        info["analog_channels"] = [channel for channel in info["channel_names"] if not channel.startswith("XD")]
+        # Digital/event channels are encoded within the digital word, so that will need more handling
+        for item in meta["niXDChans1"].split(","):
+            if ":" in item:
+                start, end = map(int, item.split(":"))
+                info["digital_channels"].extend([f"XD{i}" for i in range(start, end + 1)])
+            else:
+                info["digital_channels"].append(f"XD{int(item)}")
     return info

neo/test/rawiotest/test_spikeglxrawio.py

@@ -6,6 +6,7 @@
 
 from neo.rawio.spikeglxrawio import SpikeGLXRawIO
 from neo.test.rawiotest.common_rawio_test import BaseTestRawIO
+import numpy as np
 
 
 class TestSpikeGLXRawIO(BaseTestRawIO, unittest.TestCase):
@@ -87,6 +88,21 @@ def test_subset_with_sync(self):
         )
         assert chunk.shape[1] == 120
 
+    def test_nidq_digital_channel(self):
+        rawio_digital = SpikeGLXRawIO(self.get_local_path("spikeglx/DigitalChannelTest_g0"))
+        rawio_digital.parse_header()
+        # This data should have 8 event channels
+        assert np.shape(rawio_digital.header["event_channels"])[0] == 8
+
+        # Channel 0 in this data will have sync pulses at 1 Hz, let's confirm that
+        all_events = rawio_digital.get_event_timestamps(0, 0, 0)
+        on_events = np.where(all_events[2] == "XD0 ON")
+        on_ts = all_events[0][on_events]
+        on_ts_scaled = rawio_digital.rescale_event_timestamp(on_ts)
+        on_diff = np.diff(on_ts_scaled)
+        atol = 0.001
+        assert np.allclose(on_diff, 1, atol=atol)
+
 
 if __name__ == "__main__":
     unittest.main()