EEG analysis tool MNE tutorial

Posted by TheKiller on Tue, 11 Jan 2022 11:27:56 +0100

Source: EEG Analysis series | MNE Python summary

catalogue

1. Installation and use

2. Data structure Raw

3.Epoch

4.Evoked

5. Cases

(1) Import tool library

(2) Load data

(3) Test data

(4) Extracting EEG features

(5) Forecast

6. Read set and locs file

7. Read edf file

8. Reference electrode and Application

9. Signal space projection SSP

10. Continuous notes

11.Epoch and its visualization

12. Signal filtering

1. Installation and use

Under windows system, install MNE Python based on Anaconda. See the data for specific steps, which are not shown here.

2. Data structure Raw

Raw object is used to store continuous data, and the core data is n_channels and times, raw Info also includes ch_names, chs, etc.

Generally, the data access method of raw is as follows:

data, times = raw[picks, time_slice]

Where picks is the index selected according to conditions, time_slice is a time slice (which can control the size of EEG data time period). To access all data in raw, the code is as follows:

data,times=raw[:]

data,times=raw[:,:]

Code examples are as follows:

# Select a good signal according to the type
# Set the acquisition signal time range
picks = mne.pick_types(raw.info, meg=True, exclude='bads')
t_idx = raw.time_as_index([10., 20.])
data, times = raw[picks, t_idx[0]:t_idx[1]]
plt.plot(times,data.T)
plt.title("Sample channels")

# Set sampling frequency: sfreq
# raw returns the data and time points of the selected channel and time period,
sfreq=raw.info['sfreq']
data,times=raw[:5,int(sfreq*1):int(sfreq*3)]
plt.plot(times,data.T)
plt.title("Sample channels")

# Plot the power spectral density of each channel
raw.plot_psd()
plt.show()

# Draw SSP vector diagram
raw.plot_projs_topomap()
plt.show()

# Plot channel spectrum
raw.plot_psd_topo()
plt.show()

# Draw electrode position
raw.plot_sensors()
plt.show()

3.Epoch

Epoch: extract some signals of specific time window from continuous EEG signals

For example, suppose we have a signal x with a length of 60s and a sampling frequency of 1 Hz.
The matrix of EEG signal is expressed as 1x60 matrix. If the signal is divided into some 2s signals, there will be 30 peoch (every 2s in the signal is an epoch)

4.Evoked

Evoked potential(EP) evoked potential or evoked response refers to the potential of a specific pattern recorded from the nervous system of human or other animals, especially the specific part of the brain, after the occurrence of stimulation such as flash or pure tone. Different forms and types of stimulation will produce different types of potentials.

5. Cases

(1) Import tool library

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.datasets.sleep_physionet.age import fetch_data    # Get EEG data
from mne.time_frequency import psd_welch                   
 
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer

(2) Load data

PSG.edf contains polysomnography. Raw data from EEG helmet, Hypnogram The EDF contains comments for expert records. Merge the two objects into MNE io. In the raw object, events can be extracted according to the description of the annotation to obtain time segments (epochs).

"""
Can pass
mne.datasets.sleep_physionet.age.fetch_data(subjects,recording,path)
To get PhysioNet Polysomnogram dataset file.
subjects: Indicates which subject objects you want to use. The range of available subject objects is 0-19. 
recording: Indicates the number of the night record(Indexes),Valid values are:[1],[2]or[1,2]. 

path:PhysioNet The storage address of data. If it is not given, the default storage address of data will be loaded;
If there is no data at the address where the default data set is stored, relevant data will be downloaded from the network.
"""
# Select two subjects, Alice and Bob (the name is not the real name in the experiment, but it is a temporary name for convenience)
ALICE, BOB = 0, 1
# Load the experiment data files of Alice and Bob
[alice_files, bob_files] = fetch_data(subjects=[ALICE, BOB], recording=[1])

# Channel name mapping
mapping = {'EOG horizontal': 'eog',
           'Resp oro-nasal': 'misc',
           'EMG submental': 'misc',
           'Temp rectal': 'misc',
           'Event marker': 'misc'}

