基于Python的脑电图（EEG）信号分析（5）

背景

脑电图（EEG）和运动想象任务

        脑电图（EEG）是一种通过在头皮上放置电极来记录大脑电活动的技术，广泛应用于神经科学研究和脑机接口（BCI）开发中。运动想象任务是一种常见的实验范式，其中参与者被要求想象或执行某种运动，如握拳或抬手。通过分析这些任务期间的EEG信号，我们可以深入理解大脑的运动控制机制和神经网络活动。

数据集介绍

        本文使用的运动想象数据集来自一个参与者，他们在实验中被要求（1）用左手和右手握球（实际执行条件）以及（2）想象用左手和右手握球（想象条件）。实验数据包括固定注视任务和运动想象任务的记录。

研究目标

        在本篇文章中，我们将探讨以下两个主要内容：

1. 事件相关去同步化（Event-Related Desynchronization, ERD）
2. 使用卷积神经网络（CNN）解码运动想象

代码实现

1. 数据加载与预处理

        首先，我们加载运动想象数据集并进行必要的预处理。

import pandas as pd
import numpy as np
import matplotlib.pyplot as pltimport mne
from mne import create_info
from mne.io import RawArray
from mne import Epochs, find_events
import osimport torch
import torch.nn as nn
from torch.utils.data import TensorDataset, WeightedRandomSampler# 配置绘图尺寸
mne.set_config('MNE_BROWSE_RAW_SIZE','10,5')  # 设置设备
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
torch.cuda.get_device_name(0)# 配置参数
base_url = '../data/S020'
runs = ['04', '06', '08']
in_channels  = 2  # C3, C4
out_channels = 5  # 捕获8, 9, 10, 11, 12 Hz频率
out_size     = 2  # 左或右
kernel_size  = 250 # 1秒
BATCH_SIZE   = 9999 # 使用整个批次大小
lr           = 1e-6
weight_decay = 1e-1
num_epochs = 1000
modelpath  = '../models/motorCNN1D.pth.tar'# 加载数据
raws = []
for run in runs:filename = f'S020R{run}_raw'path_file = os.path.join(base_url, filename + '.fif')raw = mne.io.read_raw_fif(path_file , preload=True, verbose='Warning')raws.append(raw)eeg = mne.io.concatenate_raws(raws)
df = eeg.to_data_frame()
print(df.shape)
df.head()# 只保留C3和C4通道
df = df[['C3', 'C4', 'STIM MARKERS']]
df.head()# 从原始数据中删除不需要的通道
eeg = eeg.drop_channels(['Fp1', 'Fp2', 'P7', 'P8', 'O1', 'O2'])
eeg.ch_names# 检查标签分布
np.unique(df['STIM MARKERS'], return_counts=True)

2. 伪影去除与滤波

        为了确保数据质量，我们对EEG信号进行伪影去除和频率滤波处理。

# 去除50Hz电源噪声
eeg.notch_filter(50)# 滤波以捕获运动想象频率范围
eeg.filter(7, 14)
eeg.compute_psd().plot()

3. 时域切片（Epoching）

        我们对实验数据进行时域切片，以提取特定时间窗口内的信号。

def getEpochs(raw, events, event_id, tmin, tmax, picks, baseline):epochs = Epochs(raw, events=events, event_id=event_id, tmin=tmin, tmax=tmax, baseline=baseline, preload=True, picks=picks)print('样本丢失 %: ', (1 - len(epochs.events)/len(events)) * 100)return epochs# 设置事件ID和时间窗口
event_id = {'Left': 1, 'Right' : 2}
events = find_events(eeg)
events_no_fixation = mne.pick_events(events, exclude=4)  # 忽略固定注视事件
tmin = 0
tmax = 7
baseline = (0, 2)  # 基线为事件前2秒
eeg_channels = mne.pick_types(eeg.info, eeg=True)
picks = eeg_channels
epochs = getEpochs(eeg, events, event_id, tmin, tmax, picks=picks, baseline=baseline)
epochs

