From a255fcf2fe6dbf63ca73fe86f570ec2be9f496a3 Mon Sep 17 00:00:00 2001 From: Luyanda Mdanda Date: Sun, 14 Apr 2019 00:34:58 -0700 Subject: [PATCH 1/5] 1st draft of mne example --- examples/plot_mne_example.py | 223 +++++++++++++++++++++++++++++++++++ 1 file changed, 223 insertions(+) create mode 100644 examples/plot_mne_example.py diff --git a/examples/plot_mne_example.py b/examples/plot_mne_example.py new file mode 100644 index 000000000..9ed6feeb8 --- /dev/null +++ b/examples/plot_mne_example.py @@ -0,0 +1,223 @@ +""" +Using fooof with MNE +================= + +This examples illustrates how to use fooof with `MNE +`_ and create topographical plots + +This tutorial does require that you have MNE installed. If you don't already have +MNE, you can follow instructions to get it `here +`_. +""" + +################################################################################################### + + +#General imports +import numpy as np +import pandas as pd +from matplotlib import cm, colors, colorbar +import matplotlib.pyplot as plt + +# Import MNE, as well as the MNE sample dataset +import mne +from mne import io +from mne.datasets import sample +from mne.viz import plot_topomap + +# FOOOF imports +from fooof import FOOOF, FOOOFGroup +from fooof.analysis import * + + +################################################################################################### +# Load & Check MNE Data +# --------------------- +# +# First, we will load the example dataset from MNE, and have a quick look at the data. +# +# The MNE sample dataset is a combined MEG/EEG recording with an audiovisual task, as +# described `here `_. +# +# For the current example, we are going to sub-select only the EEG data, +# and analyze it as continuous (non-epoched) data. +# + +################################################################################################### + +# Get the data path for the MNE example data +raw_fname = sample.data_path() + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' +event_fname = sample.data_path() + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' + +# Load the file of example MNE data +raw = mne.io.read_raw_fif(raw_fname, preload=True) + +################################################################################################### + +# Select EEG channels from the dataset +raw = raw.pick_types(meg=False, eeg=True, eog=False, exclude='bads') + +################################################################################################### + +# Explicitly adding a default EEG average reference required for source localization +raw.set_eeg_reference() + +################################################################################################### + +# Explore the data +raw.info() + +################################################################################################### + +# Plot the channel locations and labels +raw.plot_sensors(show_names=True) + +################################################################################################### + +# Creating epochs +reject = dict(eeg=180e-6) +event_id, tmin, tmax = {'left/auditory': 1}, -0.2, 0.5 +events = mne.read_events(event_fname) +epochs = mne.Epochs(raw, events=events, event_id=event_id, + tmin=5, tmax=125, baseline=None, preload=True) + +################################################################################################### + +# Creating Power Spectra Densities +psd, freqs = mne.time_frequency.psd_welch(epochs, fmin=1., fmax=50., n_fft=2000, n_overlap=250, n_per_seg=500) + +################################################################################################### +# fooofgroup +# ------------------ +# +# The FOOOFGroup object is used to fit +# FOOOF models across the power spectra. A list of FOOOFGroup objects is returned. +# + +################################################################################################### + +# Initialize a FOOOF group object, with desired settings +fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_amplitude=0.075, + max_n_peaks=6, peak_threshold=1, verbose=False) + +################################################################################################### + +# Shaping the PSDs so that we can run fooofgroup +fooof_psd = np.squeeze(psd[4,:,:]) +n_channels, n_freq = fooof_psd.shape +num_blocks = len(mne.read_events(event_fname)) + +# Setting frequency range +freq_range = [3, 32] + +# This returns a list of FOOOFGRoup objects +fg.fit(freqs, fooof_psd, freq_range) +fg.plot() + +################################################################################################### + +# Define bands of interest +bands = {'theta': [2,7], + 'alpha': [8,14], + 'beta': [15,30]} + +# Creating dictionaries to store all the aperiodic properties across frequencies +results = {} +for band_name in bands.keys(): + results[band_name] = np.zeros(shape=[ num_blocks, n_channels, n_feats]) + +# Creating dictionaries to store all the periodic properties across frequencies +slope_results = np.zeros(shape=[num_blocks, n_channels, 2]) + +################################################################################################### + +# Populating periodic and aperiodic values +for block in range(0, num_blocks): + for ind, res in enumerate(fg): + slope_results[block, ind, :] = res.background_params + for band_label, band_range in bands.items(): + results[band_label][block, ind, :] = get_band_peak(res.peak_params, band_range, True) + +################################################################################################### +# Plotting topologies +# -------------- +# +# The following will illustrate the process of exploring the previoulsy generated +# oscillatory values + +################################################################################################### +# This is the function required to create color bars for topographies +def plot_topo_colorbar(vmin, vmax, label): + """ + Create a colorbar for the topography plots + vmin: int + vmax: int + label: str + """ + fig = plt.figure(figsize=(2, 3)) + ax1 = fig.add_axes([0.9, 0.25, 0.15, 0.9]) + + cmap = cm.viridis + norm = colors.Normalize(vmin=vmin, vmax=vmax) + + cb1 = colorbar.ColorbarBase(plt.gca(), cmap=cmap, + norm=norm, orientation='vertical')all_fg.plot() + +################################################################################################### + +# Settings: In this example, we will be plotting the Alpha Central Frequency +# and oscillatory slope +band = 'alpha' +cur_data = results[band] +feats = ["CFS", "AMPS", "BWS"] +topo_dat = np.mean(cur_data,0) + +################################################################################################### + +# Looking at the alpha Central Frequeuncy +print('CURRENT FEATURE:', feats[0]) +disp_dat = topo_dat[:,0] + +inds = np.where(np.isnan(disp_dat)) +disp_dat[inds] = np.nanmean(disp_dat) + +vbuffer = 0.1 * (disp_dat.max() - disp_dat.min()) +vmin, vmax, = disp_dat.min() - vbuffer, disp_dat.max() + vbuffer + +fig, ax = plt.subplots() +plot_topo_colorbar(vmin, vmax, feats[0]) +mne.viz.plot_topomap(disp_dat, raw.info, vmin=vmin, vmax=vmax, cmap=cm.viridis, contours=0, axes=ax) + +################################################################################################### + +# Looking at the Oscillatory +cur_data = slope_results + +topo_dat = np.mean(cur_data,0) + +print('CURRENT FEATURE:', slope_feat[1]) +disp_dat = topo_dat[:,1] + +inds = np.where(np.isnan(disp_dat)) +disp_dat[inds] = np.nanmean(disp_dat) + +vbuffer = 0.1 * (disp_dat.max() - disp_dat.min()) +vmin, vmax, = disp_dat.min() - vbuffer, disp_dat.max() + vbuffer + +fig, ax = plt.subplots() +plot_topo_colorbar(vmin, vmax, slope_feat[1]) +mne.viz.plot_topomap(disp_dat, raw.info, vmin=vmin, vmax=vmax, cmap=cm.viridis, contours=0, axes=ax) + + +################################################################################################### +# Alternative epoching methods +# -------------- +# +# Antother way to epoch your data based on time by avaeraging across events as suppose to +# selecting specifc time sections + +################################################################################################### + +epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, + reject=reject) +epochs = mne.Epochs(raw_no_ref, **epochs_params).average() \ No newline at end of file From c561026bfc6e2ea3efb5d3bfad33364c47cdcb03 Mon Sep 17 00:00:00 2001 From: Luyanda Mdanda Date: Thu, 18 Apr 2019 15:15:44 -0700 Subject: [PATCH 2/5] changes to header and aperiodic --- .gitignore | 1 + examples/plot_mne_example.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 59702341e..125bcc1bc 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ # Ignore Mac folder profile files *DS_Store* +.vs* # Ignore coverage file .coverage diff --git a/examples/plot_mne_example.py b/examples/plot_mne_example.py index 9ed6feeb8..90da71cf9 100644 --- a/examples/plot_mne_example.py +++ b/examples/plot_mne_example.py @@ -1,6 +1,6 @@ """ Using fooof with MNE -================= +======================= This examples illustrates how to use fooof with `MNE `_ and create topographical plots @@ -134,7 +134,7 @@ # Populating periodic and aperiodic values for block in range(0, num_blocks): for ind, res in enumerate(fg): - slope_results[block, ind, :] = res.background_params + slope_results[block, ind, :] = res.aperiodic_params for band_label, band_range in bands.items(): results[band_label][block, ind, :] = get_band_peak(res.peak_params, band_range, True) From acd84593c71bfb8a07760fc4c8b3bbc3eb84a4ba Mon Sep 17 00:00:00 2001 From: Luyanda Mdanda Date: Fri, 19 Apr 2019 13:22:57 -0700 Subject: [PATCH 3/5] 2nd draft of MNE example --- examples/plot_mne_example.py | 83 +++++++++++++++++------------------- 1 file changed, 39 insertions(+), 44 deletions(-) diff --git a/examples/plot_mne_example.py b/examples/plot_mne_example.py index 90da71cf9..e3888ba5b 100644 --- a/examples/plot_mne_example.py +++ b/examples/plot_mne_example.py @@ -1,11 +1,12 @@ """ Using fooof with MNE -======================= +==================== -This examples illustrates how to use fooof with `MNE -`_ and create topographical plots +This examples illustrates how to use fooof with MNE +and create topographical plots -This tutorial does require that you have MNE installed. If you don't already have +This tutorial does require that you have `MNE +`_ installed. If you don't already have MNE, you can follow instructions to get it `here `_. """ @@ -50,7 +51,7 @@ event_fname = sample.data_path() + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' # Load the file of example MNE data -raw = mne.io.read_raw_fif(raw_fname, preload=True) +raw = mne.io.read_raw_fif(raw_fname, preload=True, verbose=False) ################################################################################################### @@ -64,11 +65,6 @@ ################################################################################################### -# Explore the data -raw.info() - -################################################################################################### - # Plot the channel locations and labels raw.plot_sensors(show_names=True) @@ -76,19 +72,20 @@ # Creating epochs reject = dict(eeg=180e-6) -event_id, tmin, tmax = {'left/auditory': 1}, -0.2, 0.5 +event_id = {'left/auditory': 1} events = mne.read_events(event_fname) epochs = mne.Epochs(raw, events=events, event_id=event_id, - tmin=5, tmax=125, baseline=None, preload=True) + tmin=5, tmax=125, baseline=None, preload=True) ################################################################################################### # Creating Power Spectra Densities -psd, freqs = mne.time_frequency.psd_welch(epochs, fmin=1., fmax=50., n_fft=2000, n_overlap=250, n_per_seg=500) +psd, freqs = mne.time_frequency.psd_welch(epochs, fmin=1., fmax=50., n_fft=2000, + n_overlap=250, n_per_seg=500) ################################################################################################### # fooofgroup -# ------------------ +# ---------- # # The FOOOFGroup object is used to fit # FOOOF models across the power spectra. A list of FOOOFGroup objects is returned. @@ -97,18 +94,18 @@ ################################################################################################### # Initialize a FOOOF group object, with desired settings -fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_amplitude=0.075, - max_n_peaks=6, peak_threshold=1, verbose=False) +fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.075, + max_n_peaks=6, peak_threshold=1, verbose=False) ################################################################################################### # Shaping the PSDs so that we can run fooofgroup -fooof_psd = np.squeeze(psd[4,:,:]) +fooof_psd = np.squeeze(psd[0,:,:]) n_channels, n_freq = fooof_psd.shape num_blocks = len(mne.read_events(event_fname)) # Setting frequency range -freq_range = [3, 32] +freq_range = [3, 35] # This returns a list of FOOOFGRoup objects fg.fit(freqs, fooof_psd, freq_range) @@ -116,33 +113,30 @@ ################################################################################################### +# Periodic features +feats = ["CFS", "PWS", "BWS"] + +# Aperiodic features +aperiodic_feats = ["Offset","Exponent"] + # Define bands of interest -bands = {'theta': [2,7], - 'alpha': [8,14], +bands = {'theta': [3,7], + 'alpha': [7,14], 'beta': [15,30]} # Creating dictionaries to store all the aperiodic properties across frequencies results = {} for band_name in bands.keys(): - results[band_name] = np.zeros(shape=[ num_blocks, n_channels, n_feats]) + results[band_name] = np.zeros(shape=[ num_blocks, n_channels, len(feats)]) # Creating dictionaries to store all the periodic properties across frequencies -slope_results = np.zeros(shape=[num_blocks, n_channels, 2]) +exponent_results = np.zeros(shape=[num_blocks, n_channels, len(aperiodic_feats)]) ################################################################################################### - -# Populating periodic and aperiodic values -for block in range(0, num_blocks): - for ind, res in enumerate(fg): - slope_results[block, ind, :] = res.aperiodic_params - for band_label, band_range in bands.items(): - results[band_label][block, ind, :] = get_band_peak(res.peak_params, band_range, True) - -################################################################################################### -# Plotting topologies -# -------------- +# Plotting topographies +# --------------------- # -# The following will illustrate the process of exploring the previoulsy generated +# The following will illustrate the process of exploring the previously generated # oscillatory values ################################################################################################### @@ -161,20 +155,20 @@ def plot_topo_colorbar(vmin, vmax, label): norm = colors.Normalize(vmin=vmin, vmax=vmax) cb1 = colorbar.ColorbarBase(plt.gca(), cmap=cmap, - norm=norm, orientation='vertical')all_fg.plot() + norm=norm, orientation='vertical') ################################################################################################### -# Settings: In this example, we will be plotting the Alpha Central Frequency -# and oscillatory slope +# Settings: In this example, we will be plotting the Alpha Center Frequency +# and oscillatory exponent band = 'alpha' cur_data = results[band] -feats = ["CFS", "AMPS", "BWS"] + topo_dat = np.mean(cur_data,0) ################################################################################################### -# Looking at the alpha Central Frequeuncy +# Looking at the alpha Center Frequeuncy print('CURRENT FEATURE:', feats[0]) disp_dat = topo_dat[:,0] @@ -191,11 +185,11 @@ def plot_topo_colorbar(vmin, vmax, label): ################################################################################################### # Looking at the Oscillatory -cur_data = slope_results +cur_data = exponent_results topo_dat = np.mean(cur_data,0) -print('CURRENT FEATURE:', slope_feat[1]) +print('CURRENT FEATURE:', exponent_feat[1]) disp_dat = topo_dat[:,1] inds = np.where(np.isnan(disp_dat)) @@ -205,19 +199,20 @@ def plot_topo_colorbar(vmin, vmax, label): vmin, vmax, = disp_dat.min() - vbuffer, disp_dat.max() + vbuffer fig, ax = plt.subplots() -plot_topo_colorbar(vmin, vmax, slope_feat[1]) +plot_topo_colorbar(vmin, vmax, exponent_results[1]) mne.viz.plot_topomap(disp_dat, raw.info, vmin=vmin, vmax=vmax, cmap=cm.viridis, contours=0, axes=ax) ################################################################################################### # Alternative epoching methods -# -------------- +# ---------------------------- # # Antother way to epoch your data based on time by avaeraging across events as suppose to # selecting specifc time sections ################################################################################################### +tmin, tmax = -0.2, 0.5 epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, reject=reject) -epochs = mne.Epochs(raw_no_ref, **epochs_params).average() \ No newline at end of file +epochs = mne.Epochs(raw, **epochs_params).average() \ No newline at end of file From 58da8146aed01bafc186e168ce07e4f1a38ceaab Mon Sep 17 00:00:00 2001 From: Luyanda Mdanda Date: Fri, 19 Apr 2019 14:03:15 -0700 Subject: [PATCH 4/5] 2nd draft of MNE example, hot-fix --- examples/plot_mne_example.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/examples/plot_mne_example.py b/examples/plot_mne_example.py index e3888ba5b..02d46122f 100644 --- a/examples/plot_mne_example.py +++ b/examples/plot_mne_example.py @@ -132,6 +132,17 @@ # Creating dictionaries to store all the periodic properties across frequencies exponent_results = np.zeros(shape=[num_blocks, n_channels, len(aperiodic_feats)]) + +################################################################################################### + +# Populating periodic and aperiodic values +for block in range(0, num_blocks): + for ind, res in enumerate(fg): + exponent_results[block, ind, :] = res.aperiodic_params + for band_label, band_range in bands.items(): + results[band_label][block, ind, :] = get_band_peak(res.peak_params, band_range, True) + + ################################################################################################### # Plotting topographies # --------------------- @@ -189,7 +200,7 @@ def plot_topo_colorbar(vmin, vmax, label): topo_dat = np.mean(cur_data,0) -print('CURRENT FEATURE:', exponent_feat[1]) +print('CURRENT FEATURE:', aperiodic_feats[1]) disp_dat = topo_dat[:,1] inds = np.where(np.isnan(disp_dat)) From bf3618deba1ba06bf73ef4f7a9e006cd8f278c1d Mon Sep 17 00:00:00 2001 From: Tom Date: Fri, 19 Apr 2019 14:20:26 -0700 Subject: [PATCH 5/5] Misc style & doc updates. --- examples/plot_mne_example.py | 92 +++++++++++++----------------------- 1 file changed, 34 insertions(+), 58 deletions(-) diff --git a/examples/plot_mne_example.py b/examples/plot_mne_example.py index 02d46122f..03261f9ba 100644 --- a/examples/plot_mne_example.py +++ b/examples/plot_mne_example.py @@ -2,8 +2,7 @@ Using fooof with MNE ==================== -This examples illustrates how to use fooof with MNE -and create topographical plots +This examples illustrates how to use fooof with MNE and create topographical plots. This tutorial does require that you have `MNE `_ installed. If you don't already have @@ -13,8 +12,7 @@ ################################################################################################### - -#General imports +# General imports import numpy as np import pandas as pd from matplotlib import cm, colors, colorbar @@ -30,7 +28,6 @@ from fooof import FOOOF, FOOOFGroup from fooof.analysis import * - ################################################################################################### # Load & Check MNE Data # --------------------- @@ -74,65 +71,61 @@ reject = dict(eeg=180e-6) event_id = {'left/auditory': 1} events = mne.read_events(event_fname) -epochs = mne.Epochs(raw, events=events, event_id=event_id, - tmin=5, tmax=125, baseline=None, preload=True) +epochs = mne.Epochs(raw, events=events, event_id=event_id, tmin=5, tmax=125, + baseline=None, preload=True, verbose=False) ################################################################################################### # Creating Power Spectra Densities -psd, freqs = mne.time_frequency.psd_welch(epochs, fmin=1., fmax=50., n_fft=2000, - n_overlap=250, n_per_seg=500) +spectra, freqs = mne.time_frequency.psd_welch(epochs, fmin=1., fmax=50., n_fft=2000, + n_overlap=250, n_per_seg=500) ################################################################################################### -# fooofgroup +# FOOOFGroup # ---------- # -# The FOOOFGroup object is used to fit -# FOOOF models across the power spectra. A list of FOOOFGroup objects is returned. +# The FOOOFGroup object is used to fit FOOOF models across the power spectra. # ################################################################################################### -# Initialize a FOOOF group object, with desired settings +# Initialize a FOOOFGroup object, with desired settings fg = FOOOFGroup(peak_width_limits=[1, 6], min_peak_height=0.075, max_n_peaks=6, peak_threshold=1, verbose=False) ################################################################################################### -# Shaping the PSDs so that we can run fooofgroup -fooof_psd = np.squeeze(psd[0,:,:]) -n_channels, n_freq = fooof_psd.shape +# Selecting the first epoch of data to FOOOF +spectra = np.squeeze(spectra[0,:,:]) +n_channels, n_freq = spectra.shape num_blocks = len(mne.read_events(event_fname)) # Setting frequency range freq_range = [3, 35] -# This returns a list of FOOOFGRoup objects -fg.fit(freqs, fooof_psd, freq_range) +# Fit the FOOOF model across all channels +fg.fit(freqs, spectra, freq_range) fg.plot() ################################################################################################### -# Periodic features +# Define labels for periodic & aperiodic features feats = ["CFS", "PWS", "BWS"] - -# Aperiodic features aperiodic_feats = ["Offset","Exponent"] # Define bands of interest -bands = {'theta': [3,7], - 'alpha': [7,14], - 'beta': [15,30]} +bands = {'theta': [3, 7], + 'alpha': [7, 14], + 'beta': [15, 30]} -# Creating dictionaries to store all the aperiodic properties across frequencies +# Create dictionaries to store all the periodic properties across frequencies results = {} for band_name in bands.keys(): - results[band_name] = np.zeros(shape=[ num_blocks, n_channels, len(feats)]) + results[band_name] = np.zeros(shape=[num_blocks, n_channels, len(feats)]) -# Creating dictionaries to store all the periodic properties across frequencies +# Creating dictionaries to store all the aperiodic properties across frequencies exponent_results = np.zeros(shape=[num_blocks, n_channels, len(aperiodic_feats)]) - ################################################################################################### # Populating periodic and aperiodic values @@ -142,23 +135,18 @@ for band_label, band_range in bands.items(): results[band_label][block, ind, :] = get_band_peak(res.peak_params, band_range, True) - ################################################################################################### -# Plotting topographies +# Plotting Topographies # --------------------- # -# The following will illustrate the process of exploring the previously generated -# oscillatory values +# Now we can plot the extracted FOOOF features across all channels. +# ################################################################################################### -# This is the function required to create color bars for topographies + def plot_topo_colorbar(vmin, vmax, label): - """ - Create a colorbar for the topography plots - vmin: int - vmax: int - label: str - """ + """Helper function to create colorbars for the topography plots.""" + fig = plt.figure(figsize=(2, 3)) ax1 = fig.add_axes([0.9, 0.25, 0.15, 0.9]) @@ -168,10 +156,13 @@ def plot_topo_colorbar(vmin, vmax, label): cb1 = colorbar.ColorbarBase(plt.gca(), cmap=cmap, norm=norm, orientation='vertical') +################################################################################################### +# +# In this example, we will be plotting the alpha center frequency and oscillatory exponent. + ################################################################################################### -# Settings: In this example, we will be plotting the Alpha Center Frequency -# and oscillatory exponent +# Settings to grab the alpha center frequency band = 'alpha' cur_data = results[band] @@ -179,7 +170,7 @@ def plot_topo_colorbar(vmin, vmax, label): ################################################################################################### -# Looking at the alpha Center Frequeuncy +# Looking at the alpha center frequeuncy print('CURRENT FEATURE:', feats[0]) disp_dat = topo_dat[:,0] @@ -195,7 +186,7 @@ def plot_topo_colorbar(vmin, vmax, label): ################################################################################################### -# Looking at the Oscillatory +# Looking at the aperiodic exponent cur_data = exponent_results topo_dat = np.mean(cur_data,0) @@ -212,18 +203,3 @@ def plot_topo_colorbar(vmin, vmax, label): fig, ax = plt.subplots() plot_topo_colorbar(vmin, vmax, exponent_results[1]) mne.viz.plot_topomap(disp_dat, raw.info, vmin=vmin, vmax=vmax, cmap=cm.viridis, contours=0, axes=ax) - - -################################################################################################### -# Alternative epoching methods -# ---------------------------- -# -# Antother way to epoch your data based on time by avaeraging across events as suppose to -# selecting specifc time sections - -################################################################################################### - -tmin, tmax = -0.2, 0.5 -epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, - reject=reject) -epochs = mne.Epochs(raw, **epochs_params).average() \ No newline at end of file