用Python处理DEAP脑电数据的5个常见坑及解决方案（附代码示例）
·
从数据加载到特征工程：Python处理DEAP脑电数据的实战避坑指南
如果你刚开始接触DEAP数据集，可能会觉得它结构清晰、数据规整，上手应该不难。但真正动手用Python处理时，才会发现从.mat文件加载到最终特征提取的每一步都藏着不少“暗礁”。我见过不少研究生和算法工程师，在情绪识别项目初期，把大量时间浪费在了数据预处理的各种细节问题上，而不是模型设计本身。这篇文章就来自我处理DEAP数据时的真实笔记，记录了那些最容易让人栽跟头的环节，以及如何用Python代码优雅地跨过去。无论你是做情绪计算的研究人员，还是希望将生理信号应用于产品交互的开发者，这些经验都能帮你节省大量调试时间。
1. 数据加载与维度解析：避开第一个“维度陷阱”
拿到DEAP预处理数据（data_preprocessed_python文件夹）后，很多人会直接用scipy.io.loadmat加载.mat文件，然后看着data和labels这两个变量开始琢磨。这里第一个坑就出现了：数据维度的实际含义与直觉理解往往不一致。
官方文档说明data是40×40×8064的三维数组，但如果你直接打印shape，可能会得到(40, 40, 8064)，然后想当然地认为第一个40是试次，第二个40是通道，8064是时间点。这个理解基本正确，但实际操作时需要更小心。
1.1 正确加载与验证数据维度
首先，不要用默认的loadmat，因为MATLAB和Python在矩阵存储顺序上存在差异。MATLAB是列优先（Fortran顺序），而NumPy默认是行优先（C顺序）。虽然对于纯数据加载这可能不是大问题，但在后续处理中，不当的维度理解会导致特征提取完全错误。
import numpy as np
import scipy.io as sio
def load_deap_subject(file_path, subject_id):
"""
加载单个被试的DEAP预处理数据
参数：
file_path: 数据文件夹路径
subject_id: 被试编号（1-32）
返回：
data_array: 脑电数据，形状为(试次, 通道, 时间点)
labels_array: 标签数据，形状为(试次, 4)
channel_names: 通道名称列表（如果可用）
"""
# 构建文件名，注意DEAP文件命名从s01开始
file_name = f's{subject_id:02d}.mat'
full_path = os.path.join(file_path, file_name)
# 加载.mat文件，设置squeeze_me=True自动压缩单维度
mat_data = sio.loadmat(full_path, squeeze_me=True, struct_as_record=False)
# 提取数据和标签
data = mat_data['data']
labels = mat_data['labels']
# 验证维度
if data.ndim != 3:
raise ValueError(f"数据维度异常，期望3维，实际得到{data.ndim}维")
if data.shape[0] != 40 or data.shape[1] != 40:
print(f"警告：数据形状为{data.shape}，非标准40×40×8064")
return data, labels
注意：DEAP数据的前32个通道是EEG信号，后8个是外周生理信号（GSR、呼吸等）。如果你只关心脑电，需要明确区分。
加载后，我建议立即进行维度验证和基本统计：
# 加载第一个被试的数据
data, labels = load_deap_subject('data_preprocessed_python', 1)
print(f"数据形状: {data.shape}")
print(f"标签形状: {labels.shape}")
# 检查数据范围（EEG信号通常以微伏为单位）
print(f"数据范围: [{data.min():.2f}, {data.max():.2f}]")
print(f"数据均值: {data.mean():.4f} ± {data.std():.4f}")
# 查看前5个试次的标签（效价、唤醒度、优势度、喜欢程度）
print("\n前5个试次标签:")
for i in range(5):
print(f"试次{i+1}: 效价={labels[i, 0]:.2f}, "
f"唤醒度={labels[i, 1]:.2f}, "
f"优势度={labels[i, 2]:.2f}, "
f"喜欢程度={labels[i, 3]:.2f}")
1.2 通道选择与重参考处理
DEAP数据已经进行了预处理，包括降采样到128Hz、带通滤波和眼电伪迹去除。但原始数据是使用Biosemi系统采集的，采用的是平均参考。在某些分析中，你可能需要重新参考，比如转换为乳突参考或全脑平均参考。
这里有个细节：DEAP的32个EEG通道是按照国际10-20系统放置的，但顺序可能与你熟悉的顺序不同。我强烈建议创建一个通道位置映射表：
| 通道索引 | 通道名称 | 脑区 | 10-20系统位置 |
|---|---|---|---|
| 0 | Fp1 | 前额 | 左前额 |
| 1 | AF3 | 前额 | 左前额中部 |
| 2 | F3 | 额叶 | 左额 |
| 3 | F7 | 额叶 | 左额前部 |
| 4 | FC5 | 额中央 | 左额中央 |
| 5 | FC1 | 额中央 | 左额中央中部 |
| 6 | C3 | 中央 | 左中央 |
| 7 | T7 | 颞叶 | 左颞 |
| 8 | CP5 | 中央顶 | 左中央后部 |
| 9 | CP1 | 中央顶 | 左中央后部中部 |
| 10 | P3 | 顶叶 | 左顶 |
| 11 | P7 | 顶叶 | 左顶后部 |
| 12 | PO3 | 枕顶 | 左枕顶 |
| 13 | O1 | 枕叶 | 左枕 |
| 14 | Oz | 枕叶 | 中线枕 |
| 15 | Pz | 顶叶 | 中线顶 |
| 16 | Fp2 | 前额 | 右前额 |
| 17 | AF4 | 前额 | 右前额中部 |
| 18 | Fz | 额叶 | 中线额 |
| 19 | F4 | 额叶 | 右额 |
| 20 | F8 | 额叶 | 右额前部 |
| 21 | FC6 | 额中央 | 右额中央 |
| 22 | FC2 | 额中央 | 右额中央中部 |
| 23 | Cz | 中央 | 中线中央 |
| 24 | C4 | 中央 | 右中央 |
| 25 | T8 | 颞叶 | 右颞 |
| 26 | CP6 | 中央顶 | 右中央后部 |
| 27 | CP2 | 中央顶 | 右中央后部中部 |
| 28 | P4 | 顶叶 | 右顶 |
| 29 | P8 | 顶叶 | 右顶后部 |
| 30 | PO4 | 枕顶 | 右枕顶 |
| 31 | O2 | 枕叶 | 右枕 |
如果你需要重新参考，这里有个实用函数：
def rereference_eeg(data, ref_channels=None, method='average'):
"""
对EEG数据进行重参考
参数：
data: 原始EEG数据，形状为(试次, 通道, 时间点)
ref_channels: 参考通道索引列表，如果为None则使用所有通道
method: 参考方法，'average'为平均参考，'specific'为特定通道参考
返回：
rereferenced_data: 重参考后的数据
"""
# 只处理前32个EEG通道
eeg_data = data[:, :32, :].copy()
if method == 'average':
# 计算所有通道的平均作为参考
if ref_channels is None:
ref_signal = eeg_data.mean(axis=1, keepdims=True)
else:
ref_signal = eeg_data[:, ref_channels, :].mean(axis=1, keepdims=True)
# 减去参考信号
rereferenced = eeg_data - ref_signal
elif method == 'specific':
if ref_channels is None:
raise ValueError("使用'specific'方法时需要指定ref_channels")
# 使用特定通道作为参考
ref_signal = eeg_data[:, ref_channels, :].mean(axis=1, keepdims=True)
rereferenced = eeg_data - ref_signal
# 将处理后的EEG数据放回原数据数组
result = data.copy()
result[:, :32, :] = rereferenced
return result
2. 采样率与时间轴处理：时间对齐的精确控制
DEAP的预处理数据已经降采样到128Hz，每个试次63秒，所以总时间点是63×128=8064。这个计算看起来简单，但在实际处理中，时间轴对齐问题经常被忽视，导致后续的事件相关电位（ERP）分析或时频分析出现相位偏差。
2.1 构建精确的时间轴
首先，我们需要明确每个试次的时间零点在哪里。在DEAP实验中，每个音乐视频片段前有3秒的基线记录，然后是60秒的视频观看，所以总共63秒。但预处理后的数据是否包含了完整的63秒？时间零点对应的是视频开始还是其他时刻？
根据我的经验，预处理数据通常包含完整的63秒，其中：
- 前3秒（0-3秒）是基线期
- 后60秒（3-63秒）是刺激期
构建时间轴时：
def create_time_axis(sampling_rate=128, total_duration=63, baseline_duration=3):
"""
创建与DEAP数据对应的时间轴
参数：
sampling_rate: 采样率（Hz）
total_duration: 总时长（秒）
baseline_duration: 基线期时长（秒）
返回：
time_axis: 时间点数组（秒）
baseline_mask: 基线期布尔掩码
stimulus_mask: 刺激期布尔掩码
"""
# 计算总样本数
n_samples = int(total_duration * sampling_rate)
# 创建时间轴（从0开始）
time_axis = np.linspace(0, total_duration, n_samples, endpoint=False)
# 创建基线期和刺激期掩码
baseline_mask = time_axis < baseline_duration
stimulus_mask = time_axis >= baseline_duration
return time_axis, baseline_mask, stimulus_mask
# 使用示例
time_axis, baseline_mask, stimulus_mask = create_time_axis()
print(f"时间轴长度: {len(time_axis)} 个点")
print(f"时间范围: {time_axis[0]:.3f} 到 {time_axis[-1]:.3f} 秒")
print(f"基线期样本数: {baseline_mask.sum()}")
print(f"刺激期样本数: {stimulus_mask.sum()}")
2.2 基线校正的正确方法
基线校正是EEG处理中的关键步骤，目的是去除每个试次开始前的直流偏移和慢波漂移。DEAP数据已经进行了一定的预处理，但根据你的分析需求，可能还需要重新进行基线校正。
常见的错误是使用整个试次的均值进行校正，这会导致信号失真。正确的方法应该是使用基线期（前3秒）的均值：
def apply_baseline_correction(data, baseline_mask):
"""
应用基线校正
参数：
data: EEG数据，形状为(试次, 通道, 时间点)
baseline_mask: 基线期布尔掩码，形状为(时间点,)
返回：
corrected_data: 基线校正后的数据
"""
# 确保数据是浮点型，避免整数运算问题
if data.dtype != np.float32 and data.dtype != np.float64:
data = data.astype(np.float32)
# 计算每个试次、每个通道的基线均值
# baseline_mask需要扩展维度以匹配数据形状
baseline_data = data[:, :, baseline_mask]
baseline_mean = baseline_data.mean(axis=2, keepdims=True)
# 应用校正：减去基线均值
corrected_data = data - baseline_mean
return corrected_data
# 实际应用
corrected_data = apply_baseline_correction(data[:, :32, :], baseline_mask)
# 验证校正效果
print("基线校正验证:")
print(f"校正前基线均值: {data[:, :32, baseline_mask].mean():.4f}")
print(f"校正后基线均值: {corrected_data[:, :, baseline_mask].mean():.4f}")
提示：对于某些高级分析（如时频分析），你可能需要在变换后的域进行基线校正，而不是在时域。这取决于你的具体分析流程。
2.3 处理不完整或异常试次
虽然DEAP数据质量较高，但在实际分析中，你可能会遇到某些试次数据异常的情况。我建议在预处理流程中加入质量检查步骤：
def check_trial_quality(data, threshold_uv=200, max_bad_channels=5):
"""
检查试次数据质量
参数：
data: 单个试次的数据，形状为(通道, 时间点)
threshold_uv: 振幅阈值（微伏），超过此值视为异常
max_bad_channels: 允许的最大坏通道数
返回：
is_good: 试次是否合格（布尔值）
bad_channels: 坏通道索引列表
metrics: 质量指标字典
"""
metrics = {}
# 1. 检查振幅范围
peak_to_peak = data.max(axis=1) - data.min(axis=1)
bad_channels_amplitude = np.where(peak_to_peak > threshold_uv)[0]
# 2. 检查方差（零方差或极低方差可能表示坏通道）
channel_variances = data.var(axis=1)
variance_threshold = channel_variances.mean() * 0.01 # 低于均值1%视为异常
bad_channels_variance = np.where(channel_variances < variance_threshold)[0]
# 3. 检查NaN值
nan_channels = np.where(np.isnan(data).any(axis=1))[0]
# 合并所有坏通道
all_bad_channels = np.unique(np.concatenate([
bad_channels_amplitude,
bad_channels_variance,
nan_channels
]))
# 计算指标
metrics['peak_to_peak_mean'] = peak_to_peak.mean()
metrics['peak_to_peak_std'] = peak_to_peak.std()
metrics['variance_mean'] = channel_variances.mean()
metrics['bad_channel_count'] = len(all_bad_channels)
metrics['bad_channel_indices'] = all_bad_channels.tolist()
# 判断试次是否合格
is_good = len(all_bad_channels) <= max_bad_channels
return is_good, all_bad_channels, metrics
# 批量检查所有试次
def batch_quality_check(data, subject_id):
"""
批量检查一个被试的所有试次
"""
results = []
bad_trials = []
for trial_idx in range(data.shape[0]):
trial_data = data[trial_idx, :32, :] # 只检查EEG通道
is_good, bad_channels, metrics = check_trial_quality(trial_data)
results.append({
'subject': subject_id,
'trial': trial_idx,
'is_good': is_good,
'bad_channel_count': metrics['bad_channel_count'],
'bad_channels': metrics['bad_channel_indices'],
'peak_to_peak_mean': metrics['peak_to_peak_mean']
})
if not is_good:
bad_trials.append(trial_idx)
return results, bad_trials
3. 标签处理与情绪空间映射：超越二分类思维
DEAP数据提供了四个维度的连续标签：效价（Valence）、唤醒度（Arousal）、优势度（Dominance）和喜欢程度（Liking）。大多数研究只关注效价和唤醒度，但处理这些连续标签时，有几个关键决策点会影响你的模型性能。
3.1 连续标签与离散分类的转换
最直接的方法是将连续标签（1-9）二值化，通常以5为分界点。但这里有个微妙之处：中性地带（4.5-5.5）的样本应该如何处理？
def discretize_labels(labels, threshold=5.0, remove_ambiguous=True, ambiguous_range=(4.5, 5.5)):
"""
将连续标签离散化为二分类标签
参数：
labels: 原始标签数组，形状为(试次, 4)
threshold: 分类阈值
remove_ambiguous: 是否去除模糊样本
ambiguous_range: 模糊范围
返回：
binary_labels: 二值标签字典
valid_indices: 有效样本索引
"""
# 提取效价和唤醒度
valence = labels[:, 0]
arousal = labels[:, 1]
# 初始化有效样本掩码
valid_mask = np.ones(len(labels), dtype=bool)
if remove_ambiguous:
# 去除效价和唤醒度都在模糊范围内的样本
valence_ambiguous = (valence >= ambiguous_range[0]) & (valence <= ambiguous_range[1])
arousal_ambiguous = (arousal >= ambiguous_range[0]) & (arousal <= ambiguous_range[1])
ambiguous_mask = valence_ambiguous & arousal_ambiguous
valid_mask = ~ambiguous_mask
print(f"去除模糊样本: {ambiguous_mask.sum()} 个")
# 应用有效样本掩码
valence_valid = valence[valid_mask]
arousal_valid = arousal[valid_mask]
valid_indices = np.where(valid_mask)[0]
# 二值化
valence_binary = (valence_valid > threshold).astype(int) # 1=高效价，0=低效价
arousal_binary = (arousal_valid > threshold).astype(int) # 1=高唤醒，0=低唤醒
# 创建四象限分类（效价×唤醒度）
# 0: 低效价低唤醒，1: 低效价高唤醒，2: 高效价低唤醒，3: 高效价高唤醒
quadrant_labels = valence_binary * 2 + arousal_binary
return {
'valence': valence_binary,
'arousal': arousal_binary,
'quadrant': quadrant_labels,
'indices': valid_indices
}
# 使用示例
binary_labels = discretize_labels(labels, remove_ambiguous=True)
print(f"有效样本数: {len(binary_labels['indices'])}")
print(f"高效价样本: {(binary_labels['valence'] == 1).sum()}")
print(f"高唤醒样本: {(binary_labels['arousal'] == 1).sum()}")
print(f"四象限分布:")
for i in range(4):
count = (binary_labels['quadrant'] == i).sum()
print(f" 象限{i}: {count} 个样本")
3.2 标签分布不平衡问题
DEAP数据在情绪空间中的分布并不均匀。高效价高唤醒的样本通常较多，而低效价低唤醒的样本较少。这种不平衡会影响分类器的训练。
我常用的解决策略包括：
- 重采样技术（过采样少数类或欠采样多数类）
- 类别权重调整（在损失函数中给少数类更高权重）
- 数据增强（对EEG信号进行轻微变换生成新样本）
这里提供一个基于类别权重的处理示例：
from sklearn.utils.class_weight import compute_class_weight
def compute_class_weights(labels):
"""
计算类别权重以处理不平衡数据
参数：
labels: 类别标签数组
返回：
class_weights: 类别权重字典
"""
classes = np.unique(labels)
weights = compute_class_weight('balanced', classes=classes, y=labels)
return dict(zip(classes, weights))
# 为四象限分类计算权重
quadrant_weights = compute_class_weights(binary_labels['quadrant'])
print("四象限分类权重:")
for quadrant, weight in quadrant_weights.items():
count = (binary_labels['quadrant'] == quadrant).sum()
print(f" 象限{quadrant}: {count}个样本，权重={weight:.3f}")
# 在PyTorch或TensorFlow中使用这些权重
# PyTorch示例：
# criterion = nn.CrossEntropyLoss(weight=torch.tensor(list(quadrant_weights.values())))
3.3 连续标签的回归任务
除了分类，你还可以将情绪识别视为回归任务，直接预测效价和唤醒度的连续值。这时需要注意标签的标准化：
def prepare_regression_labels(labels, valid_indices=None):
"""
准备回归任务的标签
参数：
labels: 原始标签数组
valid_indices: 有效样本索引（如果已进行过滤）
返回：
regression_labels: 标准化后的回归标签
scaler: 用于后续逆变换的标准化器
"""
from sklearn.preprocessing import StandardScaler
if valid_indices is not None:
labels = labels[valid_indices]
# 提取效价和唤醒度
valence_arousal = labels[:, :2]
# 标准化（均值为0，方差为1）
scaler = StandardScaler()
normalized_labels = scaler.fit_transform(valence_arousal)
return normalized_labels, scaler
# 使用示例
regression_labels, label_scaler = prepare_regression_labels(labels, binary_labels['indices'])
print(f"回归标签形状: {regression_labels.shape}")
print(f"效价均值（标准化后）: {regression_labels[:, 0].mean():.4f}")
print(f"唤醒度均值（标准化后）: {regression_labels[:, 1].mean():.4f}")
4. 特征工程与提取：从时域到时频域的实用技巧
特征提取是情绪识别中最关键的环节之一。DEAP数据的时间分辨率（128Hz）和时长（63秒）为多种特征提取方法提供了可能。以下是几种我实践中发现有效的特征提取策略。
4.1 时域特征提取
时域特征计算简单，对计算资源要求低，适合作为基线特征：
def extract_time_domain_features(data, window_size=128, overlap=0.5):
"""
提取时域特征
参数：
data: EEG数据，形状为(试次, 通道, 时间点)
window_size: 滑动窗口大小（样本数）
overlap: 窗口重叠比例
返回：
features: 时域特征数组
feature_names: 特征名称列表
"""
n_trials, n_channels, n_samples = data.shape
# 计算滑动窗口参数
step_size = int(window_size * (1 - overlap))
n_windows = (n_samples - window_size) // step_size + 1
# 预分配特征数组
n_features_per_channel = 6 # 均值、方差、偏度、峰度、Hjorth参数等
features = np.zeros((n_trials, n_channels * n_features_per_channel, n_windows))
feature_names = []
for ch in range(n_channels):
for w in range(n_windows):
start_idx = w * step_size
end_idx = start_idx + window_size
window_data = data[:, ch, start_idx:end_idx]
# 1. 均值
mean_val = window_data.mean(axis=1)
# 2. 方差
var_val = window_data.var(axis=1)
# 3. 偏度（三阶中心矩）
from scipy.stats import skew
skewness = skew(window_data, axis=1)
# 4. 峰度（四阶中心矩）
from scipy.stats import kurtosis
kurt = kurtosis(window_data, axis=1)
# 5. Hjorth活动性（方差）
activity = var_val
# 6. Hjorth移动性（一阶导数方差的平方根除以活动性）
diff_signal = np.diff(window_data, axis=1)
mobility = np.sqrt(diff_signal.var(axis=1) / (var_val + 1e-10))
# 存储特征
base_idx = ch * n_features_per_channel
features[:, base_idx, w] = mean_val
features[:, base_idx + 1, w] = var_val
features[:, base_idx + 2, w] = skewness
features[:, base_idx + 3, w] = kurt
features[:, base_idx + 4, w] = activity
features[:, base_idx + 5, w] = mobility
# 生成特征名称
channel_name = f'Ch{ch:02d}'
feature_names.extend([
f'{channel_name}_mean',
f'{channel_name}_var',
f'{channel_name}_skew',
f'{channel_name}_kurt',
f'{channel_name}_activity',
f'{channel_name}_mobility'
])
# 将特征重塑为(试次, 特征×窗口)
features_reshaped = features.reshape(n_trials, -1)
return features_reshaped, feature_names
# 提取时域特征示例
time_features, time_feature_names = extract_time_domain_features(data[:, :32, :])
print(f"时域特征形状: {time_features.shape}")
print(f"特征数量: {len(time_feature_names)}")
print(f"前5个特征名: {time_feature_names[:5]}")
4.2 频域特征提取（功率谱密度）
频域特征对情绪识别特别重要，因为不同情绪状态与特定频带（α、β、θ、γ波）的功率变化相关：
def extract_frequency_domain_features(data, sampling_rate=128, band_definitions=None):
"""
提取频域特征（频带功率）
参数：
data: EEG数据，形状为(试次, 通道, 时间点)
sampling_rate: 采样率
band_definitions: 频带定义字典
返回：
band_powers: 频带功率特征
band_names: 频带名称列表
"""
from scipy.signal import welch
if band_definitions is None:
# 默认EEG频带定义
band_definitions = {
'delta': (1, 4),
'theta': (4, 8),
'alpha': (8, 13),
'beta': (13, 30),
'gamma': (30, 45)
}
n_trials, n_channels, n_samples = data.shape
n_bands = len(band_definitions)
# 预分配特征数组
band_powers = np.zeros((n_trials, n_channels * n_bands))
band_names = []
for ch in range(n_channels):
for trial_idx in range(n_trials):
# 计算功率谱密度
frequencies, psd = welch(
data[trial_idx, ch, :],
fs=sampling_rate,
nperseg=min(256, n_samples // 4),
noverlap=min(128, n_samples // 8)
)
# 计算各频带功率
for band_idx, (band_name, (low_freq, high_freq)) in enumerate(band_definitions.items()):
# 找到频带对应的频率索引
freq_mask = (frequencies >= low_freq) & (frequencies <= high_freq)
if freq_mask.any():
# 计算频带内功率的平均值（对数变换）
band_power = np.log(np.mean(psd[freq_mask]) + 1e-10)
feature_idx = ch * n_bands + band_idx
band_powers[trial_idx, feature_idx] = band_power
# 生成特征名称
for band_name in band_definitions.keys():
band_names.append(f'Ch{ch:02d}_{band_name}_power')
return band_powers, band_names
# 提取频域特征示例
freq_features, freq_feature_names = extract_frequency_domain_features(data[:, :32, :])
print(f"频域特征形状: {freq_features.shape}")
print(f"频带功率特征示例（第一个试次，前5个特征）:")
for i in range(5):
print(f" {freq_feature_names[i]}: {freq_features[0, i]:.4f}")
4.3 时频特征提取（小波变换）
时频分析能同时捕捉频率成分随时间的变化，这对非平稳的EEG信号特别有用：
def extract_wavelet_features(data, sampling_rate=128, wavelet='morl', scales=None):
"""
使用连续小波变换提取时频特征
参数：
data: EEG数据，形状为(试次, 通道, 时间点)
sampling_rate: 采样率
wavelet: 小波类型
scales: 尺度参数，如果为None则自动生成
返回：
wavelet_features: 小波特征
feature_info: 特征信息
"""
import pywt
from scipy import signal
n_trials, n_channels, n_samples = data.shape
if scales is None:
# 自动生成尺度，对应1-45Hz频率范围
# 小波中心频率与尺度的关系：f = fc * fs / scale
# 对于morlet小波，fc ≈ 0.8125
fc = pywt.central_frequency(wavelet)
# 计算对应1Hz的尺度
scale_low = fc * sampling_rate / 45 # 45Hz
scale_high = fc * sampling_rate / 1 # 1Hz
num_scales = 20
scales = np.linspace(scale_low, scale_high, num_scales)
# 预分配特征数组
# 我们取每个尺度在每个时间点的能量均值作为特征
n_scales = len(scales)
wavelet_features = np.zeros((n_trials, n_channels * n_scales))
for ch in range(n_channels):
print(f"处理通道 {ch+1}/{n_channels}", end='\r')
for trial_idx in range(n_trials):
signal_data = data[trial_idx, ch, :]
# 执行连续小波变换
coefficients, frequencies = pywt.cwt(
signal_data,
scales,
wavelet,
sampling_period=1.0/sampling_rate
)
# 计算每个尺度的平均能量
for scale_idx in range(n_scales):
energy = np.mean(np.abs(coefficients[scale_idx])**2)
feature_idx = ch * n_scales + scale_idx
wavelet_features[trial_idx, feature_idx] = np.log(energy + 1e-10)
print() # 换行
# 生成特征名称
feature_names = []
for ch in range(n_channels):
for scale_idx in range(n_scales):
# 计算该尺度对应的中心频率
fc = pywt.central_frequency(wavelet)
center_freq = fc * sampling_rate / scales[scale_idx]
feature_names.append(f'Ch{ch:02d}_wavelet_{center_freq:.1f}Hz')
return wavelet_features, {'scales': scales, 'feature_names': feature_names}
# 注意：小波变换计算量较大，建议先在小样本上测试
# 示例（注释掉以避免长时间运行）：
# wavelet_features, wavelet_info = extract_wavelet_features(data[:5, :5, :]) # 只处理前5个试次和前5个通道
4.4 特征选择与降维
提取大量特征后，下一步是选择最相关的特征。我通常采用多阶段特征选择策略：
def select_features_with_mrmr(features, labels, n_features_to_select=50):
"""
使用mRMR（最大相关最小冗余）进行特征选择
参数：
features: 特征矩阵，形状为(样本数, 特征数)
labels: 标签数组
n_features_to_select: 要选择的特征数量
返回：
selected_features: 选择后的特征
selected_indices: 选择的特征索引
feature_scores: 特征得分
"""
# 注意：这里简化实现，实际应用中可能需要安装mrmr_selection包
# 或者使用其他特征选择方法
from sklearn.feature_selection import SelectKBest, f_classif, mutual_info_classif
# 方法1：基于方差过滤（去除低方差特征）
from sklearn.feature_selection import VarianceThreshold
variance_selector = VarianceThreshold(threshold=0.01)
features_variance_filtered = variance_selector.fit_transform(features)
# 方法2：基于单变量统计检验
# 对于分类任务
if len(np.unique(labels)) < 10: # 假设是分类任务
selector = SelectKBest(score_func=f_classif, k=min(n_features_to_select*2, features_variance_filtered.shape[1]))
features_selected = selector.fit_transform(features_variance_filtered, labels)
selected_indices_step2 = selector.get_support(indices=True)
# 获取原始特征索引
variance_indices = np.where(variance_selector.get_support())[0]
selected_indices = variance_indices[selected_indices_step2]
# 方法3：基于互信息（对非线性关系更敏感）
mi_scores = mutual_info_classif(features, labels, random_state=42)
top_mi_indices = np.argsort(mi_scores)[-n_features_to_select:]
# 综合选择：取前两种方法的交集或并集
# 这里简单取并集
final_indices = np.unique(np.concatenate([selected_indices[:n_features_to_select//2], top_mi_indices]))
# 如果特征数仍过多，进一步筛选
if len(final_indices) > n_features_to_select:
# 基于互信息得分排序
final_indices_scores = mi_scores[final_indices]
top_indices = final_indices[np.argsort(final_indices_scores)[-n_features_to_select:]]
final_indices = top_indices
selected_features = features[:, final_indices]
return selected_features, final_indices, mi_scores
# 特征选择示例（假设是分类任务）
# 注意：这里需要实际的标签数据
# selected_features, selected_idx, feature_scores = select_features_with_mrmr(
# time_features,
# binary_labels['valence'],
# n_features_to_select=100
# )
5. 与EEGLAB的协同工作流：发挥工具链的最大效能
虽然Python在机器学习和深度学习方面有优势，但EEGLAB在EEG预处理和可视化方面仍然是行业标准。建立Python与EEGLAB的协同工作流可以发挥两者的优势。
5.1 数据格式转换：Python与EEGLAB/MATLAB互通
经常需要在Python和MATLAB/EEGLAB之间传递数据。以下是一些实用函数：
def save_for_eeglab(data, labels, channel_names, file_path):
"""
将Python处理的数据保存为EEGLAB可读格式
参数：
data: EEG数据，形状为(试次, 通道, 时间点)
labels: 标签数据
channel_names: 通道名称列表
file_path: 保存路径
"""
import h5py
# 将数据重塑为EEGLAB期望的格式：(通道, 时间点, 试次)
# EEGLAB通常使用(通道, 时间点)的2D数组，多个试次可以保存在结构中
eeglab_format = np.transpose(data, (1, 2, 0))
# 保存为.mat文件（使用v7.3格式以支持h5py，兼容性更好）
with h5py.File(file_path, 'w') as f:
# 保存EEG数据
f.create_dataset('EEG/data', data=eeglab_format)
# 保存通道信息
f.create_dataset('EEG/chanlocs/labels',
data=[name.encode('utf-8') for name in channel_names])
# 保存采样率
f.create_dataset('EEG/srate', data=128)
# 保存试次信息
f.create_dataset('EEG/trials', data=data.shape[0])
# 保存时间点信息
f.create_dataset('EEG/pnts', data=data.shape[2])
# 保存事件信息（将标签作为事件类型）
events = []
for trial_idx in range(data.shape[0]):
# 创建事件结构
event = {
'type': f'valence_{labels[trial_idx, 0]:.1f}_arousal_{labels[trial_idx, 1]:.1f}',
'latency': trial_idx * data.shape[2] + 1, # 试次开始位置
'duration': data.shape[2]
}
events.append(event)
# 注意：这里简化处理，实际需要更复杂的事件结构保存
print(f"数据已保存到 {file_path}，可在EEGLAB中使用load命令加载")
def load_eeglab_set(file_path):
"""
加载EEGLAB的.set文件（通过h5py）
参数：
file_path: .set文件路径
返回：
eeg_data: EEG数据数组
metadata: 元数据字典
"""
import h5py
with h5py.File(file_path, 'r') as f:
# EEGLAB的.set文件是HDF5格式
# 数据结构可能因EEGLAB版本而异
# 尝试查找数据
if 'EEG/data' in f:
data = f['EEG/data'][:]
# 转换为Python标准格式：(试次, 通道, 时间点)
if data.ndim == 3:
data = np.transpose(data, (2, 0, 1))
elif data.ndim == 2:
# 单试次数据
data = data[np.newaxis, :, :]
# 提取元数据
metadata = {}
if 'EEG/srate' in f:
metadata['srate'] = f['EEG/srate'][()]
if 'EEG/chanlocs/labels' in f:
labels = f['EEG/chanlocs/labels'][:]
metadata['channels'] = [label.decode('utf-8') for label in labels]
return data, metadata
5.2 利用EEGLAB进行ICA去伪迹
独立成分分析（ICA）是去除眼电、心电等伪迹的有效方法。虽然Python有ICA实现，但EEGLAB的ICA工具更成熟且经过充分验证。我通常的工作流是：
- 在Python中进行基本预处理（滤波、重参考）
- 将数据导出到EEGLAB进行ICA
- 识别并去除伪迹成分
- 将干净数据导回Python进行后续分析
这里提供一个自动化脚本的框架：
def run_eeglab_ica_automation(input_path, output_path, matlab_engine=None):
"""
通过MATLAB引擎调用EEGLAB进行ICA处理
参数：
input_path: 输入数据路径
output_path: 输出数据路径
matlab_engine: MATLAB引擎实例
返回：
处理状态和消息
"""
if matlab_engine is None:
# 尝试启动MATLAB引擎
try:
import matlab.engine
eng = matlab.engine.start_matlab()
except ImportError:
print("未找到MATLAB引擎，请安装MATLAB Engine API for Python")
return False, "MATLAB引擎不可用"
else:
eng = matlab_engine
# 添加EEGLAB到MATLAB路径
eeglab_path = '/path/to/eeglab' # 需要修改为实际路径
eng.addpath(eeglab_path, nargout=0)
# 构建MATLAB命令
matlab_script = f"""
% 加载数据
EEG = pop_loadset('filename', '{input_path}');
% 运行ICA
EEG = pop_runica(EEG, 'icatype', 'runica', 'extended', 1);
% 自动识别眼电成分（使用ADJUST插件）
% 需要先安装ADJUST插件
% EEG = pop_adjust(EEG);
% 或者手动选择成分后去除
% 这里假设我们已经知道要去除的成分索引
% bad_components = [1, 3, 5]; % 示例
% 去除成分
% EEG = pop_subcomp(EEG, bad_components, 0);
% 保存结果
pop_saveset(EEG, 'filename', '{output_path}');
"""
# 执行MATLAB脚本
try:
eng.eval(matlab_script, nargout=0)
return True, "ICA处理完成"
except Exception as e:
return False, f"处理失败: {str(e)}"
5.3 可视化与结果验证
可视化是验证预处理效果的关键。以下是一些有用的可视化函数：
def plot_eeg_trial(trial_data, channel_names, trial_idx=0, start_time=0, duration=5,
sampling_rate=128, figsize=(15, 10)):
"""
绘制单个试次的EEG波形
参数：
trial_data: 试次数据，形状为(试次, 通道, 时间点)
channel_names: 通道名称列表
trial_idx: 试次索引
start_time: 起始时间（秒）
duration: 显示时长（秒）
sampling_rate: 采样率
figsize: 图形大小
"""
import matplotlib.pyplot as plt
# 计算时间轴
n_samples = int(duration * sampling_rate)
start_sample = int(start_time * sampling_rate)
end_sample = start_sample + n_samples
# 提取数据
plot_data = trial_data[trial_idx, :, start_sample:end_sample]
time_axis = np.linspace(start_time, start_time + duration, n_samples, endpoint=False)
# 创建图形
fig, axes = plt.subplots(8, 4, figsize=figsize, sharex=True)
axes = axes.flatten()
# 绘制每个通道
for ch_idx in range(min(32, len(channel_names))):
ax = axes[ch_idx]
ax.plot(time_axis, plot_data[ch_idx, :], linewidth=0.8)
ax.set_ylabel(channel_names[ch_idx], fontsize=8)
ax.tick_params(labelsize=7)
# 设置y轴范围（基于数据范围）
y_min = plot_data[ch_idx, :].min()
y_max = plot_data[ch_idx, :].max()
y_padding = (y_max - y_min) * 0.1
ax.set_ylim(y_min - y_padding, y_max + y_padding)
# 设置公共标签
for i in range(32, len(axes)):
axes[i].axis('off')
fig.text(0.5, 0.04, '时间 (秒)', ha='center', fontsize=12)
fig.text(0.04, 0.5, '振幅 (μV)', va='center', rotation='vertical', fontsize=12)
fig.suptitle(f'试次 {trial_idx+1} 的EEG波形 ({start_time}-{start_time+duration}秒)', fontsize=14)
plt.tight_layout()
return fig
def plot_topographic_map(data, channel_names, trial_idx=0, time_point=None,
sampling_rate=128, figsize=(10, 8)):
"""
绘制地形图
参数：
data: EEG数据
channel_names: 通道名称列表
trial_idx: 试次索引
time_point: 时间点（秒），如果为None则使用平均值
sampling_rate: 采样率
figsize: 图形大小
"""
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# 通道位置（简化版，实际需要准确的3D坐标）
# 这里使用2D近似投影
channel_positions_2d = {
'Fp1': (-0.3, 0.8), 'Fp2': (0.3, 0.8),
'F3': (-0.5, 0.5), 'F4': (0.5, 0.5), 'Fz': (0, 0.6),
'C3': (-0.6, 0), 'C4': (0.6, 0), 'Cz': (0, 0.2),
'P3': (-0.5, -0.5), 'P4': (0.5, -0.5), 'Pz': (0, -0.3),
'O1': (-0.3, -0.8), 'O2': (0.3, -0.8), 'Oz': (0, -0.7)
# 可以添加更多通道...
}
# 提取数据
if time_point is None:
# 使用整个试次的平均值
channel_values = data[trial_idx, :32, :].mean(axis=1)
else:
# 使用特定时间点
sample_idx = int(time_point * sampling_rate)
channel_values = data[trial_idx, :32, sample_idx]
# 创建插值网格
grid_x, grid_y = np.mgrid[-1:1:100j, -1:1:100j]
# 准备已知点数据
points = []
values = []
for ch_idx, ch_name in enumerate(channel_names[:32]):
if ch_name in channel_positions_2d:
x, y = channel_positions_2d[ch_name]
points.append([x, y])
values.append(channel_values[ch_idx])
points = np.array(points)
values = np.array(values)
# 插值
grid_z = griddata(points, values, (grid_x, grid_y), method='cubic', fill_value=0)
# 绘制
fig, ax = plt.subplots(figsize=figsize)
contour = ax.contourf(grid_x, grid_y, grid_z, levels=50, cmap='RdBu_r')
# 添加通道位置点
for ch_name, (x, y) in channel_positions_2d.items():
ax.plot(x, y, 'ko', markersize=6)
ax.text(x, y + 0.05, ch_name, ha='center', fontsize=8)
# 添加颜色条
plt.colorbar(contour, ax=ax, label='振幅 (μV)')
# 设置标题
if time_point is None:
title = f'试次 {trial_idx+1} 的平均地形图'
else:
title = f'试次 {trial_idx+1} 在 {time_point:.2f} 秒的地形图'
ax.set_title(title)
ax.set_aspect('equal')
ax.axis('off')
return fig
5.4 性能评估与模型验证
最后，建立可靠的评估流程至关重要。对于情绪识别任务，我建议使用分层交叉验证，并考虑被试独立性：
def evaluate_emotion_classification(features, labels, subject_ids=None, n_splits=5):
"""
评估情绪分类性能
参数：
features: 特征矩阵
labels: 标签数组
subject_ids: 被试ID数组（用于被试独立的交叉验证）
n_splits: 交叉验证折数
返回：
results: 评估结果字典
"""
from sklearn.model_selection import StratifiedKFold, GroupKFold, cross_validate
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, f1_score, confusion_matrix
# 创建评估管道
pipelines = {
'SVM': Pipeline([
('scaler', StandardScaler()),
('classifier', SVC(kernel='rbf', C=1.0, gamma='scale', random_state=42))
]),
'RandomForest': Pipeline([
('scaler', StandardScaler()),
('classifier', RandomForestClassifier(n_estimators=100, random_state=42))
])
}
# 定义评估指标
scoring = {
'accuracy': 'accuracy',
'f1_weighted': 'f1_weighted',
'precision': 'precision_weighted',
'recall': 'recall_weighted'
}
results = {}
for model_name, pipeline in pipelines.items():
print(f"评估 {model_name}...")
if subject_ids is not None:
# 使用GroupKFold确保同一被试的数据不在训练集和测试集同时出现
cv = GroupKFold(n_splits=min(n_splits, len(np.unique(subject_ids))))
cv_splits = cv.split(features, labels, groups=subject_ids)
else:
# 使用分层K折交叉验证
cv = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=42)
cv_splits = cv.split(features, labels)
# 执行交叉验证
cv_results = cross_validate(
pipeline, features, labels,
cv=cv_splits,
scoring=scoring,
return_train_score=True,
n_jobs=-1
)
# 收集结果
results[model_name] = {
'test_accuracy': cv_results['test_accuracy'],
'test_f1_weighted': cv_results['test_f1_weighted'],
'train_accuracy': cv_results['train_accuracy'],
'fit_time': cv_results['fit_time'],
'score_time': cv_results['score_time']
}
# 打印平均性能
print(f" {model_name} 平均准确率: {cv_results['test_accuracy'].mean():.3f} ± {cv_results['test_accuracy'].std():.3f}")
print(f" {model_name} 平均F1分数: {cv_results['test_f1_weighted'].mean():.3f} ± {cv_results['test_f1_weighted'].std():.3f}")
return results
def plot_confusion_matrix(y_true, y_pred, class_names=None, figsize=(8, 6)):
"""
绘制混淆矩阵
参数：
y_true: 真实标签
y_pred: 预测标签
class_names: 类别名称列表
figsize: 图形大小
"""
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(y_true, y_pred)
if class_names is None:
class_names = [f'Class {i}' for i in range(cm.shape[0])]
fig, ax = plt.subplots(figsize=figsize)
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
xticklabels=class_names, yticklabels=class_names, ax=ax)
ax.set_xlabel('预测标签')
ax.set_ylabel('真实标签')
ax.set_title('混淆矩阵')
return fig
处理DEAP数据时，我最大的体会是：细节决定成败。那些看似微小的预处理选择——比如基线校正的时间窗口、重参考的方法、特征提取的参数——往往比模型架构的选择影响更大。特别是在情绪识别这种信噪比低的任务中，干净、一致的数据预处理流程是获得可靠结果的基石。上面这些代码片段都是我实际项目中提炼出来的，但每个项目都有其特殊性，你需要根据具体的研究问题和数据特点进行调整。最重要的是建立系统化的预处理流程，并记录每一个决策点，这样才能确保结果的可复现性。
Agent 垂直技术社区，欢迎活跃、内容共建。
更多推荐
5
0- 0
已为社区贡献2条内容
所有评论(0)