4. 事件相关去同步化（ERD）分析

        事件相关去同步化（ERD）是指在特定频率带内，EEG信号功率的相对降低，通常在实际运动执行或想象过程中观察到。

from mne.time_frequency import tfr_multitaper# 提取不同频率带的功率
freqs = np.arange(7, 14)  # 运动想象的频率范围
tfr = tfr_multitaper(epochs, freqs=freqs, n_cycles=freqs, use_fft=True,return_itc=False, average=False, decim=1)# 进行基线校正，并转换为百分比
tfr.crop(tmin, tmax).apply_baseline(baseline, mode="percent")# 转换为DataFrame进行分析
df = tfr.to_data_frame(time_format=None, long_format=True)
df.head()
np.unique(df.freq, return_counts=True)# 为频率分配频段标签
freq_bounds = {'_': 0, 'delta': 3, 'theta': 7, 'alpha': 13, 'beta': 35, 'gamma': 140}
df['band'] = pd.cut(df['freq'], list(freq_bounds.values()), labels=list(freq_bounds)[1:])
df['band'] = df['band'].cat.remove_unused_categories() 
df['band'].cat.categories# 绘制ERD/ERS图
import seaborn as sns
g = sns.FacetGrid(df, row='band', col='condition')
g.map(sns.lineplot, 'time', 'value', 'channel')
axline_kw = dict(color='black', linestyle='dashed', linewidth=0.5, alpha=0.5)
g.map(plt.axhline, y=0, **axline_kw)  # 添加水平线
g.map(plt.axvline, x=0, **axline_kw)  # 添加垂直线
g.set(ylim=(-1, None))
g.set_axis_labels("Time (s)", "ERDS (%)")
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.add_legend(ncol=2, loc='lower center')
g.fig

5. 机器学习分析

        在进行端到端神经网络训练前，我们可以尝试基于特征工程的传统机器学习方法。机器学习方法通常能在适当的特征工程下取得竞争力的性能。

from sklearn.preprocessing import LabelEncoder# 编码通道标签
labelencoder = LabelEncoder()
df['channel'] = labelencoder.fit_transform(df['channel'])# 准备训练数据和标签
X = df[['value', 'freq', 'channel', 'time']]
y = df['condition']# 分割训练集和测试集
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
X_train.shape, X_test.shape, y_train.shape, y_test.shape# 简单拟合模型
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_scoreclf = RandomForestClassifier()
score = clf.fit(X_train, y_train)y_hat = clf.predict(X_test)
print("Test acc: ", accuracy_score(y_test, y_hat))# 特征重要性分析
feature_importance = pd.DataFrame(['value', 'freq', 'channel', 'time'], columns=['features'])
feature_importance['importance'] = clf.feature_importances_ 
feature_importance = feature_importance.sort_values(by = ['importance'], ascending=True)
feature_importance.plot.barh(x='features')
plt.show();

6. 深度学习模型

        接下来，我们将构建一个端到端的卷积神经网络（CNN）模型，直接从原始数据中提取特征进行分类。

