Python-MNE全套教程(官网翻译)-连续数据的处理01:Raw数据结构

本教程详细介绍了Raw数据结构,包括如何从Raw对象加载、查询、子选择、导出和绘制数据。

老规矩,先导入:

import os

import matplotlib.pyplot as plt
import numpy as np

import mne

加载连续数据

MNE-Python数据结构基于来自Neuromag的.fif文件格式。本教程使用.fif格式的示例数据集,因此这里我们将使用函数mne.io.read_raw_fif()来加载原始数据。

​还有其他几个示例数据集,只需几行代码就可以下载,用的函数是mne.datasets。在这里我们使用mne.datasets.sample.data_path()下载Sampl”数据集,其中包含执行视听实验的受试者的EEG, MEG和结构MRI数据。下载完成后,data_path()将返回存放文件的文件夹位置。一旦有了文件路径,就可以使用read_raw_fif()加载数据。这将返回一个Raw对象,我们将其存储在一个名为Raw的变量中。

在MNE-Python中有几个示例数据集的data_path函数,如mne.datasets.kiloword.data_path()mne.datasets.spm_face.data_path()等。它们都将首先检查默认下载位置,查看数据集是否已经被下载。默认的下载位置也是可配置的,更多信息请参阅官方data_path函数文档。

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(
    sample_data_folder, "MEG", "sample", "sample_audvis_raw.fif"
)
raw = mne.io.read_raw_fif(sample_data_raw_file)

Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif…
Read a total of 3 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Range : 25800 … 192599 = 42.956 … 320.670 secs
Ready.

read_raw_fif()自动显示它正在加载的文件的一些信息。例如,这里它告诉我们,文件中有3个“projector”以及记录的数据,这些是SSP投影仪,用于计算用于从MEG信号中去除的环境噪声。除了加载过程中显示的信息,你还可以通过print raw对象来了解其基本细节:

print(raw)

<Raw | sample_audvis_raw.fif, 376 x 166800 (277.7 s), ~3.3 MB, data not loaded>

默认情况下,mne.io.read_raw_ *函数族不会将数据加载到内存中,磁盘上的数据是内存映射的。有些操作,例如筛选,需要将数据复制到物理内存中。为此,可以将preload=True参数传递给read_raw_fif(),但也可以随时使用load_data()方法将数据复制到RAM中。然而,由于是教程,我们将首先crop()原始对象到60秒,以便使用更少的内存。

raw.crop(tmax=60)

查询Raw对象

可以通过Raw类的各种属性和方法可以获得更多信息。Raw对象的一些有用属性包括通道名称列表(ch names)、以秒为单位的采样时间数组(times)和总采样次数(n次)等。

raw.info属性

raw.info中还存储了相当多的信息。它存储了一个类似于Python字典的对象(因为它有通过命名键访问的字段)。和Python字典一样,raw.info有一个.keys()方法,可以显示所有可用的字段名。与Python字典不同,打印raw.info将打印出每个字段数据的精确格式。

n_time_samps = raw.n_times
time_secs = raw.times
ch_names = raw.ch_names
n_chan = len(ch_names)  # note: there is no raw.n_channels attribute
print(
    f"the (cropped) sample data object has {n_time_samps} time samples and "
    f"{n_chan} channels."
)
print(f"The last time sample is at {time_secs[-1]} seconds.")
print("The first few channel names are {}.".format(", ".join(ch_names[:3])))
print()  # insert a blank line in the output

# some examples of raw.info:
print("bad channels:", raw.info["bads"])  # chs marked "bad" during acquisition
print(raw.info["sfreq"], "Hz")  # sampling frequency
print(raw.info["description"], "\n")  # miscellaneous acquisition info

print(raw.info)

the (cropped) sample data object has 36038 time samples and 376 channels.
The last time sample is at 60.000167471573526 seconds.
The first few channel names are MEG 0113, MEG 0112, MEG 0111.

bad channels: [‘MEG 2443’, ‘EEG 053’]
600.614990234375 Hz
acquisition (megacq) VectorView system at NMR-MGH

