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方法)。