Extracting power bands (theta, delta, alpha, beta bands) from the EEG signal

💬 Support & Discussions
1

How can I extract and store signals of the different frequency band from the EEG signal?

  • MNE-Python version: 0.19.2
  • Google Colab:

**

for each_subject in list_ids:
  
  path_0=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[0]
  path_1=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[1]
  path_2=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[2]

  print('ID_list', len(list_ids))
  
  
  df_subject_channels=pd.read_csv(path_0,
                              sep=',',
                              decimal=".")
  print('size of channel data', df_subject_channels.shape)
  #channels labels 
  ch_labels=list(df_subject_channels.labels)
  #apply montage
  internal_montage = mne.channels.make_standard_montage(MONTAGE)
  print('montage', internal_montage)
  
  
  #create info for object
  info = mne.create_info(ch_names=ch_labels, 
                            sfreq=s_freq, 
                            ch_types='eeg', 
                            montage=internal_montage)
  
  #signal
  #load and scale signal
  dat_test=np.loadtxt(path_1, delimiter=',')*SCALE
  print('size of signal data', dat_test.shape)


  #Create the MNE Raw data object
  raw = mne.io.RawArray(dat_test, info)

  #create in stimuation channel
  stim_info = mne.create_info(['stim'], s_freq, 'stim')
  #create zero signal to store stimulus
  stim_raw = mne.io.RawArray(np.zeros(shape=[1, len(raw._times)]), stim_info)

  #add stim channle to raw signal
  raw.add_channels([stim_raw], force_update_info=True)
  #events
  #read csv of events
  df_subject_event=pd.read_csv(path_2,
                                delimiter=",",
                                decimal=".")

  #fake structure of events
  evs = np.empty(shape=[0, 3])

  #from HBT, the signals were already marked each 20 seconds.
  for each_element in df_subject_event.values[1:len(df_subject_event)-1]:
      if('break cnt'!=each_element[0]):
          if(int(each_element[0])==EVENT_ID):
              evs = np.vstack((evs, np.array([each_element[1], 0,
                                                    int(each_element[0])])))

          
  #print(evs)
  # Add events to data object
  raw.add_events(evs, stim_channel='stim')

  #Check events
  #print(mne.find_events(raw))

  #detect flat channels
  flat_chans = np.mean(raw._data[:len(ch_labels), :], axis=1) == 0

  # Interpolate bad channels
  raw.info['bads'] = \
                  list(np.array(raw.ch_names[:len(ch_labels)])[flat_chans])
  print('Bad channels: ', raw.info['bads'])
  raw.interpolate_bads()

  # Get good eeg channel indices
  #eeg_chans = mne.pick_types(raw.info, meg=False, eeg=True)

  #resample to have to 250 hz, 
  #this will allow us to compare with
  #the ADHD dataset.
  raw.resample(R_S_FREQ, npad='auto')

  #set reference to Cz
  raw.set_eeg_reference(ref_channels=['Cz'])

  raw.drop_channels(['Cz'])




  events = mne.find_events(raw)

  picks = mne.pick_types(raw.info, eeg=True, eog=False, stim=True)

  #print('pick.info', picks)
  

  # Construct Epochs
  event_id, tmin, tmax = 30, -1., 3.
  baseline = (None, 0)
  epochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                    baseline=baseline,
                    preload=True)
  print('epoch shape',len(epochs))
  
  epochs.resample(200., npad='auto')  # resample to reduce computation time

 
  mne.events_from_annotations(raw)
  sf=250
  time=np.arange(len(epochs))/sf
  print('Time', time)
''' calculate the PSD and extract power bands from it'''

**.
Please help me extract features from the signals!
When I try to calculate the PSD of a signal

 from scipy import signal
  win = 4 * sf

  freqs, psd = signal.welch(raw, sf, nperseg=20)

  pritnt('Freqs shape',freqs.shape)

  pritnt('Psd shape',Psd.shape)

I get an error → ValueError: All picks must be < n_channels (111), got 111

2

an example of extracting frequency bands is here:

https://mne.tools/stable/auto_examples/time_frequency/time_frequency_global_field_power.html

3

Thanks a lot for your support @drammock !

I managed to plot the graphs for different frequency bands with the help of the link you provided. The issue with the graph is that it only shows changes in the theta band whereas gamma, beta and alpha bands do not show any changes, I tried to plot graphs for 10 different subjects, for every subject I only see changes in the theta band. Could you help me clarify the issue? Please check out the updated code for reference.

import numpy as np
import mne
import pandas as pd
import matplotlib.pyplot as plt
from mne.stats import bootstrap_confidence_interval
from mne.baseline import rescale
import configparser

PATH_SIGNALS_CSV='/content/drive/MyDrive/adhd-detector/data/'
CSV_SUBJECTS_IDS= '/content/drive/MyDrive/adhd-detector/data/subjects_id/subjects_id.csv'
FILES_USE=['RestingState_chanlocs.csv','RestingState_data.csv','RestingState_event.csv']