<Info | 21 non-empty values
acq_pars: ACQch001 110113 ACQch002 110112 ACQch003 110111 ACQch004 110122 …
bads: 2 items (MEG 2443, EEG 053)
ch_names: MEG 0113, MEG 0112, MEG 0111, MEG 0122, MEG 0123, MEG 0121, MEG …
chs: 204 Gradiometers, 102 Magnetometers, 9 Stimulus, 60 EEG, 1 EOG
custom_ref_applied: False
description: acquisition (megacq) VectorView system at NMR-MGH
dev_head_t: MEG device -> head transform
dig: 146 items (3 Cardinal, 4 HPI, 61 EEG, 78 Extra)
events: 1 item (list)
experimenter: MEG
file_id: 4 items (dict)
highpass: 0.1 Hz
hpi_meas: 1 item (list)
hpi_results: 1 item (list)
lowpass: 172.2 Hz
meas_date: 2002-12-03 19:01:10 UTC
meas_id: 4 items (dict)
nchan: 376
proj_id: 1 item (ndarray)
proj_name: test
projs: PCA-v1: off, PCA-v2: off, PCA-v3: off
sfreq: 600.6 Hz

raw.info的大多数字段反映了在获取时记录的元数据,不能更改。有一些例外,如raw.info['bads']raw.info[' projects '],但在大多数情况下,都有专用的MNE-Python函数或方法来安全地更新Info对象,如add_proj()来更新raw.info[' projects ']

时间、采样量和采样序数

Raw对象的一个常用方法是time_as_index(),它将时间以秒为单位转换为最接近该时间的采样数的整数索引。该方法也可以接受一个时间列表或数组,并返回一个索引数组。
因为是整数索引,所以time = 1和time = 2之间的样本量可能不同于time = 2和time = 3之间的样本量:

print(raw.time_as_index(20))
print(raw.time_as_index([20, 30, 40]), "\n")

print(np.diff(raw.time_as_index([1, 2, 3])))

[12012]
[12012 18018 24024]
[601 600]

修改Raw对象

Raw对象有许多方法可以就地修改。这对于method chaining很有用,例如,Raw .crop(…).pick(…).filter(…).plot(),但它也会在交互式分析期间产生问题。如果我们修改Raw对象以进行探索性绘图或分析,例如删除一些通道然后绘图,那么我们就需要重新加载数据并重复任何先前的处理步骤以撤销通道删除并尝试其他操作。因此我们要经常使用copy()方法。

虽然Raw对象本质上是将数据样本存储在NumPy数组(n_channels, n_timepoints)中,但就len()函数而言,Raw对象与NumPy数组的行为不同。len(raw)将返回时间点的数量而不是通道的数量。因此在本节中,你将看到len(raw.ch_names)获取通道的数量。

选择,删除和重新排序通道

改变一个raw对象的通道可以通过几种方式完成。作为第一个例子,我们将使用pick()方法将Raw对象限制为EEG和EOG通道:

eeg_and_eog = raw.copy().pick(picks=["eeg", "eog"])
print(len(raw.ch_names), "→", len(eeg_and_eog.ch_names))

376 → 61

此外,pick()也可用于按名称拾取通道,相应的drop_channels()方法可用于按名称删除通道:

raw_temp = raw.copy()
print("Number of channels in raw_temp:")
print(len(raw_temp.ch_names), end=" → drop two → ")
raw_temp.drop_channels(["EEG 037", "EEG 059"])
print(len(raw_temp.ch_names), end=" → pick three → ")
raw_temp.pick(["MEG 1811", "EEG 017", "EOG 061"])
print(len(raw_temp.ch_names))

Number of channels in raw_temp:
376 → drop two → 374 → pick three → 3

如果你想要通道以特定的顺序(例如,绘图时),reorder channels()就像pick()一样工作,但也会对通道进行重新排序;例如,这里我们选择EOG和额叶EEG通道,将EOG放在首位,EEG以相反的顺序排列:

channel_names = ["EOG 061", "EEG 003", "EEG 002", "EEG 001"]
eog_and_frontal_eeg = raw.copy().reorder_channels(channel_names)
print(eog_and_frontal_eeg.ch_names)

[‘EOG 061’, ‘EEG 003’, ‘EEG 002’, ‘EEG 001’]

更改通道的名称和类型

可以使用rename_channels()方法重命名通道,该方法使用Python字典将旧名称映射为新名称。您不需要一次重命名所有通道;仅为要重命名的通道提供字典项。这里有一个例子:

raw.rename_channels({"EOG 061": "blink detector"})