#Read the edf file of ALICE and its corresponding annotation file
raw_train = mne.io.read_raw_edf(alice_files[0])
annot_train = mne.read_annotations(alice_files[1])

raw_train.set_annotations(annot_train, emit_warning=False)
raw_train.set_channel_types(mapping)

# Draw the continuous channel data waveform starting from 0s with a time window length of 40s
raw_train.plot(duration=40, scalings='auto')
plt.show()

Only five phases are used here: wake-up (W), phase 1, phase 2, phase 3 / 4, and REM sleep (R). To do this, use MNE events_ from_ Event in annotations()_ ID parameter to select the events we are interested in and associate the event identifier with each event.

annotation_desc_2_event_id = {'Sleep stage W': 1,
                              'Sleep stage 1': 2,
                              'Sleep stage 2': 3,
                              'Sleep stage 3': 4,
                              'Sleep stage 4': 4,
                              'Sleep stage R': 5}

events_train, _ = mne.events_from_annotations(
    raw_train, event_id=annotation_desc_2_event_id, chunk_duration=30.)

# Create a new event_id to unify phases 3 and 4
event_id = {'Sleep stage W': 1,
            'Sleep stage 1': 2,
            'Sleep stage 2': 3,
            'Sleep stage 3/4': 4,
            'Sleep stage R': 5}

# Draw event data
mne.viz.plot_events(events_train, event_id=event_id,
                    sfreq=raw_train.info['sfreq'])

# Keep the color code for further drawing
stage_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

Create epochs (time slices) from data based on events in comments

tmax = 30. - 1. / raw_train.info['sfreq']  # tmax in included
"""
Created from tmin=0 Start, to tmax Up to epochs
"""
epochs_train = mne.Epochs(raw=raw_train, events=events_train,
                          event_id=event_id, tmin=0., tmax=tmax, baseline=None)

print(epochs_train)

(3) Test data

raw_test = mne.io.read_raw_edf(bob_files[0])
annot_test = mne.read_annotations(bob_files[1])
raw_test.set_annotations(annot_test, emit_warning=False)
raw_test.set_channel_types(mapping)
events_test, _ = mne.events_from_annotations(
    raw_test, event_id=annotation_desc_2_event_id, chunk_duration=30.)
epochs_test = mne.Epochs(raw=raw_test, events=events_test, event_id=event_id,
                         tmin=0., tmax=tmax, baseline=None)

print(epochs_test)

Using feature engineering, PSD

fig, (ax1, ax2) = plt.subplots(ncols=2)

# iterate over the subjects
stages = sorted(event_id.keys())
for ax, title, epochs in zip([ax1, ax2],
                             ['Alice', 'Bob'],
                             [epochs_train, epochs_test]):

    for stage, color in zip(stages, stage_colors):
        epochs[stage].plot_psd(area_mode=None, color=color, ax=ax,
                               fmin=0.1, fmax=20., show=False,
                               average=True, spatial_colors=False)
    ax.set(title=title, xlabel='Frequency (Hz)')
ax2.set(ylabel='uV^2/hz (dB)')
ax2.legend(ax2.lines[2::3], stages)
plt.tight_layout()
plt.show()

(4) Extracting EEG features

def eeg_power_band(epochs):
    """Feature extraction of EEG relative power band
    This function accepts a""mne.Epochs"Object,
    And based on scikit-learn Compatible relative power creation in a specific frequency band EEG features.
    Parameters
    ----------
    epochs : Epochs
        The data.
    Returns
    -------
    X : numpy array of shape [n_samples, 5]
        Transformed data.
    """
    # Specific frequency band
    FREQ_BANDS = {"delta": [0.5, 4.5],
                  "theta": [4.5, 8.5],
                  "alpha": [8.5, 11.5],
                  "sigma": [11.5, 15.5],
                  "beta": [15.5, 30]}

    psds, freqs = psd_welch(epochs, picks='eeg', fmin=0.5, fmax=30.)
    # Normalized PSDs
    psds /= np.sum(psds, axis=-1, keepdims=True)

    X = []
    for fmin, fmax in FREQ_BANDS.values():
        psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
        X.append(psds_band.reshape(len(psds), -1))

    return np.concatenate(X, axis=1)

