2nd draft of MNE example

Author: LuyandaMdanda

examples/plot_mne_example.py

@@ -1,11 +1,12 @@
 """
 Using fooof with MNE
-=======================
+====================
 
-This examples illustrates how to use fooof with `MNE
-<https://mne-tools.github.io/>`_ 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
+<https://mne-tools.github.io/>`_ installed. If you don't already have
 MNE, you can follow instructions to get it `here
 <https://mne-tools.github.io/stable/getting_started.html>`_.
 """
@@ -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,31 +65,27 @@
 
 ###################################################################################################
 
-# 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
+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,52 +94,49 @@
 ###################################################################################################
 
 # 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)
 fg.plot()
 
 ###################################################################################################
 
+# 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()
+epochs = mne.Epochs(raw, **epochs_params).average()