简介本资源为2020年研究生数学建模竞赛C题「面向康复工程的脑电信号分析和判别模型」的完整代码与数据包面向参加电赛及数学建模竞赛的研究生、本科生以及从事脑电信号处理与生物医学工程方向的学习者。包内共263个文件以Python脚本与编译文件、PNG图表、XLSX数据表、XML配置、TXT说明及DOCX赛题文档为主压缩包约118.59MB涵盖P300脑机接口数据、睡眠分类等实验模块便于复现预处理、特征提取与建模分析全流程。已有45人学习下载。读者可据此获得赛题方案、算法实现与数据组织思路对照练习信号去噪、特征工程及SVM、神经网络等分类预测方法并参考报告撰写框架适合作为竞赛复盘与课程实践的参考素材。1. 脑电波分析赛题从一份 zip 到可复现的建模流水线2020 年研究生数学建模竞赛 C 题给的是脑电波数据赛题包里通常包含原始信号、受试者标签和一份说明文档。很多人拿到这个 zip 的第一反应是「先跑起来再说」结果卡在数据读取上——脑电信号不是普通表格采样率、通道数、事件标记三者对不上后面所有特征工程都是空中楼阁。这道题真正考的不是某个分类器调得多好而是你能不能把一段非平稳、低信噪比的时序信号拆成可解释、可复现的特征再落到一个能自圆其说的模型上。适合正在做数学建模、信号处理课程设计或者想入门生理信号分析的人。下面按「读数据 → 提特征 → 建模型 → 避坑」的顺序把这条流水线讲透代码可以直接抄。2. 脑电数据到底长什么样先搞清楚三个维度再动手2.1 采样率、通道数、事件标记的三角关系脑电数据的标准结构是一个二维矩阵加一条时间轴。矩阵的行是通道电极位置列是采样点。假设采样率是 250 Hz记录 60 秒那就是 250 × 60 15000 个采样点。通道数常见的是 14、32、64 导联赛题数据一般不会超过 64 通道。事件标记trigger是另一条独立的时间序列记录某个时刻发生了什么——比如第 3.2 秒出现刺激第 5.8 秒受试者按键。这三者的关系必须对齐事件标记的时间戳要换算成采样点索引才能从原始矩阵里切出对应片段。常见错误是拿秒数直接当索引去切片结果切出来的数据完全错位。换算公式很简单# 事件时间(秒) 转 采样点索引 sample_index int(event_time_sec * sfreq) # 例如 sfreq250, event_time_sec3.2 - sample_index800参数说明sfreq是采样率单位 Hz必须从数据说明文档里确认不能猜。event_time_sec是事件发生的绝对时间。int()取整是因为索引必须是整数但要注意取整方向——用round()还是int()取决于你对齐精度要求一般int()向下取整就够。2.2 用 Python 把 zip 里的数据读进来赛题包里的数据格式通常是.matMATLAB或.csv。.mat文件用scipy.io.loadmat读.csv用pandas读。先看数据规模再决定内存策略import scipy.io as sio import numpy as np # 读取 mat 文件结构通常是 dict data sio.loadmat(eeg_data.mat) print(data.keys()) # 先看有哪些变量名 # 假设信号变量叫 EEG形状是 (通道数, 采样点) eeg data[EEG] print(f通道数: {eeg.shape[0]}, 采样点数: {eeg.shape[1]}) # 假设采样率存在 fs 变量里 sfreq int(data[fs][0][0]) print(f采样率: {sfreq} Hz, 时长: {eeg.shape[1]/sfreq:.1f} 秒)逻辑说明loadmat返回的是一个字典键名取决于原始 MATLAB 脚本怎么存的。先print(data.keys())确认变量名不要硬编码。eeg.shape返回(通道数, 采样点数)如果读出来是反的说明数据存储时转置过需要手动.T。采样率经常存在一个 1×1 的数组里用[0][0]取标量。参数说明eeg的数据类型通常是float64如果内存吃紧可以转float32精度损失对后续滤波影响很小。sfreq必须和事件标记的时间基准一致如果事件标记用的是毫秒记得先除以 1000。2.3 数据质量检查三行代码排除 80% 的坑读进来之后别急着提特征先做三件事看有没有 NaN、看幅值范围、看通道间相关性。# 检查缺失值 print(fNaN 数量: {np.isnan(eeg).sum()}) # 检查幅值范围脑电典型幅值在 ±100 μV 以内 print(f幅值范围: [{eeg.min():.2f}, {eeg.max():.2f}] μV) # 检查通道间相关性正常脑电通道间有相关性但不会完全一致 corr np.corrcoef(eeg[:5, :]) # 取前5个通道 print(f通道间平均相关: {corr[np.triu_indices(5,1)].mean():.3f})如果 NaN 数量不为零要么是采集时丢包要么是某个通道坏了。幅值超过 ±200 μV 基本可以判定是伪迹眨眼、肌肉活动。通道间相关性如果接近 1.0说明两个电极可能贴在了同一个位置或者短路了。这三步做完你对数据的底子就有数了。3. 特征工程把一段波形变成模型能吃的数字3.1 滤波带通 陷波顺序不能反脑电信号的有效频段通常在 0.5–45 Hz。工频干扰50 Hz 或 60 Hz必须去掉否则功率谱上会有一根刺。滤波顺序有讲究先带通再陷波还是先陷波再带通我一般先带通把高频噪声整体压下去再用陷波精准打掉工频。反过来做的话陷波器在宽频噪声上的瞬态响应会更难处理。from scipy.signal import butter, filtfilt, iirnotch def bandpass_filter(data, lowcut, highcut, fs, order4): nyq 0.5 * fs b, a butter(order, [lowcut/nyq, highcut/nyq], btypeband) return filtfilt(b, a, data, axis1) def notch_filter(data, freq, fs, quality30): b, a iirnotch(freq, quality, fs) return filtfilt(b, a, data, axis1) # 先带通 0.5-45 Hz再陷波 50 Hz eeg_filtered bandpass_filter(eeg, 0.5, 45, sfreq) eeg_filtered notch_filter(eeg_filtered, 50, sfreq)逻辑说明filtfilt做的是零相位滤波前后各滤波一次避免相位偏移。butter的order4是常用值阶数太高会不稳定太低则过渡带太宽。iirnotch的quality30控制陷波带宽值越大陷波越窄一般 30–50 之间。参数说明lowcut和highcut根据你的分析目标调整。如果关注慢波delta 波低频可以放到 0.5 Hz如果只看高频振荡gamma 波高频可以放到 100 Hz 但采样率要够。axis1表示沿时间轴滤波不要搞错方向。3.2 分帧与加窗把连续信号切成模型能吃的片段脑电信号是非平稳的整段做 FFT 没有意义。标准做法是分帧每帧 1–2 秒帧间重叠 50%。加窗通常是汉宁窗是为了减少频谱泄漏。def segment_signal(data, sfreq, window_sec1.0, overlap0.5): window_len int(window_sec * sfreq) step int(window_len * (1 - overlap)) segments [] for start in range(0, data.shape[1] - window_len, step): seg data[:, start:start window_len] # 加汉宁窗 window np.hanning(window_len) seg seg * window segments.append(seg) return np.array(segments) # 形状: (帧数, 通道数, 窗长) segments segment_signal(eeg_filtered, sfreq, window_sec1.0, overlap0.5) print(f分帧后形状: {segments.shape})逻辑说明window_len是每帧的采样点数step是帧间步长。重叠 50% 意味着步长是窗长的一半。加窗时seg * window是逐点相乘窗函数两端趋近于零中间为一能有效抑制截断处的频谱泄漏。参数说明window_sec选 1 秒是脑电分析的常用值太短频率分辨率不够太长则时间分辨率下降。overlap取 0.5 是经验值增加到 0.75 可以增加样本量但计算量翻倍。3.3 功率谱密度与频带能量最稳的基线特征功率谱密度PSD是脑电分析里最基础也最稳的特征。用 Welch 方法估计 PSD然后按频带积分得到 delta0.5–4 Hz、theta4–8 Hz、alpha8–13 Hz、beta13–30 Hz、gamma30–45 Hz五个频带的能量。from scipy.signal import welch def band_power(segments, sfreq, bandsNone): if bands is None: bands {delta: (0.5, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30), gamma: (30, 45)} features [] for seg in segments: freqs, psd welch(seg, sfreq, npersegseg.shape[1], axis1) band_feats [] for name, (low, high) in bands.items(): idx np.logical_and(freqs low, freqs high) band_feats.append(np.trapz(psd[:, idx], freqs[idx], axis1)) features.append(np.concatenate(band_feats)) return np.array(features) # 形状: (帧数, 通道数*频带数) features band_power(segments, sfreq) print(f特征矩阵形状: {features.shape})逻辑说明welch返回频率轴和对应的功率谱密度。np.trapz做梯形积分把频带内的 PSD 累加成能量值。每个通道每个频带一个值所以特征维度是通道数 × 5。参数说明nperseg设为窗长表示不进一步分段直接用整帧做 FFT。如果帧长超过 2 秒可以设npersegsfreq*2来增加频率分辨率。np.trapz的axis1表示沿频率轴积分。3.4 时域特征补充Hjorth 参数与过零率频域特征之外时域特征计算快、解释性强。Hjorth 参数包括活动度方差、移动度一阶导数标准差与信号标准差之比、复杂度移动度的一阶导数与移动度之比。过零率反映信号在零线附近振荡的频率。def hjorth_params(segments): activity np.var(segments, axis2) diff1 np.diff(segments, axis2) mobility np.sqrt(np.var(diff1, axis2) / (activity 1e-10)) diff2 np.diff(diff1, axis2) mobility_d1 np.sqrt(np.var(diff2, axis2) / (np.var(diff1, axis2) 1e-10)) complexity mobility_d1 / (mobility 1e-10) return activity, mobility, complexity act, mob, comp hjorth_params(segments) print(f活动度形状: {act.shape}) # (帧数, 通道数)逻辑说明np.var沿时间轴算方差。np.diff算一阶差分近似导数。分母加1e-10防止除零。三个参数分别反映信号的幅度、频率和频率变化率。参数说明Hjorth 参数对窗长不敏感1 秒和 2 秒窗算出来的值差异不大。但如果信号幅值很小接近零1e-10的保护可能不够可以改成1e-6。4. 建模与验证从特征矩阵到分类结果4.1 分类器选型为什么我优先用 SVM 而不是深度学习赛题数据量通常不大——几十个受试者每人几十到几百帧。这种规模下SVM 配合 RBF 核往往比深度学习更稳。深度学习需要大量数据才能学到有意义的表示小样本上容易过拟合。随机森林也可以但特征维度高时解释性不如 SVM 的核技巧直观。from sklearn.svm import SVC from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.model_selection import cross_val_score, StratifiedKFold # 构建流水线标准化 SVM pipe Pipeline([ (scaler, StandardScaler()), (svm, SVC(kernelrbf, C1.0, gammascale)) ]) # 5 折交叉验证 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(pipe, features, labels, cvcv, scoringaccuracy) print(f交叉验证准确率: {scores.mean():.3f} ± {scores.std():.3f})逻辑说明StandardScaler把每个特征维度标准化到零均值单位方差SVM 对尺度敏感这一步不能省。gammascale自动根据特征数和方差设置核宽度比手动调省事。StratifiedKFold保证每折的类别比例一致。参数说明C是正则化参数越大越容易过拟合越小越容易欠拟合。gamma控制核函数的宽度scale是1/(n_features * X.var())的简写。如果交叉验证准确率波动大先检查特征里有没有异常值。4.2 特征重要性分析用排列重要性看哪些特征在起作用SVM 没有直接的特征重要性输出但可以用排列重要性permutation importance来评估。思路是打乱某个特征的值看模型性能下降多少下降越多说明该特征越重要。from sklearn.inspection import permutation_importance # 先在训练集上拟合 pipe.fit(X_train, y_train) # 计算排列重要性 result permutation_importance(pipe, X_test, y_test, n_repeats10, random_state42) importance result.importances_mean # 假设特征名按通道和频带排列 feature_names [fch{i}_{band} for i in range(n_channels) for band in [delta,theta,alpha,beta,gamma]] top_idx np.argsort(importance)[-10:] for i in top_idx: print(f{feature_names[i]}: {importance[i]:.4f})逻辑说明permutation_importance对每个特征重复打乱n_repeats次计算性能下降的均值和标准差。importances_mean是平均下降幅度越大越重要。参数说明n_repeats10是精度和计算量的折中特征多时可以降到 5。random_state固定随机种子保证可复现。4.3 交叉验证的坑别用普通 KFold 做分类分类问题必须用分层交叉验证StratifiedKFold否则某折可能全是某一类准确率直接崩。另外标准化必须在每折内部做不能先标准化再交叉验证——那叫数据泄漏。# 错误做法先标准化再交叉验证 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 泄漏了测试集信息 scores_wrong cross_val_score(SVC(), X_scaled, y, cv5) # 正确做法用 Pipeline 把标准化包进去 pipe Pipeline([(scaler, StandardScaler()), (svm, SVC())]) scores_right cross_val_score(pipe, X, y, cv5)逻辑说明fit_transform在全体数据上计算均值和方差测试集的信息泄漏到了训练过程。Pipeline保证每折的标准化只在该折的训练集上拟合。参数说明如果特征维度远大于样本数标准化后 SVM 可能仍然过拟合此时可以加 PCA 降维但 PCA 也要放进 Pipeline。5. 避坑与排查脑电分析里最容易翻车的五个地方5.1 现象交叉验证准确率 99%换一批数据掉到 50%原因特征里混入了标签泄漏。比如分帧时把同一受试者的帧同时分到了训练集和测试集模型记住了受试者的个体特征而不是任务相关特征。解决按受试者划分训练集和测试集而不是按帧随机划分。用GroupKFold指定受试者 ID 作为分组变量。5.2 现象滤波后信号幅值变得极小或极大原因filtfilt的滤波器阶数太高导致数值不稳定或者lowcut设得太低比如 0.1 Hz而数据长度不够。解决降低滤波器阶数到 2–4 阶或者改用sosfiltfilt二阶节形式提高数值稳定性。lowcut不要低于 0.5 Hz。5.3 现象PSD 图上 50 Hz 处仍然有尖峰原因陷波滤波器的quality参数设得太高陷波带宽太窄没有完全覆盖工频干扰的频谱扩散。解决把quality降到 20–30或者改用自适应滤波。如果工频干扰特别严重可以在带通之前先做一次陷波。5.4 现象SVM 训练报错「收敛失败」原因特征没有标准化或者特征里有 NaN/Inf。解决先检查np.isnan(features).sum()和np.isinf(features).sum()用np.nan_to_num处理。然后确保 Pipeline 里有StandardScaler。5.5 现象分帧后样本量太少模型训不起来原因窗长太长或重叠太低。1 秒窗、50% 重叠60 秒数据只有 119 帧。解决缩短窗长到 0.5 秒或者提高重叠到 75%。如果还是不够可以用数据增强——加高斯噪声、时间平移、幅值缩放。6. 进阶技巧用协方差矩阵和 Riemannian 几何提升分类器如果你已经把上面的基线跑通了准确率卡在 70% 左右上不去可以试试 Riemannian 几何方法。这是近几年脑电分类里比较有效的进阶方案核心思路是把每个帧的协方差矩阵映射到切空间在切空间里做线性分类。from pyriemann.estimation import Covariances from pyriemann.tangentspace import TangentSpace from pyriemann.classification import MDM # 计算协方差矩阵 cov Covariances(estimatorlwf).fit_transform(segments) print(f协方差矩阵形状: {cov.shape}) # (帧数, 通道数, 通道数) # 方法一最小距离均值MDM mdm MDM() mdm.fit(cov_train, y_train) score_mdm mdm.score(cov_test, y_test) # 方法二切空间映射 逻辑回归 ts TangentSpace() X_ts ts.fit_transform(cov_train) from sklearn.linear_model import LogisticRegression lr LogisticRegression(max_iter1000) lr.fit(X_ts, y_train) score_lr lr.score(ts.transform(cov_test), y_test)逻辑说明Covariances对每个帧计算通道间的协方差矩阵estimatorlwf是 Ledoit-Wolf 收缩估计比原始协方差更稳定。TangentSpace把协方差矩阵从黎曼流形映射到欧氏空间之后就可以用普通分类器。MDM直接在黎曼流形上算距离不需要映射。参数说明Covariances的estimator可选scm样本协方差、lwfLedoit-Wolf、oasOracle 近似收缩。通道数大于帧长时用lwf或oas。TangentSpace的metric默认是riemann也可以选logeuclid。我自己的习惯是先用 PSD SVM 跑一个基线如果准确率低于 65%说明特征工程有问题回去检查滤波和分帧如果在 65–75% 之间可以试 Riemannian 方法如果已经超过 80%优先检查有没有数据泄漏。这套流程在多个生理信号数据集上验证过比盲目上深度学习靠谱得多。希望帮到你。本文还有配套的精品资源点击获取
脑电波分析实战：从数据读取到Riemannian分类的完整流水线
简介本资源为2020年研究生数学建模竞赛C题「面向康复工程的脑电信号分析和判别模型」的完整代码与数据包面向参加电赛及数学建模竞赛的研究生、本科生以及从事脑电信号处理与生物医学工程方向的学习者。包内共263个文件以Python脚本与编译文件、PNG图表、XLSX数据表、XML配置、TXT说明及DOCX赛题文档为主压缩包约118.59MB涵盖P300脑机接口数据、睡眠分类等实验模块便于复现预处理、特征提取与建模分析全流程。已有45人学习下载。读者可据此获得赛题方案、算法实现与数据组织思路对照练习信号去噪、特征工程及SVM、神经网络等分类预测方法并参考报告撰写框架适合作为竞赛复盘与课程实践的参考素材。1. 脑电波分析赛题从一份 zip 到可复现的建模流水线2020 年研究生数学建模竞赛 C 题给的是脑电波数据赛题包里通常包含原始信号、受试者标签和一份说明文档。很多人拿到这个 zip 的第一反应是「先跑起来再说」结果卡在数据读取上——脑电信号不是普通表格采样率、通道数、事件标记三者对不上后面所有特征工程都是空中楼阁。这道题真正考的不是某个分类器调得多好而是你能不能把一段非平稳、低信噪比的时序信号拆成可解释、可复现的特征再落到一个能自圆其说的模型上。适合正在做数学建模、信号处理课程设计或者想入门生理信号分析的人。下面按「读数据 → 提特征 → 建模型 → 避坑」的顺序把这条流水线讲透代码可以直接抄。2. 脑电数据到底长什么样先搞清楚三个维度再动手2.1 采样率、通道数、事件标记的三角关系脑电数据的标准结构是一个二维矩阵加一条时间轴。矩阵的行是通道电极位置列是采样点。假设采样率是 250 Hz记录 60 秒那就是 250 × 60 15000 个采样点。通道数常见的是 14、32、64 导联赛题数据一般不会超过 64 通道。事件标记trigger是另一条独立的时间序列记录某个时刻发生了什么——比如第 3.2 秒出现刺激第 5.8 秒受试者按键。这三者的关系必须对齐事件标记的时间戳要换算成采样点索引才能从原始矩阵里切出对应片段。常见错误是拿秒数直接当索引去切片结果切出来的数据完全错位。换算公式很简单# 事件时间(秒) 转 采样点索引 sample_index int(event_time_sec * sfreq) # 例如 sfreq250, event_time_sec3.2 - sample_index800参数说明sfreq是采样率单位 Hz必须从数据说明文档里确认不能猜。event_time_sec是事件发生的绝对时间。int()取整是因为索引必须是整数但要注意取整方向——用round()还是int()取决于你对齐精度要求一般int()向下取整就够。2.2 用 Python 把 zip 里的数据读进来赛题包里的数据格式通常是.matMATLAB或.csv。.mat文件用scipy.io.loadmat读.csv用pandas读。先看数据规模再决定内存策略import scipy.io as sio import numpy as np # 读取 mat 文件结构通常是 dict data sio.loadmat(eeg_data.mat) print(data.keys()) # 先看有哪些变量名 # 假设信号变量叫 EEG形状是 (通道数, 采样点) eeg data[EEG] print(f通道数: {eeg.shape[0]}, 采样点数: {eeg.shape[1]}) # 假设采样率存在 fs 变量里 sfreq int(data[fs][0][0]) print(f采样率: {sfreq} Hz, 时长: {eeg.shape[1]/sfreq:.1f} 秒)逻辑说明loadmat返回的是一个字典键名取决于原始 MATLAB 脚本怎么存的。先print(data.keys())确认变量名不要硬编码。eeg.shape返回(通道数, 采样点数)如果读出来是反的说明数据存储时转置过需要手动.T。采样率经常存在一个 1×1 的数组里用[0][0]取标量。参数说明eeg的数据类型通常是float64如果内存吃紧可以转float32精度损失对后续滤波影响很小。sfreq必须和事件标记的时间基准一致如果事件标记用的是毫秒记得先除以 1000。2.3 数据质量检查三行代码排除 80% 的坑读进来之后别急着提特征先做三件事看有没有 NaN、看幅值范围、看通道间相关性。# 检查缺失值 print(fNaN 数量: {np.isnan(eeg).sum()}) # 检查幅值范围脑电典型幅值在 ±100 μV 以内 print(f幅值范围: [{eeg.min():.2f}, {eeg.max():.2f}] μV) # 检查通道间相关性正常脑电通道间有相关性但不会完全一致 corr np.corrcoef(eeg[:5, :]) # 取前5个通道 print(f通道间平均相关: {corr[np.triu_indices(5,1)].mean():.3f})如果 NaN 数量不为零要么是采集时丢包要么是某个通道坏了。幅值超过 ±200 μV 基本可以判定是伪迹眨眼、肌肉活动。通道间相关性如果接近 1.0说明两个电极可能贴在了同一个位置或者短路了。这三步做完你对数据的底子就有数了。3. 特征工程把一段波形变成模型能吃的数字3.1 滤波带通 陷波顺序不能反脑电信号的有效频段通常在 0.5–45 Hz。工频干扰50 Hz 或 60 Hz必须去掉否则功率谱上会有一根刺。滤波顺序有讲究先带通再陷波还是先陷波再带通我一般先带通把高频噪声整体压下去再用陷波精准打掉工频。反过来做的话陷波器在宽频噪声上的瞬态响应会更难处理。from scipy.signal import butter, filtfilt, iirnotch def bandpass_filter(data, lowcut, highcut, fs, order4): nyq 0.5 * fs b, a butter(order, [lowcut/nyq, highcut/nyq], btypeband) return filtfilt(b, a, data, axis1) def notch_filter(data, freq, fs, quality30): b, a iirnotch(freq, quality, fs) return filtfilt(b, a, data, axis1) # 先带通 0.5-45 Hz再陷波 50 Hz eeg_filtered bandpass_filter(eeg, 0.5, 45, sfreq) eeg_filtered notch_filter(eeg_filtered, 50, sfreq)逻辑说明filtfilt做的是零相位滤波前后各滤波一次避免相位偏移。butter的order4是常用值阶数太高会不稳定太低则过渡带太宽。iirnotch的quality30控制陷波带宽值越大陷波越窄一般 30–50 之间。参数说明lowcut和highcut根据你的分析目标调整。如果关注慢波delta 波低频可以放到 0.5 Hz如果只看高频振荡gamma 波高频可以放到 100 Hz 但采样率要够。axis1表示沿时间轴滤波不要搞错方向。3.2 分帧与加窗把连续信号切成模型能吃的片段脑电信号是非平稳的整段做 FFT 没有意义。标准做法是分帧每帧 1–2 秒帧间重叠 50%。加窗通常是汉宁窗是为了减少频谱泄漏。def segment_signal(data, sfreq, window_sec1.0, overlap0.5): window_len int(window_sec * sfreq) step int(window_len * (1 - overlap)) segments [] for start in range(0, data.shape[1] - window_len, step): seg data[:, start:start window_len] # 加汉宁窗 window np.hanning(window_len) seg seg * window segments.append(seg) return np.array(segments) # 形状: (帧数, 通道数, 窗长) segments segment_signal(eeg_filtered, sfreq, window_sec1.0, overlap0.5) print(f分帧后形状: {segments.shape})逻辑说明window_len是每帧的采样点数step是帧间步长。重叠 50% 意味着步长是窗长的一半。加窗时seg * window是逐点相乘窗函数两端趋近于零中间为一能有效抑制截断处的频谱泄漏。参数说明window_sec选 1 秒是脑电分析的常用值太短频率分辨率不够太长则时间分辨率下降。overlap取 0.5 是经验值增加到 0.75 可以增加样本量但计算量翻倍。3.3 功率谱密度与频带能量最稳的基线特征功率谱密度PSD是脑电分析里最基础也最稳的特征。用 Welch 方法估计 PSD然后按频带积分得到 delta0.5–4 Hz、theta4–8 Hz、alpha8–13 Hz、beta13–30 Hz、gamma30–45 Hz五个频带的能量。from scipy.signal import welch def band_power(segments, sfreq, bandsNone): if bands is None: bands {delta: (0.5, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30), gamma: (30, 45)} features [] for seg in segments: freqs, psd welch(seg, sfreq, npersegseg.shape[1], axis1) band_feats [] for name, (low, high) in bands.items(): idx np.logical_and(freqs low, freqs high) band_feats.append(np.trapz(psd[:, idx], freqs[idx], axis1)) features.append(np.concatenate(band_feats)) return np.array(features) # 形状: (帧数, 通道数*频带数) features band_power(segments, sfreq) print(f特征矩阵形状: {features.shape})逻辑说明welch返回频率轴和对应的功率谱密度。np.trapz做梯形积分把频带内的 PSD 累加成能量值。每个通道每个频带一个值所以特征维度是通道数 × 5。参数说明nperseg设为窗长表示不进一步分段直接用整帧做 FFT。如果帧长超过 2 秒可以设npersegsfreq*2来增加频率分辨率。np.trapz的axis1表示沿频率轴积分。3.4 时域特征补充Hjorth 参数与过零率频域特征之外时域特征计算快、解释性强。Hjorth 参数包括活动度方差、移动度一阶导数标准差与信号标准差之比、复杂度移动度的一阶导数与移动度之比。过零率反映信号在零线附近振荡的频率。def hjorth_params(segments): activity np.var(segments, axis2) diff1 np.diff(segments, axis2) mobility np.sqrt(np.var(diff1, axis2) / (activity 1e-10)) diff2 np.diff(diff1, axis2) mobility_d1 np.sqrt(np.var(diff2, axis2) / (np.var(diff1, axis2) 1e-10)) complexity mobility_d1 / (mobility 1e-10) return activity, mobility, complexity act, mob, comp hjorth_params(segments) print(f活动度形状: {act.shape}) # (帧数, 通道数)逻辑说明np.var沿时间轴算方差。np.diff算一阶差分近似导数。分母加1e-10防止除零。三个参数分别反映信号的幅度、频率和频率变化率。参数说明Hjorth 参数对窗长不敏感1 秒和 2 秒窗算出来的值差异不大。但如果信号幅值很小接近零1e-10的保护可能不够可以改成1e-6。4. 建模与验证从特征矩阵到分类结果4.1 分类器选型为什么我优先用 SVM 而不是深度学习赛题数据量通常不大——几十个受试者每人几十到几百帧。这种规模下SVM 配合 RBF 核往往比深度学习更稳。深度学习需要大量数据才能学到有意义的表示小样本上容易过拟合。随机森林也可以但特征维度高时解释性不如 SVM 的核技巧直观。from sklearn.svm import SVC from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.model_selection import cross_val_score, StratifiedKFold # 构建流水线标准化 SVM pipe Pipeline([ (scaler, StandardScaler()), (svm, SVC(kernelrbf, C1.0, gammascale)) ]) # 5 折交叉验证 cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(pipe, features, labels, cvcv, scoringaccuracy) print(f交叉验证准确率: {scores.mean():.3f} ± {scores.std():.3f})逻辑说明StandardScaler把每个特征维度标准化到零均值单位方差SVM 对尺度敏感这一步不能省。gammascale自动根据特征数和方差设置核宽度比手动调省事。StratifiedKFold保证每折的类别比例一致。参数说明C是正则化参数越大越容易过拟合越小越容易欠拟合。gamma控制核函数的宽度scale是1/(n_features * X.var())的简写。如果交叉验证准确率波动大先检查特征里有没有异常值。4.2 特征重要性分析用排列重要性看哪些特征在起作用SVM 没有直接的特征重要性输出但可以用排列重要性permutation importance来评估。思路是打乱某个特征的值看模型性能下降多少下降越多说明该特征越重要。from sklearn.inspection import permutation_importance # 先在训练集上拟合 pipe.fit(X_train, y_train) # 计算排列重要性 result permutation_importance(pipe, X_test, y_test, n_repeats10, random_state42) importance result.importances_mean # 假设特征名按通道和频带排列 feature_names [fch{i}_{band} for i in range(n_channels) for band in [delta,theta,alpha,beta,gamma]] top_idx np.argsort(importance)[-10:] for i in top_idx: print(f{feature_names[i]}: {importance[i]:.4f})逻辑说明permutation_importance对每个特征重复打乱n_repeats次计算性能下降的均值和标准差。importances_mean是平均下降幅度越大越重要。参数说明n_repeats10是精度和计算量的折中特征多时可以降到 5。random_state固定随机种子保证可复现。4.3 交叉验证的坑别用普通 KFold 做分类分类问题必须用分层交叉验证StratifiedKFold否则某折可能全是某一类准确率直接崩。另外标准化必须在每折内部做不能先标准化再交叉验证——那叫数据泄漏。# 错误做法先标准化再交叉验证 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 泄漏了测试集信息 scores_wrong cross_val_score(SVC(), X_scaled, y, cv5) # 正确做法用 Pipeline 把标准化包进去 pipe Pipeline([(scaler, StandardScaler()), (svm, SVC())]) scores_right cross_val_score(pipe, X, y, cv5)逻辑说明fit_transform在全体数据上计算均值和方差测试集的信息泄漏到了训练过程。Pipeline保证每折的标准化只在该折的训练集上拟合。参数说明如果特征维度远大于样本数标准化后 SVM 可能仍然过拟合此时可以加 PCA 降维但 PCA 也要放进 Pipeline。5. 避坑与排查脑电分析里最容易翻车的五个地方5.1 现象交叉验证准确率 99%换一批数据掉到 50%原因特征里混入了标签泄漏。比如分帧时把同一受试者的帧同时分到了训练集和测试集模型记住了受试者的个体特征而不是任务相关特征。解决按受试者划分训练集和测试集而不是按帧随机划分。用GroupKFold指定受试者 ID 作为分组变量。5.2 现象滤波后信号幅值变得极小或极大原因filtfilt的滤波器阶数太高导致数值不稳定或者lowcut设得太低比如 0.1 Hz而数据长度不够。解决降低滤波器阶数到 2–4 阶或者改用sosfiltfilt二阶节形式提高数值稳定性。lowcut不要低于 0.5 Hz。5.3 现象PSD 图上 50 Hz 处仍然有尖峰原因陷波滤波器的quality参数设得太高陷波带宽太窄没有完全覆盖工频干扰的频谱扩散。解决把quality降到 20–30或者改用自适应滤波。如果工频干扰特别严重可以在带通之前先做一次陷波。5.4 现象SVM 训练报错「收敛失败」原因特征没有标准化或者特征里有 NaN/Inf。解决先检查np.isnan(features).sum()和np.isinf(features).sum()用np.nan_to_num处理。然后确保 Pipeline 里有StandardScaler。5.5 现象分帧后样本量太少模型训不起来原因窗长太长或重叠太低。1 秒窗、50% 重叠60 秒数据只有 119 帧。解决缩短窗长到 0.5 秒或者提高重叠到 75%。如果还是不够可以用数据增强——加高斯噪声、时间平移、幅值缩放。6. 进阶技巧用协方差矩阵和 Riemannian 几何提升分类器如果你已经把上面的基线跑通了准确率卡在 70% 左右上不去可以试试 Riemannian 几何方法。这是近几年脑电分类里比较有效的进阶方案核心思路是把每个帧的协方差矩阵映射到切空间在切空间里做线性分类。from pyriemann.estimation import Covariances from pyriemann.tangentspace import TangentSpace from pyriemann.classification import MDM # 计算协方差矩阵 cov Covariances(estimatorlwf).fit_transform(segments) print(f协方差矩阵形状: {cov.shape}) # (帧数, 通道数, 通道数) # 方法一最小距离均值MDM mdm MDM() mdm.fit(cov_train, y_train) score_mdm mdm.score(cov_test, y_test) # 方法二切空间映射 逻辑回归 ts TangentSpace() X_ts ts.fit_transform(cov_train) from sklearn.linear_model import LogisticRegression lr LogisticRegression(max_iter1000) lr.fit(X_ts, y_train) score_lr lr.score(ts.transform(cov_test), y_test)逻辑说明Covariances对每个帧计算通道间的协方差矩阵estimatorlwf是 Ledoit-Wolf 收缩估计比原始协方差更稳定。TangentSpace把协方差矩阵从黎曼流形映射到欧氏空间之后就可以用普通分类器。MDM直接在黎曼流形上算距离不需要映射。参数说明Covariances的estimator可选scm样本协方差、lwfLedoit-Wolf、oasOracle 近似收缩。通道数大于帧长时用lwf或oas。TangentSpace的metric默认是riemann也可以选logeuclid。我自己的习惯是先用 PSD SVM 跑一个基线如果准确率低于 65%说明特征工程有问题回去检查滤波和分帧如果在 65–75% 之间可以试 Riemannian 方法如果已经超过 80%优先检查有没有数据泄漏。这套流程在多个生理信号数据集上验证过比盲目上深度学习靠谱得多。希望帮到你。本文还有配套的精品资源点击获取
