From b1120471c580a9a97ea8bd2e727f312441420e1d Mon Sep 17 00:00:00 2001 From: Matteo Caffini Date: Wed, 10 May 2017 13:08:06 +0200 Subject: [PATCH 1/5] ENH: Added tutorial for fNIRS data --- tutorials/Working_with_fNIRS_data.py | 124 +++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 tutorials/Working_with_fNIRS_data.py diff --git a/tutorials/Working_with_fNIRS_data.py b/tutorials/Working_with_fNIRS_data.py new file mode 100644 index 00000000000..e2851368d4e --- /dev/null +++ b/tutorials/Working_with_fNIRS_data.py @@ -0,0 +1,124 @@ +""" +.. _tut_fnirs_data + +Use MNE with fNIRS data +======================= +# Author: Matteo Caffini +# Contact: + +""" +import numpy as np + +import mne + +raw_ndarray = np.load('/home/matteo/raw_ndarray.npy') + +############################################################################### +# MNE Python was implemented with EEG and MEG in mind, but can be useful to +# process fNIRS data as #well. A common fNIRS data set (similarly to EEG and +# MEG data) contains data recorded over time(a.k.a. samples) at many locations +# on the scalp (a.k.a. channels). +# Usually in fNIRS one looks for brain responses to external sensorial +# stimulations (visual, acoustic, etc.). In order to do so, block design +# studies are usually employed and data are later epoched to get evoked +# responses. + +############################################################################### +# fNIRS raw data come in various formats and the best way to dive into MNE is +# to deal with your own format first and fit your data into a 2D numpy array +# [channels x samples]. Once you prepared your time series for HbO and HbR +# concentrations, you can create a RawArray instance from the data array using +# mne.io.RawArray. If you have a time series of events synchronous with your +# data you can append this as a last row to your array. +# +# In this tutorial, I don't preprocess the data. If your fNIRS device outputs +# data in the form of photon fluences at detectors, you can compute HbO and HbR +# concentrations using the modified Lambert-Beer model. +# +# In the following example we have data from 20 channels plus one events +# channel and we end up # with 20 HbO time series, 20 HbR time series and one +# events time series. + +############################################################################### +# Often fNIRS data contain time series for oxyhemoglobin (HbO) and +# deoxyhemoglobin (HbR) concentrations and you want to process them both. Using +# the *channel types* field in a mne.Info class instance, you can assign either +#'hbo' or 'hbr' channel type to your channel. +nChannels = 20 # number of physical channels +sampling_frequency = 15.625 +channel_names_fnirs = ['HbO '+"%.2d" % i for i in range(1,nChannels+1)] + + ['HbR '+"%.2d" % i for i in range(1,nChannels+1)] +channel_names = channel_names_fnirs + ['Events'] +channel_types = ['hbo' for i in range(nChannels)] + + ['hbr' for i in range(nChannels)] + ['stim'] +info = mne.create_info(ch_names=channel_names, sfreq=sampling_frequency, + ch_types=channel_types, montage=None) +info['lowpass'] = 0.0 +info['highpass'] = 0.0 + +############################################################################### +# Import fNIRS data from numpy array using mne.io.RawArray. +raw_data = mne.io.RawArray(data=raw_ndarray, info=info, first_samp=0, + verbose=None) + +############################################################################### +# In our data array the last row contains event data, coded as different +# in a time series synchronous with the fNIRS recording. Let's find such events +# using mne.find_events and plot the events timeline with mne.viz.plot_events. +events = mne.find_events(raw_data, stim_channel='Events', shortest_event=1) +event_id = {'ISI':1, 'Condition1': 2, 'Condition2': 4, 'Condition3': 8, + 'Pause':16} +color = {1:'black', 2:'green', 4:'magenta', 8:'cyan', 16:'yellow'} +mne.viz.plot_events(events, raw_data.info['sfreq'], raw_data.first_samp, + color=color, event_id=event_id) + +############################################################################### +# Sometimes we want to filter out unwanted frequencies from time series. For +# example, here we bandpass filter our data between 0.01 and 2 Hz. +l_freq = 0.01 # high-pass filter cutoff ( __/¯¯¯ ) +h_freq = 2 # low-pass filter cutoff ( ¯¯¯\__ ) +raw_data.filter(l_freq, h_freq) + +############################################################################### +# In order to have a general glimpse of the dataset, I start by plotting all +# channels into a single plot window. I choose to color-code HbO, HbR and +# events data using red, blue and black respectively. I also mark color-coded +# events positions, using the same colors previously used in +# mne.viz.plot_events. +scalings = dict(hbo=10e-6, hbr=10e-6, stim=1) +fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' + + str(h_freq) + ' Hz]' +plot_colors = dict(hbo='r', hbr='b', stim='k') +raw_data.plot(title=fig_title, events=events, start=0.0, color=plot_colors, + event_color=color, duration=np.max(raw_data.times), scalings=scalings, + order='original', n_channels=len(channel_names), remove_dc=False, + highpass=None, lowpass=None) + +############################################################################### +# I can later dive into the dataset and look at the raw time series more +# closely. Same color-codes through the plots is usually a good idea. +scalings = dict(hbo=10e-6, hbr=10e-6, stim=1) +fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' + + str(h_freq) + ' Hz]' +plot_colors = dict(hbo='r', hbr='b', stim='k') +fig = raw_data.plot(title=fig_title, events=events, start=0.0, + color=plot_colors, event_color=color, scalings=scalings, order='original', + duration=36, remove_dc=False, highpass=None, lowpass=None) + +############################################################################### +# We can now select epochs (defined as windows [-2,10] s around events) and +# drop unwanted ones (in this case I drop epochs with peak-to-peak amplitude +# larger than 1.5e-5 uM). Finally, we display the epochs time series. +tmin = -2 +tmax = 10 +reject = dict(hbo=1.5e-5, hbr=1.5e-5) +epochs = mne.Epochs(raw_data, events, event_id, tmin, tmax, proj=True, + baseline=(None, 0), preload=True, reject=reject) +fig_title = 'fNIRS Epochs' +epochs.plot(title=fig_title) + +############################################################################### +# Finally, we calculate and display the evoked signals, for example the one +# corresponding to the first condition in the paradigm. +evoked_1 = epochs['Condition1'].average() +evoked_1.plot() From e480fe82084d4dcb5e37e2d6598cfb4f1ca0ac4d Mon Sep 17 00:00:00 2001 From: Matteo Caffini Date: Thu, 11 May 2017 17:59:04 +0200 Subject: [PATCH 2/5] integrate syntax comments --- tutorials/Working_with_fNIRS_data.py | 122 +++++++++++++-------------- 1 file changed, 61 insertions(+), 61 deletions(-) diff --git a/tutorials/Working_with_fNIRS_data.py b/tutorials/Working_with_fNIRS_data.py index e2851368d4e..30ee7d2512e 100644 --- a/tutorials/Working_with_fNIRS_data.py +++ b/tutorials/Working_with_fNIRS_data.py @@ -1,124 +1,124 @@ """ .. _tut_fnirs_data +======================= Use MNE with fNIRS data ======================= -# Author: Matteo Caffini -# Contact: +MNE Python was implemented with EEG and MEG in mind, but can be useful to +process fNIRS data as #well. A common fNIRS data set (similarly to EEG and +MEG data) contains data recorded over time(a.k.a. samples) at many locations +on the scalp (a.k.a. channels). +Usually in fNIRS one looks for brain responses to external sensorial +stimulations (visual, acoustic, etc.). In order to do so, block design +studies are usually employed and data are later epoched to get evoked +responses. """ +# Author: Matteo Caffini (matteo.caffini@unitn.it) +# +# License: BSD (3-clause) + import numpy as np import mne raw_ndarray = np.load('/home/matteo/raw_ndarray.npy') -############################################################################### -# MNE Python was implemented with EEG and MEG in mind, but can be useful to -# process fNIRS data as #well. A common fNIRS data set (similarly to EEG and -# MEG data) contains data recorded over time(a.k.a. samples) at many locations -# on the scalp (a.k.a. channels). -# Usually in fNIRS one looks for brain responses to external sensorial -# stimulations (visual, acoustic, etc.). In order to do so, block design -# studies are usually employed and data are later epoched to get evoked -# responses. - ############################################################################### # fNIRS raw data come in various formats and the best way to dive into MNE is # to deal with your own format first and fit your data into a 2D numpy array # [channels x samples]. Once you prepared your time series for HbO and HbR -# concentrations, you can create a RawArray instance from the data array using -# mne.io.RawArray. If you have a time series of events synchronous with your +# concentrations, you can create a :class:`mne.io.RawArray` instance from the +# data array. If you have a time series of events synchronous with your # data you can append this as a last row to your array. -# +# # In this tutorial, I don't preprocess the data. If your fNIRS device outputs # data in the form of photon fluences at detectors, you can compute HbO and HbR # concentrations using the modified Lambert-Beer model. -# +# # In the following example we have data from 20 channels plus one events # channel and we end up # with 20 HbO time series, 20 HbR time series and one -# events time series. +# events time series. The provided data come from a dataset of neural +# correlates of visual perception of biological motion in newborns. ############################################################################### # Often fNIRS data contain time series for oxyhemoglobin (HbO) and -# deoxyhemoglobin (HbR) concentrations and you want to process them both. Using -# the *channel types* field in a mne.Info class instance, you can assign either -#'hbo' or 'hbr' channel type to your channel. -nChannels = 20 # number of physical channels -sampling_frequency = 15.625 -channel_names_fnirs = ['HbO '+"%.2d" % i for i in range(1,nChannels+1)] + - ['HbR '+"%.2d" % i for i in range(1,nChannels+1)] +# deoxyhemoglobin (HbR) concentrations and you want to process them both. +# Using the *channel types* field in a :class:`mne.Info` class instance, you +# can assign either 'hbo' or 'hbr' channel type to your channel. +n_channels = 20 # number of physical channels +sfreq = 15.625 +channel_names_fnirs = ['HbO %.2d' % i for i in range(1, n_channels + 1)] + + ['HbR %.2d' % i for i in range(1, n_channels + 1)] channel_names = channel_names_fnirs + ['Events'] -channel_types = ['hbo' for i in range(nChannels)] + - ['hbr' for i in range(nChannels)] + ['stim'] -info = mne.create_info(ch_names=channel_names, sfreq=sampling_frequency, - ch_types=channel_types, montage=None) -info['lowpass'] = 0.0 -info['highpass'] = 0.0 +channel_types = ['hbo'] * n_channels + ['hbr'] * n_channels + ['stim'] +info = mne.create_info(ch_names=channel_names, sfreq=sfreq, + ch_types=channel_types, montage=None) ############################################################################### -# Import fNIRS data from numpy array using mne.io.RawArray. +# Import fNIRS data from numpy array using :class:`mne.io.RawArray`. raw_data = mne.io.RawArray(data=raw_ndarray, info=info, first_samp=0, - verbose=None) + verbose=None) ############################################################################### # In our data array the last row contains event data, coded as different -# in a time series synchronous with the fNIRS recording. Let's find such events -# using mne.find_events and plot the events timeline with mne.viz.plot_events. +# numbers in a time series synchronous with the fNIRS recording. Let's find +# such events using :func:`mne.find_events` and plot the events timeline using +# :func:`mne.viz.plot_events`. events = mne.find_events(raw_data, stim_channel='Events', shortest_event=1) -event_id = {'ISI':1, 'Condition1': 2, 'Condition2': 4, 'Condition3': 8, - 'Pause':16} -color = {1:'black', 2:'green', 4:'magenta', 8:'cyan', 16:'yellow'} +event_id = {'ISI': 1, 'Biological': 2, 'Rotation': 4, 'Scrambled': 8, + 'Pause': 16} +color = {1: 'black', 2: 'green', 4: 'magenta', 8: 'cyan', 16: 'yellow'} mne.viz.plot_events(events, raw_data.info['sfreq'], raw_data.first_samp, - color=color, event_id=event_id) + color=color, event_id=event_id) ############################################################################### # Sometimes we want to filter out unwanted frequencies from time series. For # example, here we bandpass filter our data between 0.01 and 2 Hz. -l_freq = 0.01 # high-pass filter cutoff ( __/¯¯¯ ) -h_freq = 2 # low-pass filter cutoff ( ¯¯¯\__ ) -raw_data.filter(l_freq, h_freq) +l_freq = 0.01 # high-pass filter cutoff ( __/¯¯¯ ) +h_freq = 2 # low-pass filter cutoff ( ¯¯¯\__ ) +raw_data.filter(l_freq, h_freq, fir_design='firwin') ############################################################################### -# In order to have a general glimpse of the dataset, I start by plotting all -# channels into a single plot window. I choose to color-code HbO, HbR and -# events data using red, blue and black respectively. I also mark color-coded +# In order to have a general glimpse of the dataset, we start by plotting all +# channels into a single plot window. We choose to color-code HbO, HbR and +# events data using red, blue and black respectively. We also add color-coded # events positions, using the same colors previously used in -# mne.viz.plot_events. +# :func:`mne.viz.plot_events`. scalings = dict(hbo=10e-6, hbr=10e-6, stim=1) -fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' - + str(h_freq) + ' Hz]' +fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' + + str(h_freq) + ' Hz]' plot_colors = dict(hbo='r', hbr='b', stim='k') raw_data.plot(title=fig_title, events=events, start=0.0, color=plot_colors, - event_color=color, duration=np.max(raw_data.times), scalings=scalings, - order='original', n_channels=len(channel_names), remove_dc=False, - highpass=None, lowpass=None) + event_color=color, duration=np.max(raw_data.times), + scalings=scalings, order='original', + n_channels=len(channel_names), remove_dc=False, highpass=None, + lowpass=None) ############################################################################### -# I can later dive into the dataset and look at the raw time series more +# We can later dive into the dataset and look at the raw time series more # closely. Same color-codes through the plots is usually a good idea. -scalings = dict(hbo=10e-6, hbr=10e-6, stim=1) -fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' - + str(h_freq) + ' Hz]' +fig_title = 'fNIRS raw bandpass filtered %s Hz - %s Hz' % (l_freq, h_freq) plot_colors = dict(hbo='r', hbr='b', stim='k') fig = raw_data.plot(title=fig_title, events=events, start=0.0, - color=plot_colors, event_color=color, scalings=scalings, order='original', - duration=36, remove_dc=False, highpass=None, lowpass=None) + color=plot_colors, event_color=color, scalings=scalings, + order='original', duration=36, remove_dc=False, + highpass=None, lowpass=None) ############################################################################### # We can now select epochs (defined as windows [-2,10] s around events) and -# drop unwanted ones (in this case I drop epochs with peak-to-peak amplitude +# drop unwanted ones (in this case we drop epochs with peak-to-peak amplitude # larger than 1.5e-5 uM). Finally, we display the epochs time series. tmin = -2 tmax = 10 reject = dict(hbo=1.5e-5, hbr=1.5e-5) epochs = mne.Epochs(raw_data, events, event_id, tmin, tmax, proj=True, - baseline=(None, 0), preload=True, reject=reject) + baseline=(None, 0), preload=True, reject=reject) fig_title = 'fNIRS Epochs' epochs.plot(title=fig_title) ############################################################################### -# Finally, we calculate and display the evoked signals, for example the one -# corresponding to the first condition in the paradigm. -evoked_1 = epochs['Condition1'].average() +# After getting the epochs, we calculate and display the evoked signals, for +# example the one corresponding to the first condition in the paradigm. +evoked_1 = epochs['Biological'].average() evoked_1.plot() From fcae1c5cec71d20e661fa8ad9f33dceaa3f7755d Mon Sep 17 00:00:00 2001 From: Matteo Caffini Date: Fri, 12 May 2017 10:57:58 +0200 Subject: [PATCH 3/5] mini fix syntax --- tutorials/Working_with_fNIRS_data.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tutorials/Working_with_fNIRS_data.py b/tutorials/Working_with_fNIRS_data.py index 30ee7d2512e..3e0982d277d 100644 --- a/tutorials/Working_with_fNIRS_data.py +++ b/tutorials/Working_with_fNIRS_data.py @@ -86,8 +86,7 @@ # events positions, using the same colors previously used in # :func:`mne.viz.plot_events`. scalings = dict(hbo=10e-6, hbr=10e-6, stim=1) -fig_title = 'fNIRS Raw Bandpass filtered [' + str(l_freq) + ' Hz, ' + - str(h_freq) + ' Hz]' +fig_title = 'fNIRS raw bandpass filtered %s Hz - %s Hz' % (l_freq, h_freq) plot_colors = dict(hbo='r', hbr='b', stim='k') raw_data.plot(title=fig_title, events=events, start=0.0, color=plot_colors, event_color=color, duration=np.max(raw_data.times), From 938181e69bf1a563923e8317d13faa16ab94bd5c Mon Sep 17 00:00:00 2001 From: Matteo Caffini Date: Mon, 15 May 2017 15:42:19 +0200 Subject: [PATCH 4/5] fix filename and syntax --- tutorials/{Working_with_fNIRS_data.py => plot_fnirs_data.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tutorials/{Working_with_fNIRS_data.py => plot_fnirs_data.py} (100%) diff --git a/tutorials/Working_with_fNIRS_data.py b/tutorials/plot_fnirs_data.py similarity index 100% rename from tutorials/Working_with_fNIRS_data.py rename to tutorials/plot_fnirs_data.py From c2658487c5c452fff4de11fc18839d794615b145 Mon Sep 17 00:00:00 2001 From: Matteo Caffini Date: Mon, 15 May 2017 15:47:02 +0200 Subject: [PATCH 5/5] fix filename and syntax 2 --- tutorials/plot_fnirs_data.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tutorials/plot_fnirs_data.py b/tutorials/plot_fnirs_data.py index 3e0982d277d..067085262f7 100644 --- a/tutorials/plot_fnirs_data.py +++ b/tutorials/plot_fnirs_data.py @@ -5,10 +5,10 @@ Use MNE with fNIRS data ======================= -MNE Python was implemented with EEG and MEG in mind, but can be useful to -process fNIRS data as #well. A common fNIRS data set (similarly to EEG and -MEG data) contains data recorded over time(a.k.a. samples) at many locations -on the scalp (a.k.a. channels). +MNE Python was originally implemented with EEG and MEG in mind, but can be +useful to process fNIRS data as well. A common fNIRS data set (similarly to EEG +and MEG data) contains data recorded over time(a.k.a. samples) at many +locations on the scalp (a.k.a. channels). Usually in fNIRS one looks for brain responses to external sensorial stimulations (visual, acoustic, etc.). In order to do so, block design studies are usually employed and data are later epoched to get evoked