(5) Forecast

Multi classification using scikit learn

pipe = make_pipeline(FunctionTransformer(eeg_power_band, validate=False),
                     RandomForestClassifier(n_estimators=100, random_state=42))
# train
y_train = epochs_train.events[:, 2]
pipe.fit(epochs_train, y_train)

# forecast
y_pred = pipe.predict(epochs_test)

# Assessment accuracy
y_test = epochs_test.events[:, 2]
acc = accuracy_score(y_test, y_pred)

print("Accuracy score: {}".format(acc))
print(classification_report(y_test, y_pred, target_names=event_id.keys()))

6. Read set and locs file

See information

7. Read edf file

mne.io.read_raw_edf(input_fname,  # File storage address             
     eog=None,     #The channel name or should be specified as the index list of EOG channels. The value corresponds to the electrode in the document. 
     misc=None,    #The channel name or should be specified as an indexed list of MISC channels. The value corresponds to the electrode in the document.      
     stim_channel='auto',                    
     exclude=(),        #The name of the channel to exclude.
     preload=False)     #Preload

Code routine:

# Import tool library
from mne.io import concatenate_raws, read_raw_edf
import matplotlib.pyplot as plt
import mne

# Load EDF file
raw=read_raw_edf("Affaf Ikram 20121020 1839.L1.edf",preload=False)

# Get events in raw data
events_from_annot, event_dict = mne.events_from_annotations(raw)
print(event_dict)
print(events_from_annot)

# Get the corresponding event
custom_mapping = {'stm+':5, 'stm-': 6}
(events_from_annot,
 event_dict) = mne.events_from_annotations(raw, event_id=custom_mapping)
print(event_dict)
print(events_from_annot)

# Plot events
fig = mne.viz.plot_events(events_from_annot, sfreq=raw.info['sfreq'],
                          first_samp=raw.first_samp, event_id=event_dict)
fig.subplots_adjust(right=0.7)
epochs = mne.Epochs(raw, events=events_from_annot, 
                    event_id=event_dict)

epochs.plot_image()

# Get sampling frequency
data,times=raw[:3,int(sfreq*1):int(sfreq*3)]
plt.plot(times,data.T)
plt.title("Sample channels")

8. Reference electrode and Application

See information

9. Signal space projection SSP

See information

10. Continuous notes

See information

11.Epoch and its visualization

See information

12. Signal filtering

# Trap wave
import numpy as np
import mne
from mne.datasets import sample

data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
proj_fname = data_path + '/MEG/sample/sample_audvis_eog_proj.fif'

"""
Extract data between 0 and 20 seconds
"""
tmin, tmax = 0, 20

"""
Read raw data
 Save memory by cutting raw data before loading
"""
raw = mne.io.read_raw_fif(raw_fname)
raw.crop(tmin, tmax).load_data()
raw.info['bads'] = ['MEG 2443', 'EEG 053']  # bads + 2 more

"""
Set the frequency at 2 Hz To 300 Hz between
"""
fmin, fmax = 2, 300
"""
FFT Size n_fft,Ideally, it is a power of 2
"""
n_fft = 2048

# Select a subset of channels
selection = mne.read_selection('Left-temporal')
picks = mne.pick_types(raw.info, meg='mag', eeg=False, eog=False,
                       stim=False, exclude='bads', selection=selection)
raw.plot_psd(area_mode='range', tmax=10.0, picks=picks, average=False)

# Trap wave
raw.notch_filter(np.arange(60, 241, 60), picks=picks, fir_design='firwin')
raw.plot_psd(area_mode='range', tmax=10.0, picks=picks, average=False)

# Low pass filtering below 50 Hz
raw.filter(None, 50., fir_design='firwin')
raw.plot_psd(area_mode='range', tmax=10.0, picks=picks, average=False)

# High pass filtering
raw.filter(1., None, fir_design='firwin')
raw.plot_psd(area_mode='range', tmax=10.0, picks=picks, average=False)

# Band pass filtering in the range of 1 Hz-50 Hz
raw.filter(1, 50., fir_design='firwin')
raw.plot_psd(area_mode='range', tmax=10.0, picks=picks, average=False)

Topics: Python