由于.fif文件的限制,通道名称最多限制为15个字符。

下一个示例使用Python字典推导式将通道名称中的空格替换为下划线:

print(raw.ch_names[-3:])
channel_renaming_dict = {name: name.replace(" ", "_") for name in raw.ch_names}
raw.rename_channels(channel_renaming_dict)
print(raw.ch_names[-3:])

[‘EEG 059’, ‘EEG 060’, ‘blink detector’]
[‘EEG_059’, ‘EEG_060’, ‘blink_detector’]

如果由于某种原因,Raw对象中的通道类型不准确,则可以使用set_channel_types()方法更改任何通道的类型。该方法接受一个将通道名称映射到类型的字典.改变通道类型的一个常见用例是使用额叶EEG电极作为临时EOG通道。

raw.set_channel_types({"EEG_001": "eog"})
print(raw.copy().pick(picks="eog").ch_names)

[‘EEG_001’, ‘blink_detector’]

时间阈的选择

如果您想限制Raw对象的时域,可以使用crop()方法,该方法可以就地修改Raw对象。crop()接受参数tmin和tmax,两者都以秒为单位:

raw_selection = raw.copy().crop(tmin=10, tmax=12.5)
print(raw_selection)

<Raw | sample_audvis_raw.fif, 376 x 1503 (2.5 s), ~3.3 MB, data not loaded>

Crop()还修改了first_samptimes属性,因此被裁剪对象的第一个样本现在对应于time = 0。因此,如果你想将原始选择从11秒重新裁剪到12.5秒,那么随后对crop()的调用应该得到tmin=1,而不是tmin=11),并不指定tmax以保留从tmin到对象末尾的所有内容:

print(raw_selection.times.min(), raw_selection.times.max())
raw_selection.crop(tmin=1)
print(raw_selection.times.min(), raw_selection.times.max())

0.0 2.500770084699155
0.0 1.5001290587975622

采样时间并不总是与请求的tmin或tmax值完全一致,这是裁剪文件的最大值与请求的tmax不完全匹配的原因。

raw_selection1 = raw.copy().crop(tmin=30, tmax=30.1)  # 0.1 seconds
raw_selection2 = raw.copy().crop(tmin=40, tmax=41.1)  # 1.1 seconds
raw_selection3 = raw.copy().crop(tmin=50, tmax=51.3)  # 1.3 seconds
raw_selection1.append([raw_selection2, raw_selection3])  # 2.5 seconds total
print(raw_selection1.times.min(), raw_selection1.times.max())

0.0 2.5041000049184614

在连接来自不同记录的Raw对象时要小心,特别是在保存时:append()只保留初始Raw对象。

从raw对象中提取数据

索引原始数组和索引NumPy数组有两种不同的方式。

  1. 除了请求的值之外,MNE-Python还返回一个与请求的样例对应的时间数组,以秒为单位。data数组和times数组作为元组的元素一起返回。
  2. 即使只请求单个时间采样或单个通道,数据数组也将始终是二维的。

从编号中提取数据

为了说明以上两点,让我们从第一个通道中选择几秒钟的数据:

sampling_freq = raw.info["sfreq"]
start_stop_seconds = np.array([11, 13])

start_sample, stop_sample = (start_stop_seconds * sampling_freq).astype(int)
channel_index = 0
raw_selection = raw[channel_index, start_sample:stop_sample]
print(raw_selection)

(array([[-3.85742192e-12, -3.85742192e-12, -9.64355481e-13, …,
2.89306644e-12, 3.85742192e-12, 3.85742192e-12]]), array([10.99872648, 11.00039144, 11.0020564 , …, 12.9933487 ,
12.99501366, 12.99667862]))

你可以看到它包含两个数组。这种数据和时间的组合使得绘制原始数据变得很容易,但注意,我们将数据数组转置,使每个通道都是一列而不是一行,以符合matplotlib在绘制二维y和一维x时的期望。

x = raw_selection[1]
y = raw_selection[0].T
plt.plot(x, y)

在这里插入图片描述

按名称提取通道

Raw对象也可以用通道的名称而不是它们的索引号来索引。你可以传递一个字符串来获取一个通道,也可以传递一个字符串列表来选择多个通道。与整数索引一样,这将返回一个(data_array, times_array)的元组,可以很容易地绘制。由于我们这次绘制的是2个通道,我们将为其中一个通道添加一个垂直偏移量,这样它就不会被绘制在另一个通道的正上方,更美观一些。

