MNE-python《handling bad channels》

# handling bad channels
import os
from copy import deepcopy
import numpy as np
import mne
import matplotlib.pyplot as plt

# get the sample data
sample_data_path = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_path, 'MEG', 'sample',
                                    'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False)

# marking bad channels

# 1.output the stored bad channels
print(raw.info['bads'])
# pick the 'EEG 05'channel for showing in plot
picks = mne.pick_channels_regexp(raw.ch_names, regexp='EEG 05.')
# picks:这里mne.pick_channels_regexp返回的是一个 int型array，
# 存储raw对象ch_names中同regexp value值相似或相同的channel name的Indices
# print('len of picks: ',  len(picks))  # from EEG 050-EEG059,Total number: 10

# raw.plot(order=picks, n_channels=len(picks))
# plt.show()

picks = mne.pick_channels_regexp(raw.ch_names, regexp='MEG 2..3')
# 同上，这里获取MEG 2..3，即MEG信号中四位且首位和末位分别是2和3的所有MEG信号
# raw.plot(order=picks, n_channels=len(picks))
# plt.show()

# add bad channels to raw.info['bad']
original_bads = deepcopy(raw.info['bads'])
# print(original_bads)
raw.info['bads'].append('EEG 050')  # add a single channel
raw.info['bads'].extend(['EEG 051', 'EEG 052'])  # add a list of channels
# print(raw.info['bads'])
bad_chan = raw.info['bads'].pop(-1)  # remove the last entry in the list
# print(bad_chan)
# recovery the initial bad channels
raw.info['bads'] = original_bads

# default is exclude='bads':
good_eeg = mne.pick_types(raw.info, meg=False, eeg=True)
all_eeg = mne.pick_types(raw.info, meg=False, eeg=True, exclude=[])
# print(np.setdiff1d(all_eeg, good_eeg))  # the index of 'EEG 053'channel
"""
np.setdiff1d: Find the set difference of two arrays.
def setdiff1d(ar1, ar2, assume_unique=False):
    Find the set difference of two arrays.

    Return the unique values in `ar1` that are not in `ar2`.
"""
# print(np.array(raw.ch_names)[np.setdiff1d(all_eeg, good_eeg)])

# print('raw chan: ', raw.info['nchan'])


# When to look for bad channels
raw2 = raw.copy()
# print('raw2 chan: ', raw2.info['nchan'])
raw2.info['bads'] = []
# print(raw2.info)
# print('raw2 chan: ', raw2.info['nchan'])
events = mne.find_events(raw2, stim_channel='STI 014')
# mne.Epochs(raw2, events=events)['2'] 2为event id，
# here it can be 1 or 2 or 3 or 4 or 5 or 32

# picks = mne.pick_channels_regexp(raw.ch_names, regexp='EEG 0.')
# print("len: ", len(picks))

epochs = mne.Epochs(raw2, events=events)['2'].average().plot()
# default: picks=None,but in the result the plot window contains three types channels——
# 1.EEG channels 2.Gradiometers channels 3.Magnetometers channels and not contains Stimulus channels, EOG channels and bad channels
# so in this case, the total number of channels contained in the end plot window is 366, contains 60 for EEG channels,
# 102 for Magnetometers channels, 204 for Gradiometers channels.
# but i carefully read the method'mne.Epochs',but i don't understand why the plot window doesn't contain Stimulus channels and EOG channels.
# maybe it automaticlly filters these channels for not affecting others channels.


# How interpolation works
# mne-python: interpolation——'插值'
# 引用mne官网的一段话肯翻译——
"""
In some cases simply excluding bad channels is sufficient 
(for example, if you plan only to analyze a specific sensor 
ROI, and the bad channel is outside that ROI). However, 
in cross-subject analyses it is often helpful to maintain 
the same data dimensionality for all subjects, and there is 
no guarantee that the same channels will be bad for all 
subjects. It is possible in such cases to remove each channel
that is bad for even a single subject, but that can lead to
a dramatic drop in data rank (and ends up discarding a fair
amount of clean data in the process). In such cases it is 
desirable to reconstruct bad channels by interpolating its 
signal based on the signals of the good sensors around them.

在某些情况下，仅排除坏通道就足够了（例如，如果您只计划分析特定的传感器ROI，
而坏通道在该ROI之外）。然而，在交叉学科的信号数据分析中，为所有受试者保持
相同的数据维度通常是有益的，并且不能保证相同的渠道对所有受试者都是有害的。
在这种情况下，可以删除对每个受试者都不好的单个通道，但这可能会导致数据排名
大幅下降（并最终在过程中丢弃相当数量级的干净数据）。在这种情况下，
我们希望通过基于周围良好传感器的信号插值其信号来重建不良信道。
"""

raw.crop(tmin=0, tmax=3).load_data()
eeg_data = raw.copy().pick_types(meg=False, eeg=True, exclude=[])
eeg_data_interp = eeg_data.copy().interpolate_bads(reset_bads=False)

for title, data in zip(['orig.', 'interp.'], [eeg_data, eeg_data_interp]):
    fig = data.plot(butterfly=True, color='#00000022', bad_color='r')
    fig.subplots_adjust(top=0.9)
    fig.suptitle(title, size='xx-large', weight='bold')
    plt.show()

grad_data = raw.copy().pick_types(meg='grad', exclude=[])
grad_data_interp = grad_data.copy().interpolate_bads(reset_bads=False)

for data in (grad_data, grad_data_interp):
    data.plot(butterfly=True, color='#00000009', bad_color='r')
    plt.show()