subjects_id=pd.read_csv(CSV_SUBJECTS_IDS)
list_ids=subjects_id.id.values.tolist()
i=0
for each_subject in list_ids:
      path_0=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[0]
      path_1=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[1]
      path_2=PATH_SIGNALS_CSV+each_subject+'_'+FILES_USE[2]

      dat_test=np.loadtxt(path_1, delimiter=',')
      df_subject_channels=pd.read_csv(path_0,
                              sep=',',
                              decimal=".")
      print('size of channel data', df_subject_channels.shape)
      #channels labels 
      ch_labels=list(df_subject_channels.labels)
      internal_montage = mne.channels.make_standard_montage('GSN-HydroCel-129')
      s_freq=250
      #create info for object
      info = mne.create_info(ch_names=ch_labels, 
                                sfreq=s_freq, 
                                ch_types='eeg', 
                                montage=internal_montage)
      raw= mne.io.RawArray(dat_test, info)
      # let's explore some frequency bands
      iter_freqs = [
          ('Theta', 4, 8),
          ('Alpha', 8, 16),
          ('Beta', 16, 32),
          ('Gamma', 32, 45)
      ]

      #create in stimuation channel
      s_freq=500
      r_s_freq=250
      stim_info = mne.create_info(['stim'], s_freq, 'stim')
      #create zero signal to store stimulus
      stim_raw = mne.io.RawArray(np.zeros(shape=[1, len(raw._times)]), stim_info)

      #add stim channle to raw signal
      raw.add_channels([stim_raw], force_update_info=True)

      #events
      #read csv of events
      df_subject_event=pd.read_csv(path_2,
                                    delimiter=",",
                                    decimal=".")

      #fake structure of events
      evs = np.empty(shape=[0, 3])
      # set epoching parameters
      event_id, tmin, tmax = 30, -1., 20.
      baseline = None
      #from HBT, the signals were already marked each 20 seconds.
      for each_element in df_subject_event.values[1:len(df_subject_event)-1]:
          if('break cnt'!=each_element[0]):
              if(int(each_element[0])==event_id):
                  evs = np.vstack((evs, np.array([each_element[1], 0,
                                                        int(each_element[0])])))
      print(type(evs))
      # Add events to data object
      raw.add_events(evs, stim_channel='stim')
      #Check events
      print(mne.find_events(raw))

      #detect flat channels
      flat_chans = np.mean(raw._data[:len(ch_labels), :], axis=1) == 0
      # Interpolate bad channels
      raw.info['bads'] = \
                      list(np.array(raw.ch_names[:len(ch_labels)])[flat_chans])
      print('Bad channels: ', raw.info['bads'])
      raw.interpolate_bads()


      #resample to have to 250 hz, 
      #this will allow us to compare with
      #the HDHD dataset.
      raw.resample(r_s_freq, npad='auto')

      #set reference to Cz
      raw.set_eeg_reference(ref_channels=['Cz'])

      raw.drop_channels(['Cz'])
      #return Raw object from mne class
      evs=evs.astype(int)
      frequency_map=list()
      for band, fmin, fmax in iter_freqs:
          # (re)load the data to save memory
          #raw = mne.io.read_raw_fif(raw_fname)
          raw.pick_types(eeg=True)  # we just look at gradiometers
          raw.load_data()

          # bandpass filter
          raw.filter(fmin, fmax, n_jobs=1,  # use more jobs to speed up.
                    l_trans_bandwidth=1,  # make sure filter params are the same
                    h_trans_bandwidth=1)  # in each band and skip "auto" option.
    
          # epoch
          epochs = mne.Epochs(raw, evs, event_id, tmin, tmax, baseline=baseline,
                              preload=True)
          # remove evoked response
          epochs.subtract_evoked()

          # get analytic signal (envelope)
          epochs.apply_hilbert(envelope=True)
          frequency_map.append(((band, fmin, fmax), epochs.average()))
          del epochs
      del raw
      # Helper function for plotting spread
      def stat_fun(x):
          """Return sum of squares."""
          return np.sum(x ** 2, axis=0)

      # Plot
      print(i,each_subject)
      fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True)
      colors = plt.get_cmap('winter_r')(np.linspace(0, 1, 4))
      plot_title=each_subject
      
      for ((freq_name, fmin, fmax), average), color, ax in zip(
              frequency_map, colors, axes.ravel()[::-1]):
          times = average.times * 1e3
          gfp = np.sum(average.data ** 2, axis=0)
          gfp = mne.baseline.rescale(gfp, times, baseline=(0, 0.1))
          ax.plot(times, gfp, label=freq_name, color=color, linewidth=2.5)
          
          ax.axhline(0, linestyle='--', color='grey', linewidth=2)
          ci_low, ci_up = bootstrap_confidence_interval(average.data, random_state=0,
                                                        stat_fun=stat_fun)
          ci_low = rescale(ci_low, average.times, baseline=(None, 0))
          ci_up = rescale(ci_up, average.times, baseline=(None, 0))
          ax.fill_between(times, gfp + ci_up, gfp - ci_low, color=color, alpha=0.3)
          ax.grid(True)
          
          ax.set_ylabel('GFP')
          ax.annotate('%s (%d-%dHz)' % (freq_name, fmin, fmax),
                      xy=(0.95, 0.8),
                      horizontalalignment='right',
                      xycoords='axes fraction')
          ax.set_xlim(-1000, 3000)

      i=i+1
      axes.ravel()[-1].set_xlabel('Time [ms]')
      ax.set_title(plot_title)

Here is the graph generated by the script:

plot
plot717×457 22.7 KB

Any help will be highly appreciated.

4

Inside your for-loop (for band, fmin, fmax in iter_freqs) you are filtering the data. This changes the data in-place, so after filtering once for theta, there is no energy left in the other bands on subsequent loop iterations. In the original example this was not a problem because the raw file was reloaded on each loop iteration, but you’ve commented out that line. Easiest fix is probably to change your raw.filter(...) line to say filtered_raw = raw.copy().filter(...) and then use filtered_raw in the epoching line.

5

Thank you @drammock, the issue is solved! Best wishes. :slight_smile: