MNE学习笔记（一）：数据结构与容器（上）

声明：以下内容皆摘自MNE官方文档：https://mne-tools.github.io/stable/documentation.html

并融入了一些个人理解

 

（一）数据结构与容器（上）

一、Raw数据结构

1、定义：Raw数据类型的对象用来存储连续型数据，核心数据以2维numpy array（分别是channels和samples）的形式存储，除此之外还包含了Info对象。

2、索引数据

要想访问Raw对象内的数据，必须得对其进行索引。代码如下

# Extract data from the first 5 channels, from 1 s to 3 s.
sfreq = raw.info['sfreq']
data, times = raw[:5, int(sfreq * 1):int(sfreq * 3)]
_ = plt.plot(times, data.T)
_ = plt.title('Sample channels')

其中sfreq表示采样频率，raw返回所选信道以及时间段内的数据和时间点，分别赋值给data以及times（即raw对象返回的是两个array），效果如下图：

3、选择信道和样本的子集（进阶的索引）

可以直接通过信道的名称、类型或者时间范围来提取数据

# Pull all MEG gradiometer channels:
# Make sure to use .copy() or it will overwrite the data
meg_only = raw.copy().pick_types(meg=True)
eeg_only = raw.copy().pick_types(meg=False, eeg=True)

# The MEG flag in particular lets you specify a string for more specificity
grad_only = raw.copy().pick_types(meg='grad')

# Or you can use custom channel names
pick_chans = ['MEG 0112', 'MEG 0111', 'MEG 0122', 'MEG 0123']
specific_chans = raw.copy().pick_channels(pick_chans)
print(meg_only)
print(eeg_only)
print(grad_only)
print(specific_chans)

注意：

（1）提取之前调用copy方法，以防对原始数据的永久修改；

（2）按类型提取和按信道提取所采用的方法不同，分别是pick_types以及pick_channels；

（3）提取出来的结果是和raw一样的二维数据，访问方法同2。

我们也可以使用crop方法来直接指定时间范围

raw = raw.crop(0, 50)  # in seconds
print('New time range from', raw.times.min(), 's to', raw.times.max(), 's')

以及直接通过信道名来删除信道

nchan = raw.info['nchan']
raw = raw.drop_channels(['MEG 0241', 'EEG 001'])
print('Number of channels reduced from', nchan, 'to', raw.info['nchan'])

4、连接Raw对象

Raw对象可以通过调用append方法来连接。连接的对象必须满足信道数相同，以及Info结构可兼容

# Create multiple :class:`Raw <mne.io.RawFIF>` objects
raw1 = raw.copy().crop(0, 10)
raw2 = raw.copy().crop(10, 20)
raw3 = raw.copy().crop(20, 40)

# Concatenate in time (also works without preloading)
raw1.append([raw2, raw3])
print('Time extends from', raw1.times.min(), 's to', raw1.times.max(), 's')

二、Epochs数据结构

Epochs对象是一种把连续型数据作为时间段集合的表示方法，以形如(n_events, n_channels, n_times)的array形式存储。

1、创建Epoch对象

Epochs对象可以由以下三种方式创建：

（1）通过Raw对象和事件时间点（event times）

（2）通过存储为.fif文件的现成Epochs对象

（3）通过EpochsArray方法从头创建

以下代码示例为从含有事件（events）的原始数据集中创建Epochs对象

data_path = mne.datasets.sample.data_path()
# Load a dataset that contains events
raw = mne.io.read_raw_fif(
    op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif'))

# If your raw object has a stim channel, you can construct an event array
# easily
events = mne.find_events(raw, stim_channel='STI 014')

# Show the number of events (number of rows)
print('Number of events:', len(events))

# Show all unique event codes (3rd column)
print('Unique event codes:', np.unique(events[:, 2]))

# Specify event codes of interest with descriptive labels.
# This dataset also has visual left (3) and right (4) events, but
# to save time and memory we'll just look at the auditory conditions
# for now.
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2}

上面最后一行代码是给不同的events一个描述性名称，方便今后的使用。

