blackrockrawio.py 83 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965
  1. """
  2. Module for reading data from files in the Blackrock in raw format.
  3. This work is based on:
  4. * Chris Rodgers - first version
  5. * Michael Denker, Lyuba Zehl - second version
  6. * Samuel Garcia - third version
  7. * Lyuba Zehl, Michael Denker - fourth version
  8. * Samuel Garcia, Julia Srenger - fifth version
  9. This IO supports reading only.
  10. This IO is able to read:
  11. * the nev file which contains spikes
  12. * ns1, ns2, .., ns6 files that contain signals at different sampling rates
  13. This IO can handle the following Blackrock file specifications:
  14. * 2.1
  15. * 2.2
  16. * 2.3
  17. The neural data channels are 1 - 128.
  18. The analog inputs are 129 - 144. (129 - 137 AC coupled, 138 - 144 DC coupled)
  19. spike- and event-data; 30000 Hz in NEV file.
  20. "ns1": "analog data: 500 Hz",
  21. "ns2": "analog data: 1000 Hz",
  22. "ns3": "analog data: 2000 Hz",
  23. "ns4": "analog data: 10000 Hz",
  24. "ns5": "analog data: 30000 Hz",
  25. "ns6": "analog data: 30000 Hz (no digital filter)"
  26. The possible file extensions of the Cerebus system and their content:
  27. ns1: contains analog data; sampled at 500 Hz (+ digital filters)
  28. ns2: contains analog data; sampled at 1000 Hz (+ digital filters)
  29. ns3: contains analog data; sampled at 2000 Hz (+ digital filters)
  30. ns4: contains analog data; sampled at 10000 Hz (+ digital filters)
  31. ns5: contains analog data; sampled at 30000 Hz (+ digital filters)
  32. ns6: contains analog data; sampled at 30000 Hz (no digital filters)
  33. nev: contains spike- and event-data; sampled at 30000 Hz
  34. sif: contains institution and patient info (XML)
  35. ccf: contains Cerebus configurations
  36. TODO:
  37. * videosync events (file spec 2.3)
  38. * tracking events (file spec 2.3)
  39. * buttontrigger events (file spec 2.3)
  40. * config events (file spec 2.3)
  41. * check left sweep settings of Blackrock
  42. * check nsx offsets (file spec 2.1)
  43. * add info of nev ext header (NSASEXEX) to non-neural events
  44. (file spec 2.1 and 2.2)
  45. * read sif file information
  46. * read ccf file information
  47. * fix reading of periodic sampling events (non-neural event type)
  48. (file spec 2.1 and 2.2)
  49. """
  50. import datetime
  51. import os
  52. import re
  53. import warnings
  54. import numpy as np
  55. import quantities as pq
  56. from .baserawio import (BaseRawIO, _signal_channel_dtype, _unit_channel_dtype,
  57. _event_channel_dtype)
  58. class BlackrockRawIO(BaseRawIO):
  59. """
  60. Class for reading data in from a file set recorded by the Blackrock
  61. (Cerebus) recording system.
  62. Upon initialization, the class is linked to the available set of Blackrock
  63. files.
  64. Note: This routine will handle files according to specification 2.1, 2.2,
  65. and 2.3. Recording pauses that may occur in file specifications 2.2 and
  66. 2.3 are automatically extracted and the data set is split into different
  67. segments.
  68. The Blackrock data format consists not of a single file, but a set of
  69. different files. This constructor associates itself with a set of files
  70. that constitute a common data set. By default, all files belonging to
  71. the file set have the same base name, but different extensions.
  72. However, by using the override parameters, individual filenames can
  73. be set.
  74. Args:
  75. filename (string):
  76. File name (without extension) of the set of Blackrock files to
  77. associate with. Any .nsX or .nev, .sif, or .ccf extensions are
  78. ignored when parsing this parameter.
  79. nsx_override (string):
  80. File name of the .nsX files (without extension). If None,
  81. filename is used.
  82. Default: None.
  83. nev_override (string):
  84. File name of the .nev file (without extension). If None,
  85. filename is used.
  86. Default: None.
  87. nsx_to_load (int, list, 'max', 'all' (=None)) default None:
  88. IDs of nsX file from which to load data, e.g., if set to
  89. 5 only data from the ns5 file are loaded.
  90. If 'all', then all nsX will be loaded.
  91. Contrary to previsous version of the IO (<0.7), nsx_to_load
  92. must be set at the init before parse_header().
  93. Examples:
  94. >>> reader = BlackrockRawIO(filename='FileSpec2.3001', nsx_to_load=5)
  95. >>> reader.parse_header()
  96. Inspect a set of file consisting of files FileSpec2.3001.ns5 and
  97. FileSpec2.3001.nev
  98. >>> print(reader)
  99. Display all informations about signal channels, units, segment size....
  100. """
  101. extensions = ['ns' + str(_) for _ in range(1, 7)]
  102. extensions.extend(['nev', ]) # 'sif', 'ccf' not yet supported
  103. rawmode = 'multi-file'
  104. def __init__(self, filename=None, nsx_override=None, nev_override=None,
  105. nsx_to_load=None, verbose=False):
  106. """
  107. Initialize the BlackrockIO class.
  108. """
  109. BaseRawIO.__init__(self)
  110. self.filename = filename
  111. # remove extension from base _filenames
  112. for ext in self.extensions:
  113. if self.filename.endswith(os.path.extsep + ext):
  114. self.filename = self.filename.replace(os.path.extsep + ext, '')
  115. self.nsx_to_load = nsx_to_load
  116. # remove extensions from overrides
  117. self._filenames = {}
  118. if nsx_override:
  119. self._filenames['nsx'] = re.sub(
  120. os.path.extsep + 'ns[1,2,3,4,5,6]$', '', nsx_override)
  121. else:
  122. self._filenames['nsx'] = self.filename
  123. if nev_override:
  124. self._filenames['nev'] = re.sub(
  125. os.path.extsep + 'nev$', '', nev_override)
  126. else:
  127. self._filenames['nev'] = self.filename
  128. # check which files are available
  129. self._avail_files = dict.fromkeys(self.extensions, False)
  130. self._avail_nsx = []
  131. for ext in self.extensions:
  132. if ext.startswith('ns'):
  133. file2check = ''.join(
  134. [self._filenames['nsx'], os.path.extsep, ext])
  135. else:
  136. file2check = ''.join(
  137. [self._filenames[ext], os.path.extsep, ext])
  138. if os.path.exists(file2check):
  139. self._avail_files[ext] = True
  140. if ext.startswith('ns'):
  141. self._avail_nsx.append(int(ext[-1]))
  142. if not self._avail_files['nev'] and not self._avail_nsx:
  143. raise IOError("No Blackrock files found in specified path")
  144. # These dictionaries are used internally to map the file specification
  145. # revision of the nsx and nev files to one of the reading routines
  146. # NSX
  147. self.__nsx_header_reader = {
  148. '2.1': self.__read_nsx_header_variant_a,
  149. '2.2': self.__read_nsx_header_variant_b,
  150. '2.3': self.__read_nsx_header_variant_b}
  151. self.__nsx_dataheader_reader = {
  152. '2.1': self.__read_nsx_dataheader_variant_a,
  153. '2.2': self.__read_nsx_dataheader_variant_b,
  154. '2.3': self.__read_nsx_dataheader_variant_b}
  155. self.__nsx_data_reader = {
  156. '2.1': self.__read_nsx_data_variant_a,
  157. '2.2': self.__read_nsx_data_variant_b,
  158. '2.3': self.__read_nsx_data_variant_b}
  159. self.__nsx_params = {
  160. '2.1': self.__get_nsx_param_variant_a,
  161. '2.2': self.__get_nsx_param_variant_b,
  162. '2.3': self.__get_nsx_param_variant_b}
  163. # NEV
  164. self.__nev_header_reader = {
  165. '2.1': self.__read_nev_header_variant_a,
  166. '2.2': self.__read_nev_header_variant_b,
  167. '2.3': self.__read_nev_header_variant_c}
  168. self.__nev_data_reader = {
  169. '2.1': self.__read_nev_data_variant_a,
  170. '2.2': self.__read_nev_data_variant_a,
  171. '2.3': self.__read_nev_data_variant_b}
  172. self.__waveform_size = {
  173. '2.1': self.__get_waveform_size_variant_a,
  174. '2.2': self.__get_waveform_size_variant_a,
  175. '2.3': self.__get_waveform_size_variant_b}
  176. self.__channel_labels = {
  177. '2.1': self.__get_channel_labels_variant_a,
  178. '2.2': self.__get_channel_labels_variant_b,
  179. '2.3': self.__get_channel_labels_variant_b}
  180. self.__nonneural_evdicts = {
  181. '2.1': self.__get_nonneural_evdicts_variant_a,
  182. '2.2': self.__get_nonneural_evdicts_variant_a,
  183. '2.3': self.__get_nonneural_evdicts_variant_b}
  184. self.__comment_evdict = {
  185. '2.1': self.__get_comment_evdict_variant_a,
  186. '2.2': self.__get_comment_evdict_variant_a,
  187. '2.3': self.__get_comment_evdict_variant_a}
  188. def _parse_header(self):
  189. main_sampling_rate = 30000.
  190. event_channels = []
  191. unit_channels = []
  192. sig_channels = []
  193. # Step1 NEV file
  194. if self._avail_files['nev']:
  195. # Load file spec and headers of available
  196. # read nev file specification
  197. self.__nev_spec = self.__extract_nev_file_spec()
  198. # read nev headers
  199. self.__nev_basic_header, self.__nev_ext_header = \
  200. self.__nev_header_reader[self.__nev_spec]()
  201. self.nev_data = self.__nev_data_reader[self.__nev_spec]()
  202. spikes, spike_segment_ids = self.nev_data['Spikes']
  203. # scan all channel to get number of Unit
  204. unit_channels = []
  205. self.internal_unit_ids = [] # pair of chan['packet_id'], spikes['unit_class_nb']
  206. for i in range(len(self.__nev_ext_header[b'NEUEVWAV'])):
  207. channel_id = self.__nev_ext_header[b'NEUEVWAV']['electrode_id'][i]
  208. chan_mask = (spikes['packet_id'] == channel_id)
  209. chan_spikes = spikes[chan_mask]
  210. all_unit_id = np.unique(chan_spikes['unit_class_nb'])
  211. for u, unit_id in enumerate(all_unit_id):
  212. self.internal_unit_ids.append((channel_id, unit_id))
  213. name = "ch{}#{}".format(channel_id, unit_id)
  214. _id = "Unit {}".format(1000 * channel_id + unit_id)
  215. wf_gain = self.__nev_params('digitization_factor')[channel_id] / 1000.
  216. wf_offset = 0.
  217. wf_units = 'uV'
  218. # TODO: Double check if this is the correct assumption (10 samples)
  219. # default value: threshold crossing after 10 samples of waveform
  220. wf_left_sweep = 10
  221. wf_sampling_rate = main_sampling_rate
  222. unit_channels.append((name, _id, wf_units, wf_gain,
  223. wf_offset, wf_left_sweep, wf_sampling_rate))
  224. # scan events
  225. # NonNeural: serial and digital input
  226. events_data, event_segment_ids = self.nev_data['NonNeural']
  227. ev_dict = self.__nonneural_evdicts[self.__nev_spec](events_data)
  228. if 'Comments' in self.nev_data:
  229. comments_data, comments_segment_ids = self.nev_data['Comments']
  230. ev_dict.update(self.__comment_evdict[self.__nev_spec](comments_data))
  231. for ev_name in ev_dict:
  232. event_channels.append((ev_name, '', 'event'))
  233. # TODO: TrackingEvents
  234. # TODO: ButtonTrigger
  235. # TODO: VideoSync
  236. # Step2 NSX file
  237. # Load file spec and headers of available nsx files
  238. self.__nsx_spec = {}
  239. self.__nsx_basic_header = {}
  240. self.__nsx_ext_header = {}
  241. self.__nsx_data_header = {}
  242. for nsx_nb in self._avail_nsx:
  243. spec = self.__nsx_spec[nsx_nb] = self.__extract_nsx_file_spec(nsx_nb)
  244. # read nsx headers
  245. self.__nsx_basic_header[nsx_nb], self.__nsx_ext_header[nsx_nb] = \
  246. self.__nsx_header_reader[spec](nsx_nb)
  247. # Read nsx data header(s)
  248. # for nsxdef get_analogsignal_shape(self, block_index, seg_index):
  249. self.__nsx_data_header[nsx_nb] = self.__nsx_dataheader_reader[spec](nsx_nb)
  250. # nsx_to_load can be either int, list, 'max', all' (aka None)
  251. # here make a list only
  252. if self.nsx_to_load is None or self.nsx_to_load == 'all':
  253. self.nsx_to_load = list(self._avail_nsx)
  254. elif self.nsx_to_load == 'max':
  255. if len(self._avail_nsx):
  256. self.nsx_to_load = [max(self._avail_nsx)]
  257. else:
  258. self.nsx_to_load = []
  259. elif isinstance(self.nsx_to_load, int):
  260. self.nsx_to_load = [self.nsx_to_load]
  261. elif isinstance(self.nsx_to_load, list):
  262. pass
  263. else:
  264. raise(ValueError('nsx_to_load is wrong'))
  265. assert all(nsx_nb in self._avail_nsx for nsx_nb in self.nsx_to_load),\
  266. 'nsx_to_load do not match available nsx list'
  267. # check that all files come from the same specification
  268. all_spec = [self.__nsx_spec[nsx_nb] for nsx in self.nsx_to_load]
  269. if self._avail_files['nev']:
  270. all_spec.append(self.__nev_spec)
  271. assert all(all_spec[0] == spec for spec in all_spec), \
  272. "Files don't have the same internal version"
  273. if len(self.nsx_to_load) > 0 and \
  274. self.__nsx_spec[self.nsx_to_load[0]] == '2.1' and \
  275. not self._avail_files['nev']:
  276. pass
  277. # Because rescaling to volts requires information from nev file (dig_factor)
  278. # Remove if raw loading becomes possible
  279. # raise IOError("For loading Blackrock file version 2.1 .nev files are required!")
  280. # This requires nsX to be parsed already
  281. # Needs to be called when no nsX are available as well in order to warn the user
  282. if self._avail_files['nev']:
  283. for nsx_nb in self.nsx_to_load:
  284. self.__match_nsx_and_nev_segment_ids(nsx_nb)
  285. # usefull to get local channel index in nsX from the global channel index
  286. local_sig_indexes = []
  287. self.nsx_datas = {}
  288. self.sig_sampling_rates = {}
  289. if len(self.nsx_to_load) > 0:
  290. for nsx_nb in self.nsx_to_load:
  291. spec = self.__nsx_spec[nsx_nb]
  292. self.nsx_datas[nsx_nb] = self.__nsx_data_reader[spec](nsx_nb)
  293. sr = float(main_sampling_rate / self.__nsx_basic_header[nsx_nb]['period'])
  294. self.sig_sampling_rates[nsx_nb] = sr
  295. if spec in ['2.2', '2.3']:
  296. ext_header = self.__nsx_ext_header[nsx_nb]
  297. elif spec == '2.1':
  298. ext_header = []
  299. keys = ['labels', 'units', 'min_analog_val',
  300. 'max_analog_val', 'min_digital_val', 'max_digital_val']
  301. params = self.__nsx_params[spec](nsx_nb)
  302. for i in range(len(params['labels'])):
  303. d = {}
  304. for key in keys:
  305. d[key] = params[key][i]
  306. ext_header.append(d)
  307. for i, chan in enumerate(ext_header):
  308. if spec in ['2.2', '2.3']:
  309. ch_name = chan['electrode_label'].decode()
  310. ch_id = chan['electrode_id']
  311. units = chan['units'].decode()
  312. elif spec == '2.1':
  313. ch_name = chan['labels']
  314. ch_id = self.__nsx_ext_header[nsx_nb][i]['electrode_id']
  315. units = chan['units']
  316. sig_dtype = 'int16'
  317. # max_analog_val/min_analog_val/max_digital_val/min_analog_val are int16!!!!!
  318. # dangarous situation so cast to float everyone
  319. if np.isnan(float(chan['min_analog_val'])):
  320. gain = 1
  321. offset = 0
  322. else:
  323. gain = (float(chan['max_analog_val']) - float(chan['min_analog_val'])) / \
  324. (float(chan['max_digital_val']) - float(chan['min_digital_val']))
  325. offset = -float(chan['min_digital_val']) \
  326. * gain + float(chan['min_analog_val'])
  327. group_id = nsx_nb
  328. sig_channels.append((ch_name, ch_id, sr, sig_dtype,
  329. units, gain, offset, group_id,))
  330. local_sig_indexes.extend(range(len(ext_header)))
  331. self._local_sig_indexes = np.array(local_sig_indexes)
  332. # check nb segment per nsx
  333. nb_segments_for_nsx = [len(self.nsx_datas[nsx_nb]) for nsx_nb in self.nsx_to_load]
  334. assert all(nb == nb_segments_for_nsx[0] for nb in nb_segments_for_nsx),\
  335. 'Segment nb not consistanent across nsX files'
  336. self._nb_segment = nb_segments_for_nsx[0]
  337. self.__delete_empty_segments()
  338. # t_start/t_stop for segment are given by nsx limits or nev limits
  339. self._sigs_t_starts = {nsx_nb: [] for nsx_nb in self.nsx_to_load}
  340. self._seg_t_starts, self._seg_t_stops = [], []
  341. for data_bl in range(self._nb_segment):
  342. t_stop = 0.
  343. for nsx_nb in self.nsx_to_load:
  344. length = self.nsx_datas[nsx_nb][data_bl].shape[0]
  345. if self.__nsx_data_header[nsx_nb] is None:
  346. t_start = 0.
  347. else:
  348. t_start = self.__nsx_data_header[nsx_nb][data_bl]['timestamp'] / \
  349. self.__nsx_basic_header[nsx_nb]['timestamp_resolution']
  350. t_stop = max(t_stop, t_start + length / self.sig_sampling_rates[nsx_nb])
  351. self._sigs_t_starts[nsx_nb].append(t_start)
  352. if self._avail_files['nev']:
  353. max_nev_time = 0
  354. for k, (data, ev_ids) in self.nev_data.items():
  355. segment_mask = ev_ids == data_bl
  356. if data[segment_mask].size > 0:
  357. t = data[segment_mask][-1]['timestamp'] / self.__nev_basic_header[
  358. 'timestamp_resolution']
  359. max_nev_time = max(max_nev_time, t)
  360. if max_nev_time > t_stop:
  361. t_stop = max_nev_time
  362. min_nev_time = max_nev_time
  363. for k, (data, ev_ids) in self.nev_data.items():
  364. segment_mask = ev_ids == data_bl
  365. if data[segment_mask].size > 0:
  366. t = data[segment_mask][0]['timestamp'] / self.__nev_basic_header[
  367. 'timestamp_resolution']
  368. min_nev_time = min(min_nev_time, t)
  369. if min_nev_time < t_start:
  370. t_start = min_nev_time
  371. self._seg_t_starts.append(t_start)
  372. self._seg_t_stops.append(float(t_stop))
  373. else:
  374. # When only nev is available, only segments that are documented in nev can be detected
  375. max_nev_times = {}
  376. min_nev_times = {}
  377. # Find maximal and minimal time for each nev segment
  378. for k, (data, ev_ids) in self.nev_data.items():
  379. for i in np.unique(ev_ids):
  380. mask = [ev_ids == i]
  381. curr_data = data[mask]
  382. if curr_data.size > 0:
  383. if max(curr_data['timestamp']) >= max_nev_times.get(i, 0):
  384. max_nev_times[i] = max(curr_data['timestamp'])
  385. if min(curr_data['timestamp']) <= min_nev_times.get(i,
  386. max_nev_times[i]):
  387. min_nev_times[i] = min(curr_data['timestamp'])
  388. # Calculate t_start and t_stop for each segment in seconds
  389. resolution = self.__nev_basic_header['timestamp_resolution']
  390. self._seg_t_starts = [v / float(resolution) for k, v in sorted(min_nev_times.items())]
  391. self._seg_t_stops = [v / float(resolution) for k, v in sorted(max_nev_times.items())]
  392. self._nb_segment = len(self._seg_t_starts)
  393. self._sigs_t_starts = [None] * self._nb_segment
  394. # finalize header
  395. unit_channels = np.array(unit_channels, dtype=_unit_channel_dtype)
  396. event_channels = np.array(event_channels, dtype=_event_channel_dtype)
  397. sig_channels = np.array(sig_channels, dtype=_signal_channel_dtype)
  398. self.header = {}
  399. self.header['nb_block'] = 1
  400. self.header['nb_segment'] = [self._nb_segment]
  401. self.header['signal_channels'] = sig_channels
  402. self.header['unit_channels'] = unit_channels
  403. self.header['event_channels'] = event_channels
  404. rec_datetime = self.__nev_params('rec_datetime') if self._avail_files['nev'] else None
  405. # Put annotations at some places for compatibility
  406. # with previous BlackrockIO version
  407. self._generate_minimal_annotations()
  408. block_ann = self.raw_annotations['blocks'][0]
  409. block_ann['description'] = 'Block of data from Blackrock file set.'
  410. block_ann['file_origin'] = self.filename
  411. block_ann['name'] = "Blackrock Data Block"
  412. block_ann['rec_datetime'] = rec_datetime
  413. block_ann['avail_file_set'] = [k for k, v in self._avail_files.items() if v]
  414. block_ann['avail_nsx'] = self._avail_nsx
  415. block_ann['avail_nev'] = self._avail_files['nev']
  416. # 'sif' and 'ccf' files not yet supported
  417. # block_ann['avail_sif'] = self._avail_files['sif']
  418. # block_ann['avail_ccf'] = self._avail_files['ccf']
  419. block_ann['rec_pauses'] = False
  420. for c in range(unit_channels.size):
  421. unit_ann = self.raw_annotations['unit_channels'][c]
  422. channel_id, unit_id = self.internal_unit_ids[c]
  423. unit_ann['channel_id'] = self.internal_unit_ids[c][0]
  424. unit_ann['unit_id'] = self.internal_unit_ids[c][1]
  425. unit_ann['unit_tag'] = {0: 'unclassified', 255: 'noise'}.get(unit_id, str(unit_id))
  426. unit_ann['description'] = 'Unit channel_id: {}, unit_id: {}, unit_tag: {}'.format(
  427. channel_id, unit_id, unit_ann['unit_tag'])
  428. flt_type = {0: 'None', 1: 'Butterworth'}
  429. for c in range(sig_channels.size):
  430. chidx_ann = self.raw_annotations['signal_channels'][c]
  431. if self._avail_files['nev']:
  432. neuevwav = self.__nev_ext_header[b'NEUEVWAV']
  433. if sig_channels[c]['id'] in neuevwav['electrode_id']:
  434. get_idx = list(neuevwav['electrode_id']).index(sig_channels[c]['id'])
  435. chidx_ann['connector_ID'] = neuevwav['physical_connector'][get_idx]
  436. chidx_ann['connector_pinID'] = neuevwav['connector_pin'][get_idx]
  437. chidx_ann['nev_dig_factor'] = neuevwav['digitization_factor'][get_idx]
  438. chidx_ann['nev_energy_threshold'] = neuevwav['energy_threshold'][
  439. get_idx] * chidx_ann['nev_dig_factor'] / 1000 * pq.uV
  440. chidx_ann['nev_hi_threshold'] = neuevwav['hi_threshold'][get_idx] * chidx_ann[
  441. 'nev_dig_factor'] / 1000 * pq.uV
  442. chidx_ann['nev_lo_threshold'] = neuevwav['lo_threshold'][get_idx] * chidx_ann[
  443. 'nev_dig_factor'] / 1000 * pq.uV
  444. chidx_ann['nb_sorted_units'] = neuevwav['nb_sorted_units'][get_idx]
  445. chidx_ann['waveform_size'] = self.__waveform_size[self.__nev_spec](
  446. )[sig_channels[c]['id']] * self.__nev_params('waveform_time_unit')
  447. if self.__nev_spec in ['2.2', '2.3']:
  448. neuevflt = self.__nev_ext_header[b'NEUEVFLT']
  449. get_idx = list(
  450. neuevflt['electrode_id']).index(
  451. sig_channels[c]['id'])
  452. # filter type codes (extracted from blackrock manual)
  453. chidx_ann['nev_hi_freq_corner'] = neuevflt['hi_freq_corner'][
  454. get_idx] / 1000. * pq.Hz
  455. chidx_ann['nev_hi_freq_order'] = neuevflt['hi_freq_order'][get_idx]
  456. chidx_ann['nev_hi_freq_type'] = flt_type[neuevflt['hi_freq_type'][
  457. get_idx]]
  458. chidx_ann['nev_lo_freq_corner'] = neuevflt['lo_freq_corner'][
  459. get_idx] / 1000. * pq.Hz
  460. chidx_ann['nev_lo_freq_order'] = neuevflt['lo_freq_order'][get_idx]
  461. chidx_ann['nev_lo_freq_type'] = flt_type[neuevflt['lo_freq_type'][
  462. get_idx]]
  463. if self.__nsx_spec[self.nsx_to_load[0]] in ['2.2', '2.3'] and self.__nsx_ext_header:
  464. # It does not matter which nsX file to ask for this info
  465. k = list(self.__nsx_ext_header.keys())[0]
  466. if sig_channels[c]['id'] in self.__nsx_ext_header[k]['electrode_id']:
  467. get_idx = list(
  468. self.__nsx_ext_header[k]['electrode_id']).index(
  469. sig_channels[c]['id'])
  470. chidx_ann['connector_ID'] = self.__nsx_ext_header[k]['physical_connector'][
  471. get_idx]
  472. chidx_ann['connector_pinID'] = self.__nsx_ext_header[k]['connector_pin'][
  473. get_idx]
  474. chidx_ann['nsx_hi_freq_corner'] = \
  475. self.__nsx_ext_header[k]['hi_freq_corner'][get_idx] / 1000. * pq.Hz
  476. chidx_ann['nsx_lo_freq_corner'] = \
  477. self.__nsx_ext_header[k]['lo_freq_corner'][get_idx] / 1000. * pq.Hz
  478. chidx_ann['nsx_hi_freq_order'] = self.__nsx_ext_header[k][
  479. 'hi_freq_order'][get_idx]
  480. chidx_ann['nsx_lo_freq_order'] = self.__nsx_ext_header[k][
  481. 'lo_freq_order'][get_idx]
  482. chidx_ann['nsx_hi_freq_type'] = flt_type[
  483. self.__nsx_ext_header[k]['hi_freq_type'][get_idx]]
  484. chidx_ann['nsx_lo_freq_type'] = flt_type[
  485. self.__nsx_ext_header[k]['hi_freq_type'][get_idx]]
  486. for seg_index in range(self._nb_segment):
  487. seg_ann = block_ann['segments'][seg_index]
  488. seg_ann['file_origin'] = self.filename
  489. seg_ann['name'] = "Segment {}".format(seg_index)
  490. seg_ann['description'] = "Segment containing data from t_start to t_stop"
  491. if seg_index == 0:
  492. # if more than 1 segment means pause
  493. # so datetime is valide only for seg_index=0
  494. seg_ann['rec_datetime'] = rec_datetime
  495. for c in range(sig_channels.size):
  496. nsx_nb = sig_channels['group_id'][c]
  497. anasig_an = seg_ann['signals'][c]
  498. desc = "AnalogSignal {} from channel_id: {}, label: {}, nsx: {}".format(
  499. c, sig_channels['id'][c], sig_channels['name'][c], nsx_nb)
  500. anasig_an['description'] = desc
  501. anasig_an['file_origin'] = self._filenames['nsx'] + '.ns' + str(nsx_nb)
  502. anasig_an['nsx'] = nsx_nb
  503. chidx_ann = self.raw_annotations['signal_channels'][c]
  504. chidx_ann['description'] = 'Container for Units and AnalogSignals of ' \
  505. 'one recording channel across segments.'
  506. for c in range(unit_channels.size):
  507. channel_id, unit_id = self.internal_unit_ids[c]
  508. st_ann = seg_ann['units'][c]
  509. unit_ann = self.raw_annotations['unit_channels'][c]
  510. st_ann.update(unit_ann)
  511. st_ann['description'] = 'SpikeTrain channel_id: {}, unit_id: {}'.format(
  512. channel_id, unit_id)
  513. st_ann['file_origin'] = self._filenames['nev'] + '.nev'
  514. if self._avail_files['nev']:
  515. ev_dict = self.__nonneural_evdicts[self.__nev_spec](events_data)
  516. if 'Comments' in self.nev_data:
  517. ev_dict.update(self.__comment_evdict[self.__nev_spec](comments_data))
  518. color_codes = ["#{:08X}".format(code) for code in comments_data['color']]
  519. color_codes = np.array(color_codes, dtype='S9')
  520. for c in range(event_channels.size):
  521. # Next line makes ev_ann a reference to seg_ann['events'][c]
  522. ev_ann = seg_ann['events'][c]
  523. name = event_channels['name'][c]
  524. ev_ann['description'] = ev_dict[name]['desc']
  525. ev_ann['file_origin'] = self._filenames['nev'] + '.nev'
  526. if name == 'comments':
  527. ev_ann['color_codes'] = color_codes
  528. def _source_name(self):
  529. return self.filename
  530. def _segment_t_start(self, block_index, seg_index):
  531. return self._seg_t_starts[seg_index]
  532. def _segment_t_stop(self, block_index, seg_index):
  533. return self._seg_t_stops[seg_index]
  534. def _get_nsx_and_local_indexes(self, channel_indexes):
  535. # internal helper to get nsx number and local channel index
  536. # from global channel indexes
  537. # when this is called channell_indexes are alwas in the same group_id
  538. # this is checked at BaseRaw level
  539. if channel_indexes is None:
  540. channel_indexes = slice(None)
  541. nsx_nb = self.header['signal_channels'][channel_indexes]['group_id'][0]
  542. if channel_indexes is None:
  543. local_indexes = slice(None)
  544. else:
  545. local_indexes = self._local_sig_indexes[channel_indexes]
  546. return nsx_nb, local_indexes
  547. def _get_signal_size(self, block_index, seg_index, channel_indexes):
  548. nsx_nb, local_indexes = self._get_nsx_and_local_indexes(channel_indexes)
  549. memmap_data = self.nsx_datas[nsx_nb][seg_index]
  550. return memmap_data.shape[0]
  551. def _get_signal_t_start(self, block_index, seg_index, channel_indexes):
  552. nsx_nb, local_indexes = self._get_nsx_and_local_indexes(channel_indexes)
  553. return self._sigs_t_starts[nsx_nb][seg_index]
  554. def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, channel_indexes):
  555. nsx_nb, local_indexes = self._get_nsx_and_local_indexes(channel_indexes)
  556. memmap_data = self.nsx_datas[nsx_nb][seg_index]
  557. sig_chunk = memmap_data[i_start:i_stop, local_indexes]
  558. return sig_chunk
  559. def _spike_count(self, block_index, seg_index, unit_index):
  560. channel_id, unit_id = self.internal_unit_ids[unit_index]
  561. all_spikes = self.nev_data['Spikes'][0]
  562. mask = (all_spikes['packet_id'] == channel_id) & (all_spikes['unit_class_nb'] == unit_id)
  563. if self._nb_segment == 1:
  564. # very fast
  565. nb = int(np.sum(mask))
  566. else:
  567. # must clip in time time range
  568. timestamp = all_spikes[mask]['timestamp']
  569. sl = self._get_timestamp_slice(timestamp, seg_index, None, None)
  570. timestamp = timestamp[sl]
  571. nb = timestamp.size
  572. return nb
  573. def _get_spike_timestamps(self, block_index, seg_index, unit_index, t_start, t_stop):
  574. channel_id, unit_id = self.internal_unit_ids[unit_index]
  575. all_spikes, event_segment_ids = self.nev_data['Spikes']
  576. # select by channel_id and unit_id
  577. mask = ((all_spikes['packet_id'] == channel_id) & (all_spikes['unit_class_nb'] == unit_id)
  578. & (event_segment_ids == seg_index))
  579. unit_spikes = all_spikes[mask]
  580. timestamp = unit_spikes['timestamp']
  581. sl = self._get_timestamp_slice(timestamp, seg_index, t_start, t_stop)
  582. timestamp = timestamp[sl]
  583. return timestamp
  584. def _get_timestamp_slice(self, timestamp, seg_index, t_start, t_stop):
  585. if self._nb_segment > 1:
  586. # we must clip event in seg time limits
  587. if t_start is None:
  588. t_start = self._seg_t_starts[seg_index]
  589. if t_stop is None:
  590. t_stop = self._seg_t_stops[seg_index]
  591. if t_start is None:
  592. ind_start = None
  593. else:
  594. ts = np.math.ceil(t_start * self.__nev_basic_header['timestamp_resolution'])
  595. ind_start = np.searchsorted(timestamp, ts)
  596. if t_stop is None:
  597. ind_stop = None
  598. else:
  599. ts = int(t_stop * self.__nev_basic_header['timestamp_resolution'])
  600. ind_stop = np.searchsorted(timestamp, ts) # +1
  601. sl = slice(ind_start, ind_stop)
  602. return sl
  603. def _rescale_spike_timestamp(self, spike_timestamps, dtype):
  604. spike_times = spike_timestamps.astype(dtype)
  605. spike_times /= self.__nev_basic_header['timestamp_resolution']
  606. return spike_times
  607. def _get_spike_raw_waveforms(self, block_index, seg_index, unit_index, t_start, t_stop):
  608. channel_id, unit_id = self.internal_unit_ids[unit_index]
  609. all_spikes, event_segment_ids = self.nev_data['Spikes']
  610. mask = ((all_spikes['packet_id'] == channel_id) & (all_spikes['unit_class_nb'] == unit_id)
  611. & (event_segment_ids == seg_index))
  612. unit_spikes = all_spikes[mask]
  613. wf_dtype = self.__nev_params('waveform_dtypes')[channel_id]
  614. wf_size = self.__nev_params('waveform_size')[channel_id]
  615. waveforms = unit_spikes['waveform'].flatten().view(wf_dtype)
  616. waveforms = waveforms.reshape(int(unit_spikes.size), 1, int(wf_size))
  617. timestamp = unit_spikes['timestamp']
  618. sl = self._get_timestamp_slice(timestamp, seg_index, t_start, t_stop)
  619. waveforms = waveforms[sl]
  620. return waveforms
  621. def _event_count(self, block_index, seg_index, event_channel_index):
  622. name = self.header['event_channels']['name'][event_channel_index]
  623. if name == 'comments':
  624. events_data, event_segment_ids = self.nev_data['Comments']
  625. ev_dict = self.__comment_evdict[self.__nev_spec](events_data)[name]
  626. else:
  627. events_data, event_segment_ids = self.nev_data['NonNeural']
  628. ev_dict = self.__nonneural_evdicts[self.__nev_spec](events_data)[name]
  629. mask = ev_dict['mask'] & (event_segment_ids == seg_index)
  630. if self._nb_segment == 1:
  631. # very fast
  632. nb = int(np.sum(mask))
  633. else:
  634. # must clip in time time range
  635. timestamp = events_data[ev_dict['mask']]['timestamp']
  636. sl = self._get_timestamp_slice(timestamp, seg_index, None, None)
  637. timestamp = timestamp[sl]
  638. nb = timestamp.size
  639. return nb
  640. def _get_event_timestamps(self, block_index, seg_index, event_channel_index, t_start, t_stop):
  641. name = self.header['event_channels']['name'][event_channel_index]
  642. if name == 'comments':
  643. events_data, event_segment_ids = self.nev_data['Comments']
  644. ev_dict = self.__comment_evdict[self.__nev_spec](events_data)[name]
  645. # If immediate decoding is desired:
  646. encoding = {0: 'latin_1', 1: 'utf_16', 255: 'latin_1'}
  647. labels = [data[ev_dict['field']].decode(
  648. encoding[data['char_set']]) for data in events_data]
  649. labels = np.array(labels, dtype='U')
  650. else:
  651. events_data, event_segment_ids = self.nev_data['NonNeural']
  652. ev_dict = self.__nonneural_evdicts[self.__nev_spec](events_data)[name]
  653. labels = events_data[ev_dict['field']].astype('U')
  654. mask = ev_dict['mask'] & (event_segment_ids == seg_index)
  655. timestamp = events_data[mask]['timestamp']
  656. labels = labels[mask]
  657. # time clip
  658. sl = self._get_timestamp_slice(timestamp, seg_index, t_start, t_stop)
  659. timestamp = timestamp[sl]
  660. labels = labels[sl]
  661. durations = None
  662. return timestamp, durations, labels
  663. def _rescale_event_timestamp(self, event_timestamps, dtype):
  664. ev_times = event_timestamps.astype(dtype)
  665. ev_times /= self.__nev_basic_header['timestamp_resolution']
  666. return ev_times
  667. ###################################################
  668. ###################################################
  669. # Above here code from Lyuba Zehl, Michael Denker
  670. # coming from previous BlackrockIO
  671. def __extract_nsx_file_spec(self, nsx_nb):
  672. """
  673. Extract file specification from an .nsx file.
  674. """
  675. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  676. # Header structure of files specification 2.2 and higher. For files 2.1
  677. # and lower, the entries ver_major and ver_minor are not supported.
  678. dt0 = [
  679. ('file_id', 'S8'),
  680. ('ver_major', 'uint8'),
  681. ('ver_minor', 'uint8')]
  682. nsx_file_id = np.fromfile(filename, count=1, dtype=dt0)[0]
  683. if nsx_file_id['file_id'].decode() == 'NEURALSG':
  684. spec = '2.1'
  685. elif nsx_file_id['file_id'].decode() == 'NEURALCD':
  686. spec = '{}.{}'.format(
  687. nsx_file_id['ver_major'], nsx_file_id['ver_minor'])
  688. else:
  689. raise IOError('Unsupported NSX file type.')
  690. return spec
  691. def __extract_nev_file_spec(self):
  692. """
  693. Extract file specification from an .nev file
  694. """
  695. filename = '.'.join([self._filenames['nev'], 'nev'])
  696. # Header structure of files specification 2.2 and higher. For files 2.1
  697. # and lower, the entries ver_major and ver_minor are not supported.
  698. dt0 = [
  699. ('file_id', 'S8'),
  700. ('ver_major', 'uint8'),
  701. ('ver_minor', 'uint8')]
  702. nev_file_id = np.fromfile(filename, count=1, dtype=dt0)[0]
  703. if nev_file_id['file_id'].decode() == 'NEURALEV':
  704. spec = '{}.{}'.format(
  705. nev_file_id['ver_major'], nev_file_id['ver_minor'])
  706. else:
  707. raise IOError('NEV file type {} is not supported'.format(
  708. nev_file_id['file_id']))
  709. return spec
  710. def __read_nsx_header_variant_a(self, nsx_nb):
  711. """
  712. Extract nsx header information from a 2.1 .nsx file
  713. """
  714. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  715. # basic header (file_id: NEURALCD)
  716. dt0 = [
  717. ('file_id', 'S8'),
  718. # label of sampling groun (e.g. "1kS/s" or "LFP Low")
  719. ('label', 'S16'),
  720. # number of 1/30000 seconds between data points
  721. # (e.g., if sampling rate "1 kS/s", period equals "30")
  722. ('period', 'uint32'),
  723. ('channel_count', 'uint32')]
  724. nsx_basic_header = np.fromfile(filename, count=1, dtype=dt0)[0]
  725. # "extended" header (last field of file_id: NEURALCD)
  726. # (to facilitate compatibility with higher file specs)
  727. offset_dt0 = np.dtype(dt0).itemsize
  728. shape = nsx_basic_header['channel_count']
  729. # originally called channel_id in Blackrock user manual
  730. # (to facilitate compatibility with higher file specs)
  731. dt1 = [('electrode_id', 'uint32')]
  732. nsx_ext_header = np.memmap(
  733. filename, shape=shape, offset=offset_dt0, dtype=dt1, mode='r')
  734. return nsx_basic_header, nsx_ext_header
  735. def __read_nsx_header_variant_b(self, nsx_nb):
  736. """
  737. Extract nsx header information from a 2.2 or 2.3 .nsx file
  738. """
  739. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  740. # basic header (file_id: NEURALCD)
  741. dt0 = [
  742. ('file_id', 'S8'),
  743. # file specification split into major and minor version number
  744. ('ver_major', 'uint8'),
  745. ('ver_minor', 'uint8'),
  746. # bytes of basic & extended header
  747. ('bytes_in_headers', 'uint32'),
  748. # label of the sampling group (e.g., "1 kS/s" or "LFP low")
  749. ('label', 'S16'),
  750. ('comment', 'S256'),
  751. ('period', 'uint32'),
  752. ('timestamp_resolution', 'uint32'),
  753. # time origin: 2byte uint16 values for ...
  754. ('year', 'uint16'),
  755. ('month', 'uint16'),
  756. ('weekday', 'uint16'),
  757. ('day', 'uint16'),
  758. ('hour', 'uint16'),
  759. ('minute', 'uint16'),
  760. ('second', 'uint16'),
  761. ('millisecond', 'uint16'),
  762. # number of channel_count match number of extended headers
  763. ('channel_count', 'uint32')]
  764. nsx_basic_header = np.fromfile(filename, count=1, dtype=dt0)[0]
  765. # extended header (type: CC)
  766. offset_dt0 = np.dtype(dt0).itemsize
  767. shape = nsx_basic_header['channel_count']
  768. dt1 = [
  769. ('type', 'S2'),
  770. ('electrode_id', 'uint16'),
  771. ('electrode_label', 'S16'),
  772. # used front-end amplifier bank (e.g., A, B, C, D)
  773. ('physical_connector', 'uint8'),
  774. # used connector pin (e.g., 1-37 on bank A, B, C or D)
  775. ('connector_pin', 'uint8'),
  776. # digital and analog value ranges of the signal
  777. ('min_digital_val', 'int16'),
  778. ('max_digital_val', 'int16'),
  779. ('min_analog_val', 'int16'),
  780. ('max_analog_val', 'int16'),
  781. # units of the analog range values ("mV" or "uV")
  782. ('units', 'S16'),
  783. # filter settings used to create nsx from source signal
  784. ('hi_freq_corner', 'uint32'),
  785. ('hi_freq_order', 'uint32'),
  786. ('hi_freq_type', 'uint16'), # 0=None, 1=Butterworth
  787. ('lo_freq_corner', 'uint32'),
  788. ('lo_freq_order', 'uint32'),
  789. ('lo_freq_type', 'uint16')] # 0=None, 1=Butterworth
  790. nsx_ext_header = np.memmap(
  791. filename, shape=shape, offset=offset_dt0, dtype=dt1, mode='r')
  792. return nsx_basic_header, nsx_ext_header
  793. def __read_nsx_dataheader(self, nsx_nb, offset):
  794. """
  795. Reads data header following the given offset of an nsx file.
  796. """
  797. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  798. # dtypes data header
  799. dt2 = [
  800. ('header', 'uint8'),
  801. ('timestamp', 'uint32'),
  802. ('nb_data_points', 'uint32')]
  803. return np.memmap(filename, dtype=dt2, shape=1, offset=offset, mode='r')[0]
  804. def __read_nsx_dataheader_variant_a(
  805. self, nsx_nb, filesize=None, offset=None):
  806. """
  807. Reads None for the nsx data header of file spec 2.1. Introduced to
  808. facilitate compatibility with higher file spec.
  809. """
  810. return None
  811. def __read_nsx_dataheader_variant_b(
  812. self, nsx_nb, filesize=None, offset=None, ):
  813. """
  814. Reads the nsx data header for each data block following the offset of
  815. file spec 2.2 and 2.3.
  816. """
  817. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  818. filesize = self.__get_file_size(filename)
  819. data_header = {}
  820. index = 0
  821. if offset is None:
  822. offset = self.__nsx_basic_header[nsx_nb]['bytes_in_headers']
  823. while offset < filesize:
  824. dh = self.__read_nsx_dataheader(nsx_nb, offset)
  825. data_header[index] = {
  826. 'header': dh['header'],
  827. 'timestamp': dh['timestamp'],
  828. 'nb_data_points': dh['nb_data_points'],
  829. 'offset_to_data_block': offset + dh.dtype.itemsize}
  830. # data size = number of data points * (2bytes * number of channels)
  831. # use of `int` avoids overflow problem
  832. data_size = int(dh['nb_data_points']) *\
  833. int(self.__nsx_basic_header[nsx_nb]['channel_count']) * 2
  834. # define new offset (to possible next data block)
  835. offset = data_header[index]['offset_to_data_block'] + data_size
  836. index += 1
  837. return data_header
  838. def __read_nsx_data_variant_a(self, nsx_nb):
  839. """
  840. Extract nsx data from a 2.1 .nsx file
  841. """
  842. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  843. # get shape of data
  844. shape = (
  845. self.__nsx_params['2.1'](nsx_nb)['nb_data_points'],
  846. self.__nsx_basic_header[nsx_nb]['channel_count'])
  847. offset = self.__nsx_params['2.1'](nsx_nb)['bytes_in_headers']
  848. # read nsx data
  849. # store as dict for compatibility with higher file specs
  850. data = {0: np.memmap(
  851. filename, dtype='int16', shape=shape, offset=offset, mode='r')}
  852. return data
  853. def __read_nsx_data_variant_b(self, nsx_nb):
  854. """
  855. Extract nsx data (blocks) from a 2.2 or 2.3 .nsx file. Blocks can arise
  856. if the recording was paused by the user.
  857. """
  858. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  859. data = {}
  860. for data_bl in self.__nsx_data_header[nsx_nb].keys():
  861. # get shape and offset of data
  862. shape = (
  863. self.__nsx_data_header[nsx_nb][data_bl]['nb_data_points'],
  864. self.__nsx_basic_header[nsx_nb]['channel_count'])
  865. offset = \
  866. self.__nsx_data_header[nsx_nb][data_bl]['offset_to_data_block']
  867. # read data
  868. data[data_bl] = np.memmap(
  869. filename, dtype='int16', shape=shape, offset=offset, mode='r')
  870. return data
  871. def __read_nev_header(self, ext_header_variants):
  872. """
  873. Extract nev header information from a of specific .nsx header variant
  874. """
  875. filename = '.'.join([self._filenames['nev'], 'nev'])
  876. # basic header
  877. dt0 = [
  878. # Set to "NEURALEV"
  879. ('file_type_id', 'S8'),
  880. ('ver_major', 'uint8'),
  881. ('ver_minor', 'uint8'),
  882. # Flags
  883. ('additionnal_flags', 'uint16'),
  884. # File index of first data sample
  885. ('bytes_in_headers', 'uint32'),
  886. # Number of bytes per data packet (sample)
  887. ('bytes_in_data_packets', 'uint32'),
  888. # Time resolution of time stamps in Hz
  889. ('timestamp_resolution', 'uint32'),
  890. # Sampling frequency of waveforms in Hz
  891. ('sample_resolution', 'uint32'),
  892. ('year', 'uint16'),
  893. ('month', 'uint16'),
  894. ('weekday', 'uint16'),
  895. ('day', 'uint16'),
  896. ('hour', 'uint16'),
  897. ('minute', 'uint16'),
  898. ('second', 'uint16'),
  899. ('millisecond', 'uint16'),
  900. ('application_to_create_file', 'S32'),
  901. ('comment_field', 'S256'),
  902. # Number of extended headers
  903. ('nb_ext_headers', 'uint32')]
  904. nev_basic_header = np.fromfile(filename, count=1, dtype=dt0)[0]
  905. # extended header
  906. # this consist in N block with code 8bytes + 24 data bytes
  907. # the data bytes depend on the code and need to be converted
  908. # cafilename_nsx, segse by case
  909. shape = nev_basic_header['nb_ext_headers']
  910. offset_dt0 = np.dtype(dt0).itemsize
  911. # This is the common structure of the beginning of extended headers
  912. dt1 = [
  913. ('packet_id', 'S8'),
  914. ('info_field', 'S24')]
  915. raw_ext_header = np.memmap(
  916. filename, offset=offset_dt0, dtype=dt1, shape=shape, mode='r')
  917. nev_ext_header = {}
  918. for packet_id in ext_header_variants.keys():
  919. mask = (raw_ext_header['packet_id'] == packet_id)
  920. dt2 = self.__nev_ext_header_types()[packet_id][
  921. ext_header_variants[packet_id]]
  922. nev_ext_header[packet_id] = raw_ext_header.view(dt2)[mask]
  923. return nev_basic_header, nev_ext_header
  924. def __read_nev_header_variant_a(self):
  925. """
  926. Extract nev header information from a 2.1 .nev file
  927. """
  928. ext_header_variants = {
  929. b'NEUEVWAV': 'a',
  930. b'ARRAYNME': 'a',
  931. b'ECOMMENT': 'a',
  932. b'CCOMMENT': 'a',
  933. b'MAPFILE': 'a',
  934. b'NSASEXEV': 'a'}
  935. return self.__read_nev_header(ext_header_variants)
  936. def __read_nev_header_variant_b(self):
  937. """
  938. Extract nev header information from a 2.2 .nev file
  939. """
  940. ext_header_variants = {
  941. b'NEUEVWAV': 'b',
  942. b'ARRAYNME': 'a',
  943. b'ECOMMENT': 'a',
  944. b'CCOMMENT': 'a',
  945. b'MAPFILE': 'a',
  946. b'NEUEVLBL': 'a',
  947. b'NEUEVFLT': 'a',
  948. b'DIGLABEL': 'a',
  949. b'NSASEXEV': 'a'}
  950. return self.__read_nev_header(ext_header_variants)
  951. def __read_nev_header_variant_c(self):
  952. """
  953. Extract nev header information from a 2.3 .nev file
  954. """
  955. ext_header_variants = {
  956. b'NEUEVWAV': 'b',
  957. b'ARRAYNME': 'a',
  958. b'ECOMMENT': 'a',
  959. b'CCOMMENT': 'a',
  960. b'MAPFILE': 'a',
  961. b'NEUEVLBL': 'a',
  962. b'NEUEVFLT': 'a',
  963. b'DIGLABEL': 'a',
  964. b'VIDEOSYN': 'a',
  965. b'TRACKOBJ': 'a'}
  966. return self.__read_nev_header(ext_header_variants)
  967. def __read_nev_data(self, nev_data_masks, nev_data_types):
  968. """
  969. Extract nev data from a 2.1 or 2.2 .nev file
  970. """
  971. filename = '.'.join([self._filenames['nev'], 'nev'])
  972. data_size = self.__nev_basic_header['bytes_in_data_packets']
  973. header_size = self.__nev_basic_header['bytes_in_headers']
  974. # read all raw data packets and markers
  975. dt0 = [
  976. ('timestamp', 'uint32'),
  977. ('packet_id', 'uint16'),
  978. ('value', 'S{}'.format(data_size - 6))]
  979. raw_data = np.memmap(filename, offset=header_size, dtype=dt0, mode='r')
  980. masks = self.__nev_data_masks(raw_data['packet_id'])
  981. types = self.__nev_data_types(data_size)
  982. event_segment_ids = self.__get_event_segment_ids(raw_data, masks, nev_data_masks)
  983. data = {}
  984. for k, v in nev_data_masks.items():
  985. mask = masks[k][v]
  986. data[k] = (raw_data.view(types[k][nev_data_types[k]])[mask], event_segment_ids[mask])
  987. return data
  988. def __get_reset_event_mask(self, raw_event_data, masks, nev_data_masks):
  989. """
  990. Extract mask for reset comment events in 2.3 .nev file
  991. """
  992. restart_mask = np.logical_and(masks['Comments'][nev_data_masks['Comments']],
  993. raw_event_data['value']
  994. == b'\x00\x00\x00\x00\x00\x00critical load restart')
  995. # TODO: Fix hardcoded number of bytes
  996. return restart_mask
  997. def __get_event_segment_ids(self, raw_event_data, masks, nev_data_masks):
  998. """
  999. Construct array of corresponding segment ids for each event for nev version 2.3
  1000. """
  1001. if self.__nev_spec in ['2.1', '2.2']:
  1002. # No pause or reset mechanism present for file version 2.1 and 2.2
  1003. return np.zeros(len(raw_event_data), dtype=int)
  1004. elif self.__nev_spec == '2.3':
  1005. reset_ev_mask = self.__get_reset_event_mask(raw_event_data, masks, nev_data_masks)
  1006. reset_ev_ids = np.where(reset_ev_mask)[0]
  1007. # consistency check for monotone increasing time stamps
  1008. # explicitely converting to int to allow for negative diff values
  1009. jump_ids = \
  1010. np.where(np.diff(np.asarray(raw_event_data['timestamp'], dtype=int)) < 0)[0] + 1
  1011. overlap = np.in1d(jump_ids, reset_ev_ids)
  1012. if not all(overlap):
  1013. # additional resets occurred without a reset event being stored
  1014. additional_ids = jump_ids[np.invert(overlap)]
  1015. warnings.warn('Detected {} undocumented segments within '
  1016. 'nev data after timestamps {}.'
  1017. ''.format(len(additional_ids), additional_ids))
  1018. reset_ev_ids = sorted(np.unique(np.concatenate((reset_ev_ids, jump_ids))))
  1019. event_segment_ids = np.zeros(len(raw_event_data), dtype=int)
  1020. for reset_event_id in reset_ev_ids:
  1021. event_segment_ids[reset_event_id:] += 1
  1022. self._nb_segment_nev = len(reset_ev_ids) + 1
  1023. return event_segment_ids
  1024. def __match_nsx_and_nev_segment_ids(self, nsx_nb):
  1025. """
  1026. Ensure matching ids of segments detected in nsx and nev file for version 2.3
  1027. """
  1028. # NSX required for matching, if not available, warn the user
  1029. if not self._avail_nsx:
  1030. warnings.warn("No nsX available so it cannot be checked whether "
  1031. "the segments in nev are all correct. Most importantly, "
  1032. "recording pauses will not be detected", UserWarning)
  1033. return
  1034. # Only needs to be done for nev version 2.3
  1035. if self.__nev_spec == '2.3':
  1036. nsx_offset = self.__nsx_data_header[nsx_nb][0]['timestamp']
  1037. # Multiples of 1/30.000s that pass between two nsX samples
  1038. nsx_period = self.__nsx_basic_header[nsx_nb]['period']
  1039. # NSX segments needed as dict and list
  1040. nonempty_nsx_segments = {}
  1041. list_nonempty_nsx_segments = []
  1042. # Counts how many segments CAN be created from nev
  1043. nb_possible_nev_segments = self._nb_segment_nev
  1044. # Nonempty segments are those containing at least 2 samples
  1045. # These have to be able to be mapped to nev
  1046. for k, v in sorted(self.__nsx_data_header[nsx_nb].items()):
  1047. if v['nb_data_points'] > 1:
  1048. nonempty_nsx_segments[k] = v
  1049. list_nonempty_nsx_segments.append(v)
  1050. # Account for paused segments
  1051. # This increases nev event segment ids if from the nsx an additional segment is found
  1052. # If one new segment, i.e. that could not be determined from the nev was found,
  1053. # all following ids need to be increased to account for the additional segment before
  1054. for k, (data, ev_ids) in self.nev_data.items():
  1055. # Check all nonempty nsX segments
  1056. for i, seg in enumerate(list_nonempty_nsx_segments[:]):
  1057. # Last timestamp in this nsX segment
  1058. # Not subtracting nsX offset from end because spike extraction might continue
  1059. end_of_current_nsx_seg = seg['timestamp'] + \
  1060. seg['nb_data_points'] * self.__nsx_basic_header[nsx_nb]['period']
  1061. mask_after_seg = (ev_ids == i) & \
  1062. (data['timestamp'] > end_of_current_nsx_seg + nsx_period)
  1063. # Show warning if spikes do not fit any segment (+- 1 sampling 'tick')
  1064. # Spike should belong to segment before
  1065. mask_outside = (ev_ids == i) & \
  1066. (data['timestamp'] < int(seg['timestamp']) - nsx_offset - nsx_period)
  1067. if len(data[mask_outside]) > 0:
  1068. warnings.warn("Spikes outside any segment. Detected on segment #{}".
  1069. format(i))
  1070. ev_ids[mask_outside] -= 1
  1071. # If some nev data are outside of this nsX segment, increase their segment ids
  1072. # and the ids of all following segments. They are checked for the next nsX
  1073. # segment then. If they do not fit any of them,
  1074. # a warning will be shown, indicating how far outside the segment spikes are
  1075. # If they fit the next segment, more segments are possible in nev,
  1076. # because a new one has been discovered
  1077. if len(data[mask_after_seg]) > 0:
  1078. # Warning if spikes are after last segment
  1079. if i == len(list_nonempty_nsx_segments) - 1:
  1080. timestamp_resolution = self.__nsx_params[self.__nsx_spec[
  1081. nsx_nb]]('timestamp_resolution', nsx_nb)
  1082. time_after_seg = (data[mask_after_seg]['timestamp'][-1]
  1083. - end_of_current_nsx_seg) / timestamp_resolution
  1084. warnings.warn("Spikes {}s after last segment.".format(time_after_seg))
  1085. # Break out of loop because it's the last iteration
  1086. # and the spikes should stay connected to last segment
  1087. break
  1088. # If reset and no segment detected in nev, then these segments cannot be
  1089. # distinguished in nev, which is a big problem
  1090. # XXX 96 is an arbitrary number based on observations in available files
  1091. elif list_nonempty_nsx_segments[i + 1]['timestamp'] - nsx_offset <= 96:
  1092. # If not all definitely belong to the next segment,
  1093. # then it cannot be distinguished where some belong
  1094. if len(data[ev_ids == i]) != len(data[mask_after_seg]):
  1095. raise ValueError("Some segments in nsX cannot be detected in nev")
  1096. # Actual processing if no problem has occurred
  1097. nb_possible_nev_segments += 1
  1098. ev_ids[ev_ids > i] += 1
  1099. ev_ids[mask_after_seg] += 1
  1100. # consistency check: same number of segments for nsx and nev data
  1101. assert nb_possible_nev_segments == len(nonempty_nsx_segments), \
  1102. ('Inconsistent ns{0} and nev file. {1} segments present in .nev file, but {2} in '
  1103. 'ns{0} file.'.format(nsx_nb, nb_possible_nev_segments,
  1104. len(nonempty_nsx_segments)))
  1105. new_nev_segment_id_mapping = dict(zip(range(nb_possible_nev_segments),
  1106. sorted(list(nonempty_nsx_segments))))
  1107. # replacing event ids by matched event ids in place
  1108. for k, (data, ev_ids) in self.nev_data.items():
  1109. if len(ev_ids):
  1110. ev_ids[:] = np.vectorize(new_nev_segment_id_mapping.__getitem__)(ev_ids)
  1111. def __read_nev_data_variant_a(self):
  1112. """
  1113. Extract nev data from a 2.1 & 2.2 .nev file
  1114. """
  1115. nev_data_masks = {
  1116. 'NonNeural': 'a',
  1117. 'Spikes': 'a'}
  1118. nev_data_types = {
  1119. 'NonNeural': 'a',
  1120. 'Spikes': 'a'}
  1121. return self.__read_nev_data(nev_data_masks, nev_data_types)
  1122. def __read_nev_data_variant_b(self):
  1123. """
  1124. Extract nev data from a 2.3 .nev file
  1125. """
  1126. nev_data_masks = {
  1127. 'NonNeural': 'a',
  1128. 'Spikes': 'b',
  1129. 'Comments': 'a',
  1130. 'VideoSync': 'a',
  1131. 'TrackingEvents': 'a',
  1132. 'ButtonTrigger': 'a',
  1133. 'ConfigEvent': 'a'}
  1134. nev_data_types = {
  1135. 'NonNeural': 'b',
  1136. 'Spikes': 'a',
  1137. 'Comments': 'a',
  1138. 'VideoSync': 'a',
  1139. 'TrackingEvents': 'a',
  1140. 'ButtonTrigger': 'a',
  1141. 'ConfigEvent': 'a'}
  1142. return self.__read_nev_data(nev_data_masks, nev_data_types)
  1143. def __nev_ext_header_types(self):
  1144. """
  1145. Defines extended header types for different .nev file specifications.
  1146. """
  1147. nev_ext_header_types = {
  1148. b'NEUEVWAV': {
  1149. # Version>=2.1
  1150. 'a': [
  1151. ('packet_id', 'S8'),
  1152. ('electrode_id', 'uint16'),
  1153. ('physical_connector', 'uint8'),
  1154. ('connector_pin', 'uint8'),
  1155. ('digitization_factor', 'uint16'),
  1156. ('energy_threshold', 'uint16'),
  1157. ('hi_threshold', 'int16'),
  1158. ('lo_threshold', 'int16'),
  1159. ('nb_sorted_units', 'uint8'),
  1160. # number of bytes per waveform sample
  1161. ('bytes_per_waveform', 'uint8'),
  1162. ('unused', 'S10')],
  1163. # Version>=2.3
  1164. 'b': [
  1165. ('packet_id', 'S8'),
  1166. ('electrode_id', 'uint16'),
  1167. ('physical_connector', 'uint8'),
  1168. ('connector_pin', 'uint8'),
  1169. ('digitization_factor', 'uint16'),
  1170. ('energy_threshold', 'uint16'),
  1171. ('hi_threshold', 'int16'),
  1172. ('lo_threshold', 'int16'),
  1173. ('nb_sorted_units', 'uint8'),
  1174. # number of bytes per waveform sample
  1175. ('bytes_per_waveform', 'uint8'),
  1176. # number of samples for each waveform
  1177. ('spike_width', 'uint16'),
  1178. ('unused', 'S8')]},
  1179. b'ARRAYNME': {
  1180. 'a': [
  1181. ('packet_id', 'S8'),
  1182. ('electrode_array_name', 'S24')]},
  1183. b'ECOMMENT': {
  1184. 'a': [
  1185. ('packet_id', 'S8'),
  1186. ('extra_comment', 'S24')]},
  1187. b'CCOMMENT': {
  1188. 'a': [
  1189. ('packet_id', 'S8'),
  1190. ('continued_comment', 'S24')]},
  1191. b'MAPFILE': {
  1192. 'a': [
  1193. ('packet_id', 'S8'),
  1194. ('mapFile', 'S24')]},
  1195. b'NEUEVLBL': {
  1196. 'a': [
  1197. ('packet_id', 'S8'),
  1198. ('electrode_id', 'uint16'),
  1199. # label of this electrode
  1200. ('label', 'S16'),
  1201. ('unused', 'S6')]},
  1202. b'NEUEVFLT': {
  1203. 'a': [
  1204. ('packet_id', 'S8'),
  1205. ('electrode_id', 'uint16'),
  1206. ('hi_freq_corner', 'uint32'),
  1207. ('hi_freq_order', 'uint32'),
  1208. # 0=None 1=Butterworth
  1209. ('hi_freq_type', 'uint16'),
  1210. ('lo_freq_corner', 'uint32'),
  1211. ('lo_freq_order', 'uint32'),
  1212. # 0=None 1=Butterworth
  1213. ('lo_freq_type', 'uint16'),
  1214. ('unused', 'S2')]},
  1215. b'DIGLABEL': {
  1216. 'a': [
  1217. ('packet_id', 'S8'),
  1218. # Read name of digital
  1219. ('label', 'S16'),
  1220. # 0=serial, 1=parallel
  1221. ('mode', 'uint8'),
  1222. ('unused', 'S7')]},
  1223. b'NSASEXEV': {
  1224. 'a': [
  1225. ('packet_id', 'S8'),
  1226. # Read frequency of periodic packet generation
  1227. ('frequency', 'uint16'),
  1228. # Read if digital input triggers events
  1229. ('digital_input_config', 'uint8'),
  1230. # Read if analog input triggers events
  1231. ('analog_channel_1_config', 'uint8'),
  1232. ('analog_channel_1_edge_detec_val', 'uint16'),
  1233. ('analog_channel_2_config', 'uint8'),
  1234. ('analog_channel_2_edge_detec_val', 'uint16'),
  1235. ('analog_channel_3_config', 'uint8'),
  1236. ('analog_channel_3_edge_detec_val', 'uint16'),
  1237. ('analog_channel_4_config', 'uint8'),
  1238. ('analog_channel_4_edge_detec_val', 'uint16'),
  1239. ('analog_channel_5_config', 'uint8'),
  1240. ('analog_channel_5_edge_detec_val', 'uint16'),
  1241. ('unused', 'S6')]},
  1242. b'VIDEOSYN': {
  1243. 'a': [
  1244. ('packet_id', 'S8'),
  1245. ('video_source_id', 'uint16'),
  1246. ('video_source', 'S16'),
  1247. ('frame_rate', 'float32'),
  1248. ('unused', 'S2')]},
  1249. b'TRACKOBJ': {
  1250. 'a': [
  1251. ('packet_id', 'S8'),
  1252. ('trackable_type', 'uint16'),
  1253. ('trackable_id', 'uint16'),
  1254. ('point_count', 'uint16'),
  1255. ('video_source', 'S16'),
  1256. ('unused', 'S2')]}}
  1257. return nev_ext_header_types
  1258. def __nev_data_masks(self, packet_ids):
  1259. """
  1260. Defines data masks for different .nev file specifications depending on
  1261. the given packet identifiers.
  1262. """
  1263. __nev_data_masks = {
  1264. 'NonNeural': {
  1265. 'a': (packet_ids == 0)},
  1266. 'Spikes': {
  1267. # Version 2.1 & 2.2
  1268. 'a': (0 < packet_ids) & (packet_ids <= 255),
  1269. # Version>=2.3
  1270. 'b': (0 < packet_ids) & (packet_ids <= 2048)},
  1271. 'Comments': {
  1272. 'a': (packet_ids == 0xFFFF)},
  1273. 'VideoSync': {
  1274. 'a': (packet_ids == 0xFFFE)},
  1275. 'TrackingEvents': {
  1276. 'a': (packet_ids == 0xFFFD)},
  1277. 'ButtonTrigger': {
  1278. 'a': (packet_ids == 0xFFFC)},
  1279. 'ConfigEvent': {
  1280. 'a': (packet_ids == 0xFFFB)}}
  1281. return __nev_data_masks
  1282. def __nev_data_types(self, data_size):
  1283. """
  1284. Defines data types for different .nev file specifications depending on
  1285. the given packet identifiers.
  1286. """
  1287. __nev_data_types = {
  1288. 'NonNeural': {
  1289. # Version 2.1 & 2.2
  1290. 'a': [
  1291. ('timestamp', 'uint32'),
  1292. ('packet_id', 'uint16'),
  1293. ('packet_insertion_reason', 'uint8'),
  1294. ('reserved', 'uint8'),
  1295. ('digital_input', 'uint16'),
  1296. ('analog_input_channel_1', 'int16'),
  1297. ('analog_input_channel_2', 'int16'),
  1298. ('analog_input_channel_3', 'int16'),
  1299. ('analog_input_channel_4', 'int16'),
  1300. ('analog_input_channel_5', 'int16'),
  1301. ('unused', 'S{}'.format(data_size - 20))],
  1302. # Version>=2.3
  1303. 'b': [
  1304. ('timestamp', 'uint32'),
  1305. ('packet_id', 'uint16'),
  1306. ('packet_insertion_reason', 'uint8'),
  1307. ('reserved', 'uint8'),
  1308. ('digital_input', 'uint16'),
  1309. ('unused', 'S{}'.format(data_size - 10))]},
  1310. 'Spikes': {
  1311. 'a': [
  1312. ('timestamp', 'uint32'),
  1313. ('packet_id', 'uint16'),
  1314. ('unit_class_nb', 'uint8'),
  1315. ('reserved', 'uint8'),
  1316. ('waveform', 'S{}'.format(data_size - 8))]},
  1317. 'Comments': {
  1318. 'a': [
  1319. ('timestamp', 'uint32'),
  1320. ('packet_id', 'uint16'),
  1321. ('char_set', 'uint8'),
  1322. ('flag', 'uint8'),
  1323. ('color', 'uint32'),
  1324. ('comment', 'S{}'.format(data_size - 12))]},
  1325. 'VideoSync': {
  1326. 'a': [
  1327. ('timestamp', 'uint32'),
  1328. ('packet_id', 'uint16'),
  1329. ('video_file_nb', 'uint16'),
  1330. ('video_frame_nb', 'uint32'),
  1331. ('video_elapsed_time', 'uint32'),
  1332. ('video_source_id', 'uint32'),
  1333. ('unused', 'int8', (data_size - 20,))]},
  1334. 'TrackingEvents': {
  1335. 'a': [
  1336. ('timestamp', 'uint32'),
  1337. ('packet_id', 'uint16'),
  1338. ('parent_id', 'uint16'),
  1339. ('node_id', 'uint16'),
  1340. ('node_count', 'uint16'),
  1341. ('point_count', 'uint16'),
  1342. ('tracking_points', 'uint16', ((data_size - 14) // 2,))]},
  1343. 'ButtonTrigger': {
  1344. 'a': [
  1345. ('timestamp', 'uint32'),
  1346. ('packet_id', 'uint16'),
  1347. ('trigger_type', 'uint16'),
  1348. ('unused', 'int8', (data_size - 8,))]},
  1349. 'ConfigEvent': {
  1350. 'a': [
  1351. ('timestamp', 'uint32'),
  1352. ('packet_id', 'uint16'),
  1353. ('config_change_type', 'uint16'),
  1354. ('config_changed', 'S{}'.format(data_size - 8))]}}
  1355. return __nev_data_types
  1356. def __nev_params(self, param_name):
  1357. """
  1358. Returns wanted nev parameter.
  1359. """
  1360. nev_parameters = {
  1361. 'bytes_in_data_packets':
  1362. self.__nev_basic_header['bytes_in_data_packets'],
  1363. 'rec_datetime': datetime.datetime(
  1364. year=self.__nev_basic_header['year'],
  1365. month=self.__nev_basic_header['month'],
  1366. day=self.__nev_basic_header['day'],
  1367. hour=self.__nev_basic_header['hour'],
  1368. minute=self.__nev_basic_header['minute'],
  1369. second=self.__nev_basic_header['second'],
  1370. microsecond=self.__nev_basic_header['millisecond']),
  1371. 'max_res': self.__nev_basic_header['timestamp_resolution'],
  1372. 'channel_ids': self.__nev_ext_header[b'NEUEVWAV']['electrode_id'],
  1373. 'channel_labels': self.__channel_labels[self.__nev_spec](),
  1374. 'event_unit': pq.CompoundUnit("1.0/{} * s".format(
  1375. self.__nev_basic_header['timestamp_resolution'])),
  1376. 'nb_units': dict(zip(
  1377. self.__nev_ext_header[b'NEUEVWAV']['electrode_id'],
  1378. self.__nev_ext_header[b'NEUEVWAV']['nb_sorted_units'])),
  1379. 'digitization_factor': dict(zip(
  1380. self.__nev_ext_header[b'NEUEVWAV']['electrode_id'],
  1381. self.__nev_ext_header[b'NEUEVWAV']['digitization_factor'])),
  1382. 'data_size': self.__nev_basic_header['bytes_in_data_packets'],
  1383. 'waveform_size': self.__waveform_size[self.__nev_spec](),
  1384. 'waveform_dtypes': self.__get_waveforms_dtype(),
  1385. 'waveform_sampling_rate':
  1386. self.__nev_basic_header['sample_resolution'] * pq.Hz,
  1387. 'waveform_time_unit': pq.CompoundUnit("1.0/{} * s".format(
  1388. self.__nev_basic_header['sample_resolution'])),
  1389. 'waveform_unit': pq.uV}
  1390. return nev_parameters[param_name]
  1391. def __get_file_size(self, filename):
  1392. """
  1393. Returns the file size in bytes for the given file.
  1394. """
  1395. filebuf = open(filename, 'rb')
  1396. filebuf.seek(0, os.SEEK_END)
  1397. file_size = filebuf.tell()
  1398. filebuf.close()
  1399. return file_size
  1400. def __get_min_time(self):
  1401. """
  1402. Returns the smallest time that can be determined from the recording for
  1403. use as the lower bound n in an interval [n,m).
  1404. """
  1405. tp = []
  1406. if self._avail_files['nev']:
  1407. tp.extend(self.__get_nev_rec_times()[0])
  1408. for nsx_i in self._avail_nsx:
  1409. tp.extend(self.__nsx_rec_times[self.__nsx_spec[nsx_i]](nsx_i)[0])
  1410. return min(tp)
  1411. def __get_max_time(self):
  1412. """
  1413. Returns the largest time that can be determined from the recording for
  1414. use as the upper bound m in an interval [n,m).
  1415. """
  1416. tp = []
  1417. if self._avail_files['nev']:
  1418. tp.extend(self.__get_nev_rec_times()[1])
  1419. for nsx_i in self._avail_nsx:
  1420. tp.extend(self.__nsx_rec_times[self.__nsx_spec[nsx_i]](nsx_i)[1])
  1421. return max(tp)
  1422. def __get_nev_rec_times(self):
  1423. """
  1424. Extracts minimum and maximum time points from a nev file.
  1425. """
  1426. filename = '.'.join([self._filenames['nev'], 'nev'])
  1427. dt = [('timestamp', 'uint32')]
  1428. offset = \
  1429. self.__get_file_size(filename) - \
  1430. self.__nev_params('bytes_in_data_packets')
  1431. last_data_packet = np.memmap(filename, offset=offset, dtype=dt, mode='r')[0]
  1432. n_starts = [0 * self.__nev_params('event_unit')]
  1433. n_stops = [
  1434. last_data_packet['timestamp'] * self.__nev_params('event_unit')]
  1435. return n_starts, n_stops
  1436. def __get_waveforms_dtype(self):
  1437. """
  1438. Extracts the actual waveform dtype set for each channel.
  1439. """
  1440. # Blackrock code giving the approiate dtype
  1441. conv = {0: 'int8', 1: 'int8', 2: 'int16', 4: 'int32'}
  1442. # get all electrode ids from nev ext header
  1443. all_el_ids = self.__nev_ext_header[b'NEUEVWAV']['electrode_id']
  1444. # get the dtype of waveform (this is stupidly complicated)
  1445. if self.__is_set(
  1446. np.array(self.__nev_basic_header['additionnal_flags']), 0):
  1447. dtype_waveforms = {k: 'int16' for k in all_el_ids}
  1448. else:
  1449. # extract bytes per waveform
  1450. waveform_bytes = \
  1451. self.__nev_ext_header[b'NEUEVWAV']['bytes_per_waveform']
  1452. # extract dtype for waveforms fro each electrode
  1453. dtype_waveforms = dict(zip(all_el_ids, conv[waveform_bytes]))
  1454. return dtype_waveforms
  1455. def __get_channel_labels_variant_a(self):
  1456. """
  1457. Returns labels for all channels for file spec 2.1
  1458. """
  1459. elids = self.__nev_ext_header[b'NEUEVWAV']['electrode_id']
  1460. labels = []
  1461. for elid in elids:
  1462. if elid < 129:
  1463. labels.append('chan%i' % elid)
  1464. else:
  1465. labels.append('ainp%i' % (elid - 129 + 1))
  1466. return dict(zip(elids, labels))
  1467. def __get_channel_labels_variant_b(self):
  1468. """
  1469. Returns labels for all channels for file spec 2.2 and 2.3
  1470. """
  1471. elids = self.__nev_ext_header[b'NEUEVWAV']['electrode_id']
  1472. labels = self.__nev_ext_header[b'NEUEVLBL']['label']
  1473. return dict(zip(elids, labels)) if len(labels) > 0 else None
  1474. def __get_waveform_size_variant_a(self):
  1475. """
  1476. Returns wavform sizes for all channels for file spec 2.1 and 2.2
  1477. """
  1478. wf_dtypes = self.__get_waveforms_dtype()
  1479. nb_bytes_wf = self.__nev_basic_header['bytes_in_data_packets'] - 8
  1480. wf_sizes = {
  1481. ch: int(nb_bytes_wf / np.dtype(dt).itemsize) for ch, dt in
  1482. wf_dtypes.items()}
  1483. return wf_sizes
  1484. def __get_waveform_size_variant_b(self):
  1485. """
  1486. Returns wavform sizes for all channels for file spec 2.3
  1487. """
  1488. elids = self.__nev_ext_header[b'NEUEVWAV']['electrode_id']
  1489. spike_widths = self.__nev_ext_header[b'NEUEVWAV']['spike_width']
  1490. return dict(zip(elids, spike_widths))
  1491. def __get_left_sweep_waveforms(self):
  1492. """
  1493. Returns left sweep of waveforms for each channel. Left sweep is defined
  1494. as the time from the beginning of the waveform to the trigger time of
  1495. the corresponding spike.
  1496. """
  1497. # TODO: Double check if this is the actual setting for Blackrock
  1498. wf_t_unit = self.__nev_params('waveform_time_unit')
  1499. all_ch = self.__nev_params('channel_ids')
  1500. # TODO: Double check if this is the correct assumption (10 samples)
  1501. # default value: threshold crossing after 10 samples of waveform
  1502. wf_left_sweep = {ch: 10 * wf_t_unit for ch in all_ch}
  1503. # non-default: threshold crossing at center of waveform
  1504. # wf_size = self.__nev_params('waveform_size')
  1505. # wf_left_sweep = dict(
  1506. # [(ch, (wf_size[ch] / 2) * wf_t_unit) for ch in all_ch])
  1507. return wf_left_sweep
  1508. def __get_nsx_param_variant_a(self, nsx_nb):
  1509. """
  1510. Returns parameter (param_name) for a given nsx (nsx_nb) for file spec
  1511. 2.1.
  1512. """
  1513. # Here, min/max_analog_val and min/max_digital_val are not available in
  1514. # the nsx, so that we must estimate these parameters from the
  1515. # digitization factor of the nev (information by Kian Torab, Blackrock
  1516. # Microsystems). Here dig_factor=max_analog_val/max_digital_val. We set
  1517. # max_digital_val to 1000, and max_analog_val=dig_factor. dig_factor is
  1518. # given in nV by definition, so the units turn out to be uV.
  1519. labels = []
  1520. dig_factor = []
  1521. for elid in self.__nsx_ext_header[nsx_nb]['electrode_id']:
  1522. if self._avail_files['nev']:
  1523. # This is a workaround for the DigitalFactor overflow in NEV
  1524. # files recorded with buggy Cerebus system.
  1525. # Fix taken from: NMPK toolbox by Blackrock,
  1526. # file openNEV, line 464,
  1527. # git rev. d0a25eac902704a3a29fa5dfd3aed0744f4733ed
  1528. df = self.__nev_params('digitization_factor')[elid]
  1529. if df == 21516:
  1530. df = 152592.547
  1531. dig_factor.append(df)
  1532. else:
  1533. dig_factor.append(float('nan'))
  1534. if elid < 129:
  1535. labels.append('chan%i' % elid)
  1536. else:
  1537. labels.append('ainp%i' % (elid - 129 + 1))
  1538. filename = '.'.join([self._filenames['nsx'], 'ns%i' % nsx_nb])
  1539. bytes_in_headers = self.__nsx_basic_header[nsx_nb].dtype.itemsize + \
  1540. self.__nsx_ext_header[nsx_nb].dtype.itemsize * \
  1541. self.__nsx_basic_header[nsx_nb]['channel_count']
  1542. if np.isnan(dig_factor[0]):
  1543. units = ''
  1544. warnings.warn("Cannot rescale to voltage, raw data will be returned.", UserWarning)
  1545. else:
  1546. units = 'uV'
  1547. nsx_parameters = {
  1548. 'nb_data_points': int(
  1549. (self.__get_file_size(filename) - bytes_in_headers)
  1550. / (2 * self.__nsx_basic_header[nsx_nb]['channel_count']) - 1),
  1551. 'labels': labels,
  1552. 'units': np.array([units] * self.__nsx_basic_header[nsx_nb]['channel_count']),
  1553. 'min_analog_val': -1 * np.array(dig_factor),
  1554. 'max_analog_val': np.array(dig_factor),
  1555. 'min_digital_val': np.array(
  1556. [-1000] * self.__nsx_basic_header[nsx_nb]['channel_count']),
  1557. 'max_digital_val': np.array([1000] * self.__nsx_basic_header[nsx_nb]['channel_count']),
  1558. 'timestamp_resolution': 30000,
  1559. 'bytes_in_headers': bytes_in_headers,
  1560. 'sampling_rate': 30000 / self.__nsx_basic_header[nsx_nb]['period'] * pq.Hz,
  1561. 'time_unit': pq.CompoundUnit("1.0/{}*s".format(
  1562. 30000 / self.__nsx_basic_header[nsx_nb]['period']))}
  1563. # Returns complete dictionary because then it does not need to be called so often
  1564. return nsx_parameters
  1565. def __get_nsx_param_variant_b(self, param_name, nsx_nb):
  1566. """
  1567. Returns parameter (param_name) for a given nsx (nsx_nb) for file spec
  1568. 2.2 and 2.3.
  1569. """
  1570. nsx_parameters = {
  1571. 'labels':
  1572. self.__nsx_ext_header[nsx_nb]['electrode_label'],
  1573. 'units':
  1574. self.__nsx_ext_header[nsx_nb]['units'],
  1575. 'min_analog_val':
  1576. self.__nsx_ext_header[nsx_nb]['min_analog_val'],
  1577. 'max_analog_val':
  1578. self.__nsx_ext_header[nsx_nb]['max_analog_val'],
  1579. 'min_digital_val':
  1580. self.__nsx_ext_header[nsx_nb]['min_digital_val'],
  1581. 'max_digital_val':
  1582. self.__nsx_ext_header[nsx_nb]['max_digital_val'],
  1583. 'timestamp_resolution':
  1584. self.__nsx_basic_header[nsx_nb]['timestamp_resolution'],
  1585. 'bytes_in_headers':
  1586. self.__nsx_basic_header[nsx_nb]['bytes_in_headers'],
  1587. 'sampling_rate':
  1588. self.__nsx_basic_header[nsx_nb]['timestamp_resolution']
  1589. / self.__nsx_basic_header[nsx_nb]['period'] * pq.Hz,
  1590. 'time_unit': pq.CompoundUnit("1.0/{}*s".format(
  1591. self.__nsx_basic_header[nsx_nb]['timestamp_resolution']
  1592. / self.__nsx_basic_header[nsx_nb]['period']))}
  1593. return nsx_parameters[param_name]
  1594. def __get_nonneural_evdicts_variant_a(self, data):
  1595. """
  1596. Defines event types and the necessary parameters to extract them from
  1597. a 2.1 and 2.2 nev file.
  1598. """
  1599. # TODO: add annotations of nev ext header (NSASEXEX) to event types
  1600. # digital events
  1601. event_types = {
  1602. 'digital_input_port': {
  1603. 'name': 'digital_input_port',
  1604. 'field': 'digital_input',
  1605. 'mask': data['packet_insertion_reason'] == 1,
  1606. 'desc': "Events of the digital input port"},
  1607. 'serial_input_port': {
  1608. 'name': 'serial_input_port',
  1609. 'field': 'digital_input',
  1610. 'mask': data['packet_insertion_reason'] == 129,
  1611. 'desc': "Events of the serial input port"}}
  1612. # analog input events via threshold crossings
  1613. for ch in range(5):
  1614. event_types.update({
  1615. 'analog_input_channel_{}'.format(ch + 1): {
  1616. 'name': 'analog_input_channel_{}'.format(ch + 1),
  1617. 'field': 'analog_input_channel_{}'.format(ch + 1),
  1618. 'mask': self.__is_set(
  1619. data['packet_insertion_reason'], ch + 1),
  1620. 'desc': "Values of analog input channel {} in mV "
  1621. "(+/- 5000)".format(ch + 1)}})
  1622. # TODO: define field and desc
  1623. event_types.update({
  1624. 'periodic_sampling_events': {
  1625. 'name': 'periodic_sampling_events',
  1626. 'field': 'digital_input',
  1627. 'mask': self.__is_set(data['packet_insertion_reason'], 6),
  1628. 'desc': 'Periodic sampling event of a certain frequency'}})
  1629. return event_types
  1630. def __delete_empty_segments(self):
  1631. """
  1632. If there are empty segments (e.g. due to a reset or clock synchronization across
  1633. two systems), these can be discarded.
  1634. If this is done, all the data and data_headers need to be remapped to become a range
  1635. starting from 0 again. Nev data are mapped accordingly to stay with their corresponding
  1636. segment in the nsX data.
  1637. """
  1638. # Discard empty segments
  1639. removed_seg = []
  1640. for data_bl in range(self._nb_segment):
  1641. keep_seg = True
  1642. for nsx_nb in self.nsx_to_load:
  1643. length = self.nsx_datas[nsx_nb][data_bl].shape[0]
  1644. keep_seg = keep_seg and (length >= 2)
  1645. if not keep_seg:
  1646. removed_seg.append(data_bl)
  1647. for nsx_nb in self.nsx_to_load:
  1648. self.nsx_datas[nsx_nb].pop(data_bl)
  1649. self.__nsx_data_header[nsx_nb].pop(data_bl)
  1650. # Keys need to be increasing from 0 to maximum in steps of 1
  1651. # To ensure this after removing empty segments, some keys need to be re mapped
  1652. for i in removed_seg[::-1]:
  1653. for j in range(i + 1, self._nb_segment):
  1654. # remap nsx seg index
  1655. for nsx_nb in self.nsx_to_load:
  1656. data = self.nsx_datas[nsx_nb].pop(j)
  1657. self.nsx_datas[nsx_nb][j - 1] = data
  1658. data_header = self.__nsx_data_header[nsx_nb].pop(j)
  1659. self.__nsx_data_header[nsx_nb][j - 1] = data_header
  1660. # Also remap nev data, ev_ids are the equivalent to keys above
  1661. if self._avail_files['nev']:
  1662. for k, (data, ev_ids) in self.nev_data.items():
  1663. ev_ids[ev_ids == j] -= 1
  1664. self._nb_segment -= 1
  1665. def __get_nonneural_evdicts_variant_b(self, data):
  1666. """
  1667. Defines event types and the necessary parameters to extract them from
  1668. a 2.3 nev file.
  1669. """
  1670. # digital events
  1671. if not np.all(np.in1d(data['packet_insertion_reason'], [1, 129])):
  1672. # Blackrock spec gives reason==64 means PERIODIC, but never seen this live
  1673. warnings.warn("Unknown event codes found", RuntimeWarning)
  1674. event_types = {
  1675. 'digital_input_port': {
  1676. 'name': 'digital_input_port',
  1677. 'field': 'digital_input',
  1678. 'mask': self.__is_set(data['packet_insertion_reason'], 0)
  1679. & ~self.__is_set(data['packet_insertion_reason'], 7),
  1680. 'desc': "Events of the digital input port"},
  1681. 'serial_input_port': {
  1682. 'name': 'serial_input_port',
  1683. 'field': 'digital_input',
  1684. 'mask': self.__is_set(data['packet_insertion_reason'], 0)
  1685. & self.__is_set(data['packet_insertion_reason'], 7),
  1686. 'desc': "Events of the serial input port"}}
  1687. return event_types
  1688. def __get_comment_evdict_variant_a(self, data):
  1689. return {
  1690. 'comments': {
  1691. 'name': 'comments',
  1692. 'field': 'comment',
  1693. 'mask': data['packet_id'] == 65535,
  1694. 'desc': 'Comments'
  1695. }
  1696. }
  1697. def __is_set(self, flag, pos):
  1698. """
  1699. Checks if bit is set at the given position for flag. If flag is an
  1700. array, an array will be returned.
  1701. """
  1702. return flag & (1 << pos) > 0