脑电波分析实战：从原始EEG信号到特征提取与分类的完整链路

简介这份资源是2020年研究生数学建模竞赛C题的完整备赛包聚焦面向康复工程的脑电信号分析与判别建模适合参加电赛、数模竞赛的研究生及对生物医学信号处理感兴趣的读者。压缩包共263个文件约118.59MB包含22个Python脚本、83个pyc编译文件、25个xlsx数据表、80张png结果图、20个xml配置及若干txt说明与docx报告覆盖数据预处理、特征提取、模型训练到结果可视化的完整流程。资源围绕P300脑机接口数据与睡眠分类任务提供可运行的代码与配套数据便于读者复现线性回归、聚类、主成分分析及支持向量机、神经网络等建模方法理解脑电信号去噪、频域变换与分类预测的实现细节。目前已有45人学习下载适合希望系统掌握脑电分析建模思路、对照赛题查漏补缺的参赛者参考。1. 脑电波分析赛题从原始信号到可分特征的那条链路2020年研究生数学建模竞赛C题给了一批脑电波数据要求做特征提取与分类判别。很多人拿到数据第一反应是直接上深度学习结果被采样率不统一、通道数不一致、标签对齐错位这三座大山挡住连baseline都跑不出来。这道题真正的价值不在模型多深而在于把「原始EEG信号→预处理→特征工程→分类器」这条链路走通。适合正在做数学建模、信号处理课程设计或者手上有脑电数据但不知道怎么下手的同学。代码和数据拿到手只是起点关键是理解每一步为什么这么做、参数怎么定、哪里容易翻车。下面按我实际复现这套方案的顺序展开中间会给出可直接运行的Python代码和参数说明。2. 数据到手先别急着建模EEG文件结构与预处理链路2.1 先搞清楚数据长什么样脑电波分析的第一步不是写模型而是把数据读进来看看它到底长什么样。常见的EEG数据格式有.mat、.csv、.edf、.set几种。数学建模竞赛给的数据通常是.mat或.csv因为方便参赛者用MATLAB或Python直接加载。拿到数据后先确认三件事采样率是多少、有多少个通道、每个被试的记录时长和标签怎么对应。我一般会先跑一段探查代码把数据的形状、采样率、通道名称、标签分布全部打印出来。这一步看起来简单但很多翻车案例都是因为跳过了这步——比如把采样率480Hz的数据当成250Hz处理后面滤波器的截止频率全错特征完全不可用。import scipy.io as sio import numpy as np # 加载.mat文件squeeze_meTrue可以去掉多余维度 data sio.loadmat(eeg_data.mat, squeeze_meTrue) print(文件中的变量名, [k for k in data.keys() if not k.startswith(__)]) # 假设数据存在EEG变量中结构为[通道数 x 采样点] eeg data[EEG] print(数据形状, eeg.shape) print(数据类型, eeg.dtype) # 如果采样率单独存储 if fs in data: fs data[fs] print(f采样率{fs} Hz) print(f记录时长{eeg.shape[1] / fs:.2f} 秒)这段代码的逻辑很直接先看文件里有哪些变量再确认EEG矩阵的维度含义。squeeze_meTrue是为了避免MATLAB保存时多出来的冗余维度。参数方面eeg.shape如果是(通道数, 采样点)那shape[0]就是通道数如果反过来后面所有处理都要转置。采样率fs决定了后续滤波器的归一化频率这个参数绝对不能猜。2.2 预处理滤波、去伪迹、分段原始EEG信号里混着工频干扰、眼电伪迹、肌电噪声。不做预处理直接提特征等于在噪声里找信号。标准流程是带通滤波→陷波滤波→去伪迹→分段。带通滤波保留0.5-45Hz或根据任务调整目的是去掉基线漂移和大部分高频噪声。陷波滤波去掉50Hz工频国内电网频率。去伪迹常用方法有ICA独立成分分析或者简单的幅值阈值法。分段则根据实验范式把连续信号切成试次。from scipy.signal import butter, filtfilt, iirnotch def bandpass_filter(signal, lowcut, highcut, fs, order4): 带通滤波器signal形状为[通道数 x 采样点] nyq fs / 2.0 b, a butter(order, [lowcut/nyq, highcut/nyq], btypeband) # filtfilt零相位滤波避免时间偏移 return filtfilt(b, a, signal, axis1) def notch_filter(signal, freq, fs, quality30): 陷波滤波器去除工频干扰 b, a iirnotch(freq, quality, fs) return filtfilt(b, a, signal, axis1) # 实际调用 fs 500 # 假设采样率500Hz eeg_filtered bandpass_filter(eeg, 0.5, 45, fs) eeg_filtered notch_filter(eeg_filtered, 50, fs) print(滤波后数据范围, eeg_filtered.min(), ~, eeg_filtered.max())butter函数的order4是常用值阶数太高会导致数值不稳定太低则滤波效果差。filtfilt做双向滤波相位不偏移这对后续分段和特征提取很关键。陷波滤波器的quality30控制带宽值越大陷波越窄。滤波后检查数据范围如果出现异常大的值说明滤波器不稳定或者数据本身有极端伪迹。去伪迹这一步如果数据质量还行可以用简单阈值法超过±100μV的片段标记为伪迹并剔除。如果伪迹严重就得上ICA。ICA的代码稍长核心思路是把多通道信号分解成独立成分手动或自动识别眼电成分并置零再重构回通道空间。分段则根据标签文件把连续数据切成固定长度的epoch。比如每个试次从刺激呈现后0.5秒到2.5秒那就是2秒的数据段。分段后得到三维矩阵[试次数 x 通道数 x 采样点]这是后续特征提取的标准输入格式。注意分段时一定要核对标签对齐。常见错误是标签文件里的顺序和数据段顺序不一致导致训练出来的模型准确率还不如随机猜。3. 特征工程从时域、频域到空域的三条路3.1 时域特征均值、方差、Hjorth参数时域特征是最直观的。对每个通道的每个epoch计算均值、方差、偏度、峰度、过零率。这些特征计算快物理意义明确。Hjorth参数是EEG分析里常用的三个时域指标活动度Activity即方差、移动度Mobility一阶导数标准差与原始标准差之比、复杂度Complexity移动度的一阶导数与移动度之比。from scipy.stats import skew, kurtosis def time_domain_features(epoch): epoch形状为[通道数 x 采样点]返回每个通道的时域特征 features [] for ch in range(epoch.shape[0]): sig epoch[ch, :] feat { mean: np.mean(sig), var: np.var(sig), skew: skew(sig), kurt: kurtosis(sig), zcr: np.sum(np.diff(np.sign(sig)) ! 0) / len(sig) } # Hjorth参数 d1 np.diff(sig) d2 np.diff(d1) activity np.var(sig) mobility np.sqrt(np.var(d1) / activity) if activity 0 else 0 complexity np.sqrt(np.var(d2) / np.var(d1)) / mobility if np.var(d1) 0 and mobility 0 else 0 feat[hjorth_activity] activity feat[hjorth_mobility] mobility feat[hjorth_complexity] complexity features.append(feat) return features这段代码对每个通道独立计算特征。zcr过零率反映信号振荡频率的粗略估计。Hjorth参数里activity就是方差mobility衡量信号斜率的变化率complexity衡量信号与纯正弦波的偏离程度。参数方面没有需要调的但要注意如果activity为0信号完全平坦要做除零保护。时域特征的优点是计算快、可解释性强缺点是对噪声敏感且无法捕捉频率信息。在脑电分析里时域特征通常作为辅助主力还是频域特征。3.2 频域特征功率谱密度与各频带能量比脑电的节律性是其核心特征。δ波0.5-4Hz、θ波4-8Hz、α波8-13Hz、β波13-30Hz、γ波30-45Hz各自对应不同的认知状态。频域特征就是计算每个通道在各频带的能量或功率谱密度。from scipy.signal import welch def frequency_domain_features(epoch, fs): 计算各频带能量占比 bands { delta: (0.5, 4), theta: (4, 8), alpha: (8, 13), beta: (13, 30), gamma: (30, 45) } features [] for ch in range(epoch.shape[0]): sig epoch[ch, :] # 计算功率谱密度 freqs, psd welch(sig, fs, npersegmin(256, len(sig))) total_power np.trapz(psd, freqs) feat {} for band_name, (low, high) in bands.items(): idx np.where((freqs low) (freqs high))[0] if len(idx) 0: band_power np.trapz(psd[idx], freqs[idx]) feat[f{band_name}_power] band_power feat[f{band_name}_ratio] band_power / total_power if total_power 0 else 0 else: feat[f{band_name}_power] 0 feat[f{band_name}_ratio] 0 features.append(feat) return featureswelch方法用分段平均来估计功率谱比直接FFT更平滑。nperseg是每段长度太小则频率分辨率低太大则平滑效果差。一般取256或512。np.trapz做梯形积分把功率谱密度在频带内积分得到频带能量。ratio是频带能量占总能量的比例这个特征比绝对能量更鲁棒因为不同被试的绝对能量差异很大。频域特征的关键参数是频带划分。标准划分如上但如果任务涉及特定频率比如SSVEP需要根据刺激频率调整。另外如果数据分段较短比如1秒频率分辨率只有1Hz低频段的能量估计会不准。这时候要么加长分段要么用多窗谱估计如DPSS。3.3 空域特征通道间相关性与共空间模式空域特征利用多通道之间的空间关系。最简单的做法是计算通道两两之间的相关系数矩阵取上三角作为特征。更高级的是共空间模式CSP它通过最大化两类信号的方差比来构造空间滤波器是运动想象脑电分类的经典方法。from sklearn.covariance import LedoitWolf def channel_correlation_features(epoch): 计算通道间相关系数矩阵的上三角 n_channels epoch.shape[0] corr_matrix np.corrcoef(epoch) # 取上三角不含对角线 triu_indices np.triu_indices(n_channels, k1) return corr_matrix[triu_indices] # CSP的简化实现二分类 def csp_fit(X_class1, X_class2): X形状为[试次数 x 通道数 x 采样点] # 计算两类平均协方差矩阵 cov1 np.mean([np.cov(x) for x in X_class1], axis0) cov2 np.mean([np.cov(x) for x in X_class2], axis0) # 广义特征值分解 eigvals, eigvecs np.linalg.eig(np.linalg.inv(cov1 np.eye(cov1.shape[0])*1e-6) cov2) # 按特征值排序 idx np.argsort(eigvals)[::-1] eigvecs eigvecs[:, idx] return eigvecs def csp_transform(epoch, filters): 用CSP滤波器投影 return filters.T epoch通道相关性特征计算简单但维度是n_channels*(n_channels-1)/2通道多的时候维度爆炸。CSP的核心是广义特征值分解np.linalg.inv(cov1 I*1e-6)里的正则项防止协方差矩阵奇异。特征值越大说明该空间滤波器对两类信号的区分度越高。通常取前m个和后m个特征向量m2或3作为滤波器。CSP的坑在于它是有监督的必须用训练数据拟合滤波器再应用到测试数据。如果直接用全部数据拟合就是数据泄露交叉验证结果会虚高。另外CSP对伪迹很敏感预处理没做好CSP滤波器会学到伪迹而不是脑电特征。4. 分类器选型与交叉验证别在验证集上翻车4.1 从SVM到随机森林小样本下的选择脑电数据通常是「小样本高维度」几十个被试每个被试几十个试次特征维度可能上百。这种数据下SVM尤其是线性SVM和正则化逻辑回归是首选因为它们的泛化能力在样本少时比深度学习更稳。随机森林可以作为baseline但要注意它容易过拟合。from sklearn.svm import SVC from sklearn.ensemble import RandomForestClassifier from sklearn.linear_model import LogisticRegression from sklearn.preprocessing import StandardScaler from sklearn.pipeline import Pipeline from sklearn.model_selection import StratifiedKFold, cross_val_score # 构建pipeline标准化 分类器 def build_classifier(clf_typesvm): if clf_type svm: clf SVC(kernelrbf, C1.0, gammascale, probabilityTrue) elif clf_type rf: clf RandomForestClassifier(n_estimators100, max_depth5, random_state42) else: clf LogisticRegression(C1.0, max_iter1000) return Pipeline([ (scaler, StandardScaler()), (clf, clf) ]) # 交叉验证 def evaluate_classifier(X, y, clf_typesvm): X形状为[样本数 x 特征数]y为标签 clf build_classifier(clf_type) cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) scores cross_val_score(clf, X, y, cvcv, scoringaccuracy) return scores.mean(), scores.std() # 示例 # X np.array(features) # 特征矩阵 # y np.array(labels) # 标签 # mean_acc, std_acc evaluate_classifier(X, y, svm) # print(fSVM准确率{mean_acc:.3f} ± {std_acc:.3f})StandardScaler必须放在pipeline里不能提前对全部数据做标准化否则测试集的均值和方差信息泄露到训练中。SVM的C参数控制正则化强度C越大越容易过拟合。gammascale是sklearn的默认值等于1/(n_features * X.var())。随机森林的max_depth5是保守设置防止树太深记住噪声。交叉验证用StratifiedKFold保证每折的类别比例一致。shuffleTrue打乱顺序避免数据按类别排序导致的偏差。如果数据有被试间差异应该用「留一被试交叉验证」LOSO即每次留一个被试做测试其余被试做训练。这样评估的是模型跨被试的泛化能力比随机划分更严格。4.2 交叉验证的三种划分策略与适用场景随机K折、留一被试、时间序列划分这三种策略对应不同的评估目标。随机K折适合被试内分类同一个人的数据分训练测试留一被试适合被试间分类跨人泛化时间序列划分适合在线系统用过去预测未来。from sklearn.model_selection import LeaveOneGroupOut def loso_cv(X, y, groups): 留一被试交叉验证groups为被试ID logo LeaveOneGroupOut() clf build_classifier(svm) scores [] for train_idx, test_idx in logo.split(X, y, groups): clf.fit(X[train_idx], y[train_idx]) score clf.score(X[test_idx], y[test_idx]) scores.append(score) return np.mean(scores), np.std(scores)LeaveOneGroupOut的groups参数指定每个样本属于哪个被试。如果10个被试就训练10次每次留一个被试。这种评估方式得到的准确率通常比随机K折低10-20个百分点但更接近实际应用场景。如果竞赛要求报告被试间分类结果必须用这种划分。注意如果数据里每个被试的试次数差异很大LOSO的每折测试集大小不一计算平均准确率时最好按试次数加权而不是简单平均。5. 避坑与排查脑电分析里最容易翻车的五个地方5.1 采样率不统一导致滤波器失效现象滤波后信号出现异常振荡或者频带能量全为0。原因不同被试的采样率不同但代码里用了固定的fs。解决在读数据时逐个文件检查采样率统一重采样到同一频率如250Hz或500Hz。重采样用scipy.signal.resample或mne的resample函数。5.2 标签对齐错位导致准确率异常现象模型准确率在50%左右徘徊或者交叉验证方差极大。原因标签文件里的顺序和数据段顺序不一致或者标签编码时0/1反了。解决打印前10个样本的标签和对应的数据段人工核对几个。另外用np.unique(y, return_countsTrue)检查类别分布如果严重不平衡准确率就不是好指标要看AUC或F1。5.3 数据泄露导致交叉验证虚高现象交叉验证准确率95%但换一批数据就掉到60%。原因标准化、特征选择、CSP拟合在交叉验证之前对全部数据做了。解决把所有预处理和特征提取步骤放进pipeline或者严格在每折训练集上拟合再应用到测试集。特征选择尤其容易泄露SelectKBest必须放在pipeline里。5.4 伪迹未剔除导致CSP学到噪声现象CSP滤波器的空间模式看起来像眼电或肌电的拓扑图。原因伪迹幅度远大于脑电CSP会优先最大化伪迹的方差比。解决预处理阶段用ICA或阈值法剔除伪迹。如果伪迹比例超过20%考虑降低阈值或手动检查。5.5 频带划分与任务不匹配现象频域特征区分度低分类器学不到东西。原因用了标准频带划分但任务相关的频率成分不在这些频带里。解决先画功率谱图看两类信号在哪些频率上有差异再据此调整频带。比如SSVEP任务要看刺激频率及其谐波运动想象要看mu波8-13Hz和beta波13-30Hz的事件相关去同步。6. 把特征工程做扎实一个可复现的完整流程与调参技巧整套流程跑通后真正决定分类效果的不是分类器多复杂而是特征工程做得多扎实。我一般会按这个顺序迭代先跑时域频域特征用线性SVM看baseline如果准确率低于70%检查预处理和标签对齐如果70-85%加空域特征通道相关性或CSP如果85%以上再考虑调分类器参数。调参方面SVM的C和gamma用网格搜索但搜索范围不要太大。C取[0.1, 1, 10, 100]gamma取[0.001, 0.01, 0.1, 1]。随机森林的n_estimators取100-500max_depth取3-10。注意调参必须在交叉验证内部做不能在全量数据上调完再交叉验证否则又是数据泄露。from sklearn.model_selection import GridSearchCV def tune_svm(X, y, groupsNone): 网格搜索SVM参数如果提供groups则用LOSO param_grid { clf__C: [0.1, 1, 10, 100], clf__gamma: [0.001, 0.01, 0.1, 1] } clf build_classifier(svm) if groups is not None: cv LeaveOneGroupOut() grid GridSearchCV(clf, param_grid, cvcv, scoringaccuracy, n_jobs-1) grid.fit(X, y, groupsgroups) else: cv StratifiedKFold(n_splits5, shuffleTrue, random_state42) grid GridSearchCV(clf, param_grid, cvcv, scoringaccuracy, n_jobs-1) grid.fit(X, y) print(f最佳参数{grid.best_params_}) print(f最佳交叉验证准确率{grid.best_score_:.3f}) return grid.best_estimator_GridSearchCV的n_jobs-1用全部CPU核心并行。如果数据量大可以用RandomizedSearchCV随机采样参数组合速度更快。groups参数在LOSO时传入被试IDgrid.fit(X, y, groupsgroups)会自动按组划分。最后说一个我踩过的坑有次跑完流程发现准确率只有55%排查了半天以为是特征问题最后发现是数据加载时把两个被试的文件顺序搞反了标签全错位。从那以后我养成了一个习惯——每次加载完数据先打印前5个样本的标签和对应的原始信号均值人工确认一眼。这个习惯帮我省了至少十次返工。希望帮到你。本文还有配套的精品资源点击获取