From 16cf86139dbdaa2007cee54d975343dec0c8c4da Mon Sep 17 00:00:00 2001 From: Giovanni Date: Tue, 19 Apr 2022 15:37:57 +0200 Subject: [PATCH 1/3] Added FOOOF section on 01r_n170_viz --- examples/visual_n170/01r__n170_viz.py | 85 +++++++++++++++++++++------ 1 file changed, 67 insertions(+), 18 deletions(-) diff --git a/examples/visual_n170/01r__n170_viz.py b/examples/visual_n170/01r__n170_viz.py index 9aa48f7d..b0850d73 100644 --- a/examples/visual_n170/01r__n170_viz.py +++ b/examples/visual_n170/01r__n170_viz.py @@ -23,17 +23,24 @@ # --------------------- # Some standard pythonic imports +from eegnb.analysis.utils import load_data, plot_conditions +from mne import Epochs, find_events +import matplotlib.pyplot as plt +from fooof.plts import plot_spectra +from fooof import FOOOF +from eegnb.datasets import fetch_dataset +from mne.time_frequency import psd_welch import os from collections import OrderedDict import warnings warnings.filterwarnings('ignore') # MNE functions -from mne import Epochs,find_events # EEG-Notebooks functions -from eegnb.analysis.utils import load_data,plot_conditions -from eegnb.datasets import fetch_dataset + +# FOOOF functions + # sphinx_gallery_thumbnail_number = 3 @@ -49,18 +56,19 @@ ################################################################################################### -eegnb_data_path = os.path.join(os.path.expanduser('~/'),'.eegnb', 'data') +eegnb_data_path = os.path.join(os.path.expanduser('~/'), '.eegnb', 'data') n170_data_path = os.path.join(eegnb_data_path, 'visual-N170', 'eegnb_examples') -# If dataset hasn't been downloaded yet, download it +# If dataset hasn't been downloaded yet, download it if not os.path.isdir(n170_data_path): - fetch_dataset(data_dir=eegnb_data_path, experiment='visual-N170', site='eegnb_examples'); + fetch_dataset(data_dir=eegnb_data_path, + experiment='visual-N170', site='eegnb_examples') subject = 1 session = 1 -raw = load_data(subject,session, +raw = load_data(subject, session, experiment='visual-N170', site='eegnb_examples', device_name='muse2016_bfn', - data_dir = eegnb_data_path) + data_dir=eegnb_data_path) ################################################################################################### # Visualize the power spectrum @@ -72,8 +80,8 @@ # Filtering # ---------------------------- -raw.filter(1,30, method='iir') -raw.plot_psd(fmin=1, fmax=30); +raw.filter(1, 30, method='iir') +raw.plot_psd(fmin=1, fmax=30) ################################################################################################### # Epoching @@ -84,10 +92,10 @@ event_id = {'House': 1, 'Face': 2} # Create an MNE Epochs object representing all the epochs around stimulus presentation -epochs = Epochs(raw, events=events, event_id=event_id, +epochs = Epochs(raw, events=events, event_id=event_id, tmin=-0.1, tmax=0.6, baseline=None, - reject={'eeg': 5e-5}, preload=True, - verbose=False, picks=[0,1,2,3]) + reject={'eeg': 5e-5}, preload=True, + verbose=False, picks=[0, 1, 2, 3]) print('sample drop %: ', (1 - len(epochs.events)/len(events)) * 100) epochs @@ -99,13 +107,54 @@ conditions['House'] = [1] conditions['Face'] = [2] -fig, ax = plot_conditions(epochs, conditions=conditions, +fig, ax = plot_conditions(epochs, conditions=conditions, ci=97.5, n_boot=1000, title='', - diff_waveform=None, #(1, 2)) - channel_order=[1,0,2,3]) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes + diff_waveform=None, # (1, 2)) + channel_order=[1, 0, 2, 3]) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes # Manually adjust the ylims -for i in [0,2]: ax[i].set_ylim([-0.5,0.5]) -for i in [1,3]: ax[i].set_ylim([-1.5,2.5]) +for i in [0, 2]: + ax[i].set_ylim([-0.5, 0.5]) +for i in [1, 3]: + ax[i].set_ylim([-1.5, 2.5]) + +################################################################################################### +# Frequency Oscillations & One Over F (FOOOF) +# ---------------------------- + +# Averaging epochs by event type (House, Face) +av_epochs = epochs.average(picks=[0, 1, 2, 3], by_event_type=True) + +fms = list() +labels = ["House", "Face"] + +# Computing PSD for each event type and evaluating FOOOF +for i, e in enumerate(av_epochs): + power, freq = psd_welch(e, n_fft=len(e.times), n_overlap=len(e.times)/2, + proj=False, average='mean', window='hamming') + + fms.append(FOOOF(aperiodic_mode='knee', verbose=False)) + fms[i].fit(freq, power[3], freq_range=[1, 30]) + + fig, ax = plt.subplots(figsize=[15, 5]) + fms[i].plot(plot_peaks='line-shade', ax=ax) + plt.title(f"FOOOF - {labels[i]}") + +# Visualization of spectra parametrization +fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots( + 2, 2, sharex=True, figsize=[15, 10]) +plot_spectra(fms[0].freqs, [fms[0].power_spectrum, + fms[1].power_spectrum], ax=ax1, labels=labels) +ax1.set_title("Original Spectra") + +plot_spectra(fms[0].freqs, [fms[0]._spectrum_flat, + fms[1]._spectrum_flat], ax=ax2, labels=labels) +ax2.set_title("Flat Spectrum (NO AP)") +plot_spectra(fms[0].freqs, [fms[0]._peak_fit, + fms[1]._peak_fit], ax=ax3, labels=labels) +ax3.set_title("Peak Fit") +plot_spectra(fms[0].freqs, [fms[0]._ap_fit, fms[1]._ap_fit], + ax=ax4, labels=labels) +ax4.set_title("AP Fit") From 6b90e628505045f8aa385aa7b23076af598471a1 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 3 Oct 2022 19:35:01 +0200 Subject: [PATCH 2/3] Adding a new example file for FOOOF analysis. Reverting previous editing on first example. --- examples/visual_n170/01r__n170_viz.py | 118 ++++++++------------ examples/visual_n170/03r_n170_fooof.py | 148 +++++++++++++++++++++++++ 2 files changed, 193 insertions(+), 73 deletions(-) create mode 100644 examples/visual_n170/03r_n170_fooof.py diff --git a/examples/visual_n170/01r__n170_viz.py b/examples/visual_n170/01r__n170_viz.py index b0850d73..3c7c96a9 100644 --- a/examples/visual_n170/01r__n170_viz.py +++ b/examples/visual_n170/01r__n170_viz.py @@ -23,24 +23,18 @@ # --------------------- # Some standard pythonic imports -from eegnb.analysis.utils import load_data, plot_conditions -from mne import Epochs, find_events -import matplotlib.pyplot as plt -from fooof.plts import plot_spectra -from fooof import FOOOF -from eegnb.datasets import fetch_dataset -from mne.time_frequency import psd_welch import os from collections import OrderedDict import warnings -warnings.filterwarnings('ignore') + +warnings.filterwarnings("ignore") # MNE functions +from mne import Epochs, find_events # EEG-Notebooks functions - -# FOOOF functions - +from eegnb.analysis.utils import load_data, plot_conditions +from eegnb.datasets import fetch_dataset # sphinx_gallery_thumbnail_number = 3 @@ -56,19 +50,25 @@ ################################################################################################### -eegnb_data_path = os.path.join(os.path.expanduser('~/'), '.eegnb', 'data') -n170_data_path = os.path.join(eegnb_data_path, 'visual-N170', 'eegnb_examples') +eegnb_data_path = os.path.join(os.path.expanduser("~/"), ".eegnb", "data") +n170_data_path = os.path.join(eegnb_data_path, "visual-N170", "eegnb_examples") # If dataset hasn't been downloaded yet, download it if not os.path.isdir(n170_data_path): - fetch_dataset(data_dir=eegnb_data_path, - experiment='visual-N170', site='eegnb_examples') + fetch_dataset( + data_dir=eegnb_data_path, experiment="visual-N170", site="eegnb_examples" + ) subject = 1 session = 1 -raw = load_data(subject, session, - experiment='visual-N170', site='eegnb_examples', device_name='muse2016_bfn', - data_dir=eegnb_data_path) +raw = load_data( + subject, + session, + experiment="visual-N170", + site="eegnb_examples", + device_name="muse2016_bfn", + data_dir=eegnb_data_path, +) ################################################################################################### # Visualize the power spectrum @@ -80,7 +80,7 @@ # Filtering # ---------------------------- -raw.filter(1, 30, method='iir') +raw.filter(1, 30, method="iir") raw.plot_psd(fmin=1, fmax=30) ################################################################################################### @@ -89,14 +89,22 @@ # Create an array containing the timestamps and type of each stimulus (i.e. face or house) events = find_events(raw) -event_id = {'House': 1, 'Face': 2} +event_id = {"House": 1, "Face": 2} # Create an MNE Epochs object representing all the epochs around stimulus presentation -epochs = Epochs(raw, events=events, event_id=event_id, - tmin=-0.1, tmax=0.6, baseline=None, - reject={'eeg': 5e-5}, preload=True, - verbose=False, picks=[0, 1, 2, 3]) -print('sample drop %: ', (1 - len(epochs.events)/len(events)) * 100) +epochs = Epochs( + raw, + events=events, + event_id=event_id, + tmin=-0.1, + tmax=0.6, + baseline=None, + reject={"eeg": 5e-5}, + preload=True, + verbose=False, + picks=[0, 1, 2, 3], +) +print("sample drop %: ", (1 - len(epochs.events) / len(events)) * 100) epochs ################################################################################################### @@ -104,57 +112,21 @@ # ---------------------------- conditions = OrderedDict() -conditions['House'] = [1] -conditions['Face'] = [2] - -fig, ax = plot_conditions(epochs, conditions=conditions, - ci=97.5, n_boot=1000, title='', - diff_waveform=None, # (1, 2)) - channel_order=[1, 0, 2, 3]) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes +conditions["House"] = [1] +conditions["Face"] = [2] + +fig, ax = plot_conditions( + epochs, + conditions=conditions, + ci=97.5, + n_boot=1000, + title="", + diff_waveform=None, # (1, 2)) + channel_order=[1, 0, 2, 3], +) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes # Manually adjust the ylims for i in [0, 2]: ax[i].set_ylim([-0.5, 0.5]) for i in [1, 3]: ax[i].set_ylim([-1.5, 2.5]) - -################################################################################################### -# Frequency Oscillations & One Over F (FOOOF) -# ---------------------------- - -# Averaging epochs by event type (House, Face) -av_epochs = epochs.average(picks=[0, 1, 2, 3], by_event_type=True) - -fms = list() -labels = ["House", "Face"] - -# Computing PSD for each event type and evaluating FOOOF -for i, e in enumerate(av_epochs): - power, freq = psd_welch(e, n_fft=len(e.times), n_overlap=len(e.times)/2, - proj=False, average='mean', window='hamming') - - fms.append(FOOOF(aperiodic_mode='knee', verbose=False)) - fms[i].fit(freq, power[3], freq_range=[1, 30]) - - fig, ax = plt.subplots(figsize=[15, 5]) - fms[i].plot(plot_peaks='line-shade', ax=ax) - plt.title(f"FOOOF - {labels[i]}") - -# Visualization of spectra parametrization -fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots( - 2, 2, sharex=True, figsize=[15, 10]) -plot_spectra(fms[0].freqs, [fms[0].power_spectrum, - fms[1].power_spectrum], ax=ax1, labels=labels) -ax1.set_title("Original Spectra") - -plot_spectra(fms[0].freqs, [fms[0]._spectrum_flat, - fms[1]._spectrum_flat], ax=ax2, labels=labels) -ax2.set_title("Flat Spectrum (NO AP)") - -plot_spectra(fms[0].freqs, [fms[0]._peak_fit, - fms[1]._peak_fit], ax=ax3, labels=labels) -ax3.set_title("Peak Fit") - -plot_spectra(fms[0].freqs, [fms[0]._ap_fit, fms[1]._ap_fit], - ax=ax4, labels=labels) -ax4.set_title("AP Fit") diff --git a/examples/visual_n170/03r_n170_fooof.py b/examples/visual_n170/03r_n170_fooof.py new file mode 100644 index 00000000..9cde85bf --- /dev/null +++ b/examples/visual_n170/03r_n170_fooof.py @@ -0,0 +1,148 @@ +""" +N170 FOOOF Analysis +=============================== + +This example execute a Frequency Over One Over F (FOOOF) analysis +on the N170 faces/houses dataset, and compares it with respect to a classical +power spectra analysis. + +The data used is exactly the same as in the N170 `load_and_visualize` example. + +""" + +################################################################################################### +# Setup +# --------------------- + +# Some standard pythonic imports +import os + +# EEG-Notebooks functions +from eegnb.analysis.utils import load_data, plot_conditions +from eegnb.datasets import fetch_dataset + +# MNE functions +from mne import Epochs, find_events +from mne.time_frequency import psd_welch + +# FOOOF functions +from fooof import FOOOF +from fooof.plts import plot_spectra + +# Matplotlib for visualization +import matplotlib.pyplot as plt + +################################################################################################### +# Load Data +# --------------------- +# +# ( See the n170 `load_and_visualize` example for further description of this) +# + +eegnb_data_path = os.path.join(os.path.expanduser("~/"), ".eegnb", "data") +n170_data_path = os.path.join(eegnb_data_path, "visual-N170", "eegnb_examples") + +# If dataset hasn't been downloaded yet, download it +if not os.path.isdir(n170_data_path): + fetch_dataset( + data_dir=eegnb_data_path, experiment="visual-N170", site="eegnb_examples" + ) + +subject = 1 +session = 1 +raw = load_data( + subject, + session, + experiment="visual-N170", + site="eegnb_examples", + device_name="muse2016_bfn", + data_dir=eegnb_data_path, +) + +################################################################################################### +# Filtering +# ---------------------------- + +raw.filter(1, 30, method="iir") +raw.plot_psd(fmin=1, fmax=30) + +################################################################################################### +# Epoching +# ---------------------------- + +# Create an array containing the timestamps and type of each stimulus (i.e. face or house) +events = find_events(raw) +event_id = {"House": 1, "Face": 2} + +# Create an MNE Epochs object representing all the epochs around stimulus presentation +epochs = Epochs( + raw, + events=events, + event_id=event_id, + tmin=-0.1, + tmax=0.6, + baseline=None, + reject={"eeg": 5e-5}, + preload=True, + verbose=False, + picks=[0, 1, 2, 3], +) +print("sample drop %: ", (1 - len(epochs.events) / len(events)) * 100) +epochs + +################################################################################################### + +# Averaging epochs by event type (House, Face) +av_epochs = epochs.average(picks=[0, 1, 2, 3], by_event_type=True) + +fms = list() +labels = ["House", "Face"] + +# Computing PSD for each event type and evaluating FOOOF +for i, e in enumerate(av_epochs): + power, freq = psd_welch( + e, + n_fft=len(e.times), + n_overlap=len(e.times) / 2, + proj=False, + average="mean", + window="hamming", + ) + + # A full exponential equation to fit the aperiodic component is chosen, thus an aperiodic "knee" mode is requested + fms.append(FOOOF(aperiodic_mode="knee", verbose=False)) + fms[i].fit(freq, power[3], freq_range=[1, 30]) + + # Plotting the results of FOOOF fitting spectra in the range 1-30 Hz + fig, ax = plt.subplots(figsize=[15, 5]) + fms[i].plot(plot_peaks="line-shade", ax=ax) + plt.title(f"FOOOF - {labels[i]}") + plt.show(block=False) + +################################################################################################### +# Results Visualization +# ---------------------------- +# +# Interestingly, it is easily possible to see in the Peak Fit subplot how "Face" subgroup presents actual lower power with +# respect to "Face" one althought in the Original Spectra it could be assessed the opposite. This can be justified by the +# higher aperiodic contribute of Face subgroup that does not reflect however the oscillation activity power. +# +# + +fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=[15, 10]) +plot_spectra( + fms[0].freqs, [fms[0].power_spectrum, fms[1].power_spectrum], ax=ax1, labels=labels +) +ax1.set_title("Original Spectra") + +plot_spectra( + fms[0].freqs, [fms[0]._spectrum_flat, fms[1]._spectrum_flat], ax=ax2, labels=labels +) +ax2.set_title("Flat Spectrum (NO AP)") + +plot_spectra(fms[0].freqs, [fms[0]._peak_fit, fms[1]._peak_fit], ax=ax3, labels=labels) +ax3.set_title("Peak Fit") + +plot_spectra(fms[0].freqs, [fms[0]._ap_fit, fms[1]._ap_fit], ax=ax4, labels=labels) +ax4.set_title("AP Fit") +plt.show() From 4e0507daf3a55a22638a17a8e9d65b8edf2f9bca Mon Sep 17 00:00:00 2001 From: Giovanni Date: Tue, 4 Oct 2022 14:38:08 +0200 Subject: [PATCH 3/3] reverting original format to 01r_n179_viz.py --- examples/visual_n170/01r__n170_viz.py | 79 ++++++++++----------------- 1 file changed, 29 insertions(+), 50 deletions(-) diff --git a/examples/visual_n170/01r__n170_viz.py b/examples/visual_n170/01r__n170_viz.py index 3c7c96a9..9aa48f7d 100644 --- a/examples/visual_n170/01r__n170_viz.py +++ b/examples/visual_n170/01r__n170_viz.py @@ -26,14 +26,13 @@ import os from collections import OrderedDict import warnings - -warnings.filterwarnings("ignore") +warnings.filterwarnings('ignore') # MNE functions -from mne import Epochs, find_events +from mne import Epochs,find_events # EEG-Notebooks functions -from eegnb.analysis.utils import load_data, plot_conditions +from eegnb.analysis.utils import load_data,plot_conditions from eegnb.datasets import fetch_dataset # sphinx_gallery_thumbnail_number = 3 @@ -50,25 +49,18 @@ ################################################################################################### -eegnb_data_path = os.path.join(os.path.expanduser("~/"), ".eegnb", "data") -n170_data_path = os.path.join(eegnb_data_path, "visual-N170", "eegnb_examples") +eegnb_data_path = os.path.join(os.path.expanduser('~/'),'.eegnb', 'data') +n170_data_path = os.path.join(eegnb_data_path, 'visual-N170', 'eegnb_examples') -# If dataset hasn't been downloaded yet, download it +# If dataset hasn't been downloaded yet, download it if not os.path.isdir(n170_data_path): - fetch_dataset( - data_dir=eegnb_data_path, experiment="visual-N170", site="eegnb_examples" - ) + fetch_dataset(data_dir=eegnb_data_path, experiment='visual-N170', site='eegnb_examples'); subject = 1 session = 1 -raw = load_data( - subject, - session, - experiment="visual-N170", - site="eegnb_examples", - device_name="muse2016_bfn", - data_dir=eegnb_data_path, -) +raw = load_data(subject,session, + experiment='visual-N170', site='eegnb_examples', device_name='muse2016_bfn', + data_dir = eegnb_data_path) ################################################################################################### # Visualize the power spectrum @@ -80,8 +72,8 @@ # Filtering # ---------------------------- -raw.filter(1, 30, method="iir") -raw.plot_psd(fmin=1, fmax=30) +raw.filter(1,30, method='iir') +raw.plot_psd(fmin=1, fmax=30); ################################################################################################### # Epoching @@ -89,22 +81,14 @@ # Create an array containing the timestamps and type of each stimulus (i.e. face or house) events = find_events(raw) -event_id = {"House": 1, "Face": 2} +event_id = {'House': 1, 'Face': 2} # Create an MNE Epochs object representing all the epochs around stimulus presentation -epochs = Epochs( - raw, - events=events, - event_id=event_id, - tmin=-0.1, - tmax=0.6, - baseline=None, - reject={"eeg": 5e-5}, - preload=True, - verbose=False, - picks=[0, 1, 2, 3], -) -print("sample drop %: ", (1 - len(epochs.events) / len(events)) * 100) +epochs = Epochs(raw, events=events, event_id=event_id, + tmin=-0.1, tmax=0.6, baseline=None, + reject={'eeg': 5e-5}, preload=True, + verbose=False, picks=[0,1,2,3]) +print('sample drop %: ', (1 - len(epochs.events)/len(events)) * 100) epochs ################################################################################################### @@ -112,21 +96,16 @@ # ---------------------------- conditions = OrderedDict() -conditions["House"] = [1] -conditions["Face"] = [2] - -fig, ax = plot_conditions( - epochs, - conditions=conditions, - ci=97.5, - n_boot=1000, - title="", - diff_waveform=None, # (1, 2)) - channel_order=[1, 0, 2, 3], -) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes +conditions['House'] = [1] +conditions['Face'] = [2] + +fig, ax = plot_conditions(epochs, conditions=conditions, + ci=97.5, n_boot=1000, title='', + diff_waveform=None, #(1, 2)) + channel_order=[1,0,2,3]) # reordering of epochs.ch_names according to [[0,2],[1,3]] of subplot axes # Manually adjust the ylims -for i in [0, 2]: - ax[i].set_ylim([-0.5, 0.5]) -for i in [1, 3]: - ax[i].set_ylim([-1.5, 2.5]) +for i in [0,2]: ax[i].set_ylim([-0.5,0.5]) +for i in [1,3]: ax[i].set_ylim([-1.5,2.5]) + +