Source: EEG Analysis series | MNE Python summary
catalogue
8. Reference electrode and Application
9. Signal space projection SSP
11.Epoch and its visualization
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) #PreloadCode 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)