# 获取数据和标签
X = epochs.get_data()
y = epochs.events[:, -1]# 数据标准化
from sklearn.preprocessing import MinMaxScalerscalers = {}
for i in range(X.shape[1]):scalers[i] = MinMaxScaler(feature_range=(0, 1))X[:, i, :] = scalers[i].fit_transform(X[:, i, :]) # 转换为PyTorch张量
X_train = torch.FloatTensor(X)
y_train = torch.LongTensor(y) - 1  # 标签从0开始# 创建数据加载器
ds_train = TensorDataset(X_train, y_train)
train_loader = torch.utils.data.DataLoader(dataset=ds_train, batch_size=BATCH_SIZE)# 定义CNN模型
class motorCNN1D(nn.Module):def __init__(self, in_channels, out_channels, kernel_size, out_size):super().__init__()self.conv1d = nn.Sequential(nn.Conv1d(in_channels, out_channels, kernel_size=kernel_size),nn.ELU(inplace=True),nn.AvgPool1d(10, stride=2),nn.BatchNorm1d(out_channels),nn.Dropout(0.50))self.fc = nn.Linear(linear_shape, out_size)def forward(self, seq):out = self.conv1d(seq)out = out.reshape(out.size(0), -1)out = self.fc(out)return outmodel = motorCNN1D(in_channels, out_channels, kernel_size, out_size).to(device)# 定义损失函数和优化器
criterion = nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(), lr = lr, weight_decay = lr)# 训练模型
best_valid_loss = float('inf')
model.train()for i in range(num_epochs):train_total = 0train_correct = 0train_acc   = 0for X_train, y_train in train_loader:X_train = X_train.to(device)y_train = y_train.to(device)yhat_train = model(X_train)# 计算训练准确率_, predicted = torch.max(yhat_train.data, 1)train_total += y_train.size(0)train_correct += (predicted == y_train).sum().item()train_acc = 100 * (train_correct / train_total)# 计算损失并更新模型train_loss = criterion(yhat_train, y_train)optimizer.zero_grad()train_loss.backward()optimizer.step()# 保存最佳模型if train_loss < best_valid_loss:best_valid_loss = train_losstorch.save(model.state_dict(), modelpath)if i % 100 == 0:print(f"Epoch: {i:4.0f} | Train acc: {train_acc: 3.2f} | " +f"loss: {train_loss:3.2f}")

7. 模型测试

        最后，我们使用测试集评估模型的性能。

model = motorCNN1D(in_channels, out_channels, kernel_size, out_size).to(device)
model.load_state_dict(torch.load(modelpath))
model.eval()with torch.no_grad():total = 0correct = 0for X_test, y_test in test_loader:X_test = X_test.to(device)y_test = y_test.to(device)predictions = model(X_test)_, predicted = torch.max(predictions.data, 1)total += y_test.size(0)correct += (predicted == y_test).sum().item()acc = 100 * (correct / total)print(f"Accuracy: {acc:2.3f}")

结语

本文结语

        在本文中，我们探讨了基于运动想象数据集的事件相关去同步化（ERD）分析和卷积神经网络（CNN）模型的构建与训练。通过深入的特征提取和信号处理，我们展示了如何从原始EEG数据中提取有意义的特征，并使用这些特征进行分类和预测。

系列结语

        通过这个系列的五篇文章，我们深入探讨了EEG信号分析的各个方面：

1. 基于时间、频率和成分的EEG信号分析：介绍了基础的EEG信号处理和分析方法。
2. 1维卷积神经网络（Conv1D）在EEG信号分析中的应用：展示了如何使用卷积神经网络处理和分析EEG数据。
3. EEG情感识别：使用DEAP数据集进行情感识别，并通过频谱分析和非对称性分析提取特征。
4. SSVEP数据集上的典型相关分析（CCA）：探索了如何使用CCA技术分析稳态视觉诱发电位数据。
5. 运动想象数据集上的事件相关去同步化（ERD）分析和深度学习模型：结合ERD和CNN模型，解码运动想象任务中的脑电信号。

        这些分析和方法不仅为我们理解大脑活动提供了新的视角，也为脑机接口的开发和应用奠定了坚实的基础。我们相信，通过不断的探索和创新，能够在脑科学和神经技术领域取得更大的突破。

        希望本系列文章能为读者提供有价值的参考，并激发更多关于脑电图信号分析的研究和应用。感谢您的阅读！

如果你觉得这篇博文对你有帮助，请点赞、收藏、关注我，并且可以打赏支持我！

欢迎关注我的后续博文，我将分享更多关于人工智能、自然语言处理和计算机视觉的精彩内容。

谢谢大家的支持！