注：这种方式构造的epochs的数据是无法使用的，除非通过get_data方法将数据显式读入内存；或者在读入raw数据时将preload参数设为True、

以下代码的作用：Expose the raw data as epochs, cut from -0.1 s to 1.0 s relative to the event onsets.

这句话我的理解是，把原始数据看成若干个epoch（时间片）的集合，以每个提前标注好的event（某种波形、病人的动作、症状发作点）为参照，取tmin到tmax这一段时间作为一个epoch来切分原数据。下面这段代码就是取每个event的前0.1秒和后1秒共1.1秒的长度作为一个epoch，为后续同类型的event的epochs叠加（average）分析做准备。

epochs = mne.Epochs(raw, events, event_id, tmin=-0.1, tmax=1,
                    baseline=(None, 0), preload=True)
print(epochs)

注：事件（event）描述的是某一种波形（症状）的起始点，它是一个三元组，第一个元素是以整数来描述的事件起始采样点（根据采样频率可以换算出时间）；第二个元素对应的是当前事件来源的刺激通道（stimulus channel）的先前值（previous value），大多数情况是0；第三个元素则是该事件的id。

事件的输出格式如下：

[[28771     0     1]
 [30450     0     1]
 [32101     0     1]
 [33712     0     1]
 [35428     0     1]]

---

 [[27977     0     2]
 [28771     0     1]
 [29652     0     2]
 [30450     0     1]
 [31240     0     2]]


At sample 28771 stim channel went from 0 to 1
At sample 30450 stim channel went from 0 to 1
At sample 32101 stim channel went from 0 to 1
At sample 33712 stim channel went from 0 to 1
At sample 35428 stim channel went from 0 to 1

2、查看Epoch对象

对于一个Eopchs对象，可以通过以下两种形式来查看当前epoch内的event相关信息

print(epochs.events[:3])
print(epochs.event_id)

输出如下

[[27977     0     2]
 [28771     0     1]
 [29652     0     2]]
{'Auditory/Left': 1, 'Auditory/Right': 2}

类似的，可以通过指定python的列表切片方式访问epochs对象内的events，或者直接通过events的描述性名称直接访问：

print(epochs[1:5])
print(epochs['Auditory/Right'])

注：'/'符号用来划分标签（tag），每个被'/'划分开的单词都可以作为检索的字符串，如直接搜索epochs['Right']也能得到结果

3、删除某个epoch

可以使用epochs.drop()方法或者epochs.drop_bad()来删除某个epoch

epochs.drop([0], reason='User reason')
epochs.drop_bad(reject=dict(grad=2500e-13, mag=4e-12, eog=200e-6), flat=None)
print(epochs.drop_log)
epochs.plot_drop_log()
print('Selection from original events:\n%s' % epochs.selection)
print('Removed events (from numpy setdiff1d):\n%s'
      % (np.setdiff1d(np.arange(len(events)), epochs.selection).tolist(),))
print('Removed events (from list comprehension -- should match!):\n%s'
      % ([li for li, log in enumerate(epochs.drop_log) if len(log) > 0]))

4、epoch平均叠加

使用Evoked类型对象可以实现epoch的叠加。通过调用mne.Epochs.average()方法可返回Evoked对象，average()方法可以通过参数指定需要的信道。

代码如下：

ev_left = epochs['Auditory/Left'].average()
ev_right = epochs['Auditory/Right'].average()

f, axs = plt.subplots(3, 2, figsize=(10, 5))
_ = f.suptitle('Left / Right auditory', fontsize=20)
_ = ev_left.plot(axes=axs[:, 0], show=False, time_unit='s')
_ = ev_right.plot(axes=axs[:, 1], show=False, time_unit='s')
plt.tight_layout()

以左上角图为例，该图提取了原数据中带有'Auditory/Left'事件标记的epochs，同时抽取了所有EEG信道进行平均值叠加。实际操作中可以指定更为详细的信道进行叠加（通过pick方法提取信道并传入average方法）。