This is a first working version of the converter incl. LFP generation

code/convert_to_nix.py

@@ -3,31 +3,124 @@ This script loads a complete session in the blackrock format and converts it to
 """
 import os
 
+import numpy as np
+import quantities as pq
+
 import neo
 from neo.test.tools import (assert_same_sub_schema,
                             assert_same_annotations,
                             assert_same_array_annotations)
+from elephant.signal_processing import butter
+
 from reachgraspio import reachgraspio
 
 
 # Choose which session you want to convert into a nix file
-# session = "l101210-001"
 session = "i140703-001"
+# session = "l101210-001"
 
-dataset_dir = '../datasets'
+# Input data. i.e., original Blackrock files and odML
+dataset_dir = '../datasets_blackrock'
 session_path = f'{dataset_dir}/{session}'
+odml_dir = os.path.join('..', 'datasets_blackrock')
+
+# Output for the nix files
+nix_dataset_dir = '../datasets_nix'
+nix_session_path = f'{nix_dataset_dir}/{session}'
+
 
 ##### LOAD BLACKROCK FILES ############
 session = reachgraspio.ReachGraspIO(
     filename=session_path,
-    odml_directory=dataset_dir,
+    odml_directory=odml_dir,
     verbose=False)
 
 block = session.read_block(lazy=False, load_waveforms=True)
 
+# =============================================================================
+# Create offline filtered LFP
+#
+# Here, we construct one offline filtered LFP from each ns5 (monkey L) or ns6
+# (monkey N) raw recording trace. For monkey N, this filtered LFP can be
+# compared to the LFPs in the ns2 file (note that monkey L contains only
+# behavioral signals in the ns2 file). Also, we assign telling names to each
+# Neo AnalogSignal, which is used for plotting later on in this script.
+# =============================================================================
+
+nsx_to_anasig_name = {2: 'LFP signal (online filtered)',
+                      5: 'raw signal',
+                      6: 'raw signal'}
+
+# this factor was experimentally determined as an approximate shift introduced
+# by the online filtering. Here, we integrate this shift such that the offline
+# filtered signal is aligned to the online filtered signal (if that were
+# available)
+time_shift_factor = -0.42*pq.ms
+
+filtered_anasig = None
+raw_anasig = None
+
+# identify neuronal signals and provide labels for plotting
+for anasig in block.segments[0].analogsignals:
+    # skip non-neuronal signals
+    if not anasig.annotations['neural_signal']:
+        continue
+
+    # identify nsx source of signals in this AnalogSignal object
+    if 'nsx' in anasig.annotations:
+        nsx = anasig.annotations['nsx']
+    else:
+        nsx = np.unique(anasig.array_annotations['nsx'])
+        assert len(nsx) == 1, 'Different nsx sources in AnalogSignal'
+        nsx = nsx[0]
+
+    if nsx == 2:
+        # AnalogSignal is LFP from ns2
+        filtered_anasig = anasig
+    elif nsx in [5, 6]:
+        # AnalogSignal is raw signal from ns5 or ns6
+        raw_anasig = anasig
+
+if filtered_anasig is None:
+    print("Filtering raw time series to obtain LFP")
+    for f in range(raw_anasig.shape[1]):
+        # filtering must be done channel by channel for memory reasons (requires approx. 32 GB RAM)
+        print(f"Processing channel {f}")
+
+        filtered_signal = butter(
+            raw_anasig[:, f],
+            highpass_freq=None,
+            lowpass_freq=250 * pq.Hz,
+            filter_function='sosfiltfilt',
+            order=4)
+
+        downsampled_signal=filtered_signal.downsample(30).time_shift(time_shift_factor)
+
+        # first run? Create a new Analogsignal
+        if f == 0:
+            offline_filtered_anasig = neo.AnalogSignal(
+                np.zeros((downsampled_signal.shape[0], raw_anasig.shape[1])) *\
+                    downsampled_signal.units,
+                t_start=downsampled_signal.t_start,
+                sampling_rate=downsampled_signal.sampling_rate)
+
+        offline_filtered_anasig[:, f] = downsampled_signal
+
+    if 'nsx' in anasig.annotations:
+        nsx = anasig.annotations['nsx']
+    else:
+        nsx = anasig.array_annotations["nsx"][0]
+
+    offline_filtered_anasig.name = f"NeuralTimeSeriesDownsampled"
+    offline_filtered_anasig.description = "Downsampled continuous neuronal recordings, where the downsampling was " \
+                           "performed off-line during post-processing"
+
+    # Attach all offline filtered LFPs to the segment of data
+    block.segments[0].analogsignals.append(offline_filtered_anasig)
+
 
 ##### SAVE NIX FILE ###################
-nix_filename = session_path + '.nix'
+nix_filename = nix_session_path + '.nix'
 if os.path.exists(nix_filename):
     print('Nix file already exists and will not be overwritten.')
 else:
@@ -45,8 +138,6 @@ with neo.NixIO(nix_filename, mode='ro') as io:
         for anasig_old, anasig_new in zip(seg_old.analogsignals, seg_new.analogsignals):
             # ignoring differences in the file_origin attribute
             anasig_old.file_origin = anasig_new.file_origin
-            # ignoring nix_name annotation generated by nix file
-            anasig_new.annotations.pop('nix_name')
 
             assert_same_sub_schema(anasig_old, anasig_new)
             assert_same_annotations(anasig_old, anasig_new)
@@ -57,8 +148,6 @@ with neo.NixIO(nix_filename, mode='ro') as io:
         for st_old, st_new in zip(seg_old.spiketrains, seg_new.spiketrains):
             # ignoring differences in the file_origin attribute
             st_old.file_origin = st_new.file_origin
-            # ignoring nix_name annotation generated by nix file
-            st_new.annotations.pop('nix_name')
 
             assert_same_sub_schema(st_old, st_new)
             assert_same_annotations(st_old, st_new)
@@ -69,8 +158,7 @@ with neo.NixIO(nix_filename, mode='ro') as io:
         for ev_old, ev_new in zip(seg_old.events, seg_new.events):
             # ignoring differences in the file_origin attribute
             ev_old.file_origin = ev_new.file_origin
-            # ignoring nix_name annotation generated by nix file
-            ev_new.annotations.pop('nix_name')
+
             # ignore list-array type changes
             if 'color_codes' in ev_old.annotations:
                 ev_old.annotations['color_codes'] = list(ev_old.annotations['color_codes'])
@@ -84,8 +172,6 @@ with neo.NixIO(nix_filename, mode='ro') as io:
         for ep_old, ep_new in zip(seg_old.epochs, seg_new.epochs):
             # ignoring differences in the file_origin attribute
             ep_old.file_origin = ep_new.file_origin
-            # ignoring nix_name annotation generated by nix file
-            ep_new.annotations.pop('nix_name')
 
             assert_same_sub_schema(ep_old, ep_new)
             assert_same_annotations(ep_old, ep_new)