channel_names = ["MEG_0712", "MEG_1022"]
two_meg_chans = raw[channel_names, start_sample:stop_sample]
# 第一个通道设置偏移量,第二个不设置
y_offset = np.array([5e-11, 0])  # just enough to separate the channel traces
x = two_meg_chans[1]
y = two_meg_chans[0].T + y_offset
lines = plt.plot(x, y)
plt.legend(lines, channel_names)

在这里插入图片描述

按类型提取通道

有几种方法可以从raw对象中选择给定类型的所有通道。最安全的方法是使用mne.pick_types()来获得所需通道的整数索引,然后将这些索引与上面所示的方括号索引方法一起使用。pick_types()函数使用raw对象的Info属性来确定通道类型,并接受布尔或字符串参数来表示保留哪一种类型。meg参数默认为True,其他参数默认为False,因此为了获取EEG通道,我们传入EEG =Truemeg=False

eeg_channel_indices = mne.pick_types(raw.info, meg=False, eeg=True)
eeg_data, times = raw[eeg_channel_indices]
print(eeg_data.shape)

(58, 36038)

mne.pick_types()的一些参数既可以接受字符串参数,也可以接受布尔值。例如,meg参数可以取值maggradplanar1planar2来选择磁力仪、所有梯度仪或特定类型的梯度仪。

Raw.get_data() 方法

如果你只想要数据,而不是相应的时间数组,Raw对象有一个get_data()方法。在不指定参数的情况下,它将从所有通道中提取所有数据,存储在一个(n_channels, n_timepoints) 的NumPy数组中:

data = raw.get_data()
print(data.shape)

(376, 36038)

如果你想要时间数组,get_data()有一个可选的return_times参数:

data, times = raw.get_data(return_times=True)
print(data.shape)
print(times.shape)

(376, 36038)
(36038,)

get_data()方法还可以通过其pickstartstop参数来提取特定的通道和样本范围。pick参数接受通道索引通道名称通道类型,并保留请求的通道顺序作为其pick参数。

first_channel_data = raw.get_data(picks=0)
eeg_and_eog_data = raw.get_data(picks=["eeg", "eog"])
two_meg_chans_data = raw.get_data(picks=["MEG_0712", "MEG_1022"], start=1000, stop=2000)

print(first_channel_data.shape)
print(eeg_and_eog_data.shape)
print(two_meg_chans_data.shape)

(1, 36038)
(61, 36038)
(2, 1000)

从raw对象提取数据的总结

在这里插入图片描述

导出和保存raw对象

raw对象有一个内置的save()方法,可以使用该方法将经过部分处理的raw对象作为.fif文件写入磁盘,以便稍后可以无损地重新加载它。
还有其他一些方法可以从原始对象中只导出传感器数据。一种是使用索引或get_data()方法提取数据,并使用numpy.save()保存数据数组:

data = raw.get_data()
np.save(file="my_data.npy", arr=data)

也可以将数据导出到Pandas_DataFrame对象中,并使用Pandas提供的保存方法。raw对象的to_data_frame()方法与get_data()类似,它有一个用于限制导出通道的picks参数,以及时间的startstop参数。默认情况下,times将被四舍五入到最接近的毫秒,并用作DataFrame的索引。

sampling_freq = raw.info["sfreq"]
start_end_secs = np.array([10, 13])
start_sample, stop_sample = (start_end_secs * sampling_freq).astype(int)
df = raw.to_data_frame(picks=["eeg"], start=start_sample, stop=stop_sample)
# then save using df.to_csv(...), df.to_hdf(...), etc
print(df.head())
    time       EEG_002  ...       EEG_059       EEG_060
    0   9.999750  3.847885e+08  ...  5.243790e+08  6.952283e+08
    1  10.001415  3.612900e+08  ...  5.200958e+08  7.069226e+08
    2  10.003080  3.589401e+08  ...  5.176483e+08  7.080921e+08
    3  10.004745  3.724518e+08  ...  5.182602e+08  7.010755e+08
    4  10.006410  3.847885e+08  ...  5.286621e+08  7.069226e+08
    
    [5 rows x 60 columns]
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值