从EEGLAB工具箱到实战：MATLAB脑电（EEG）数据处理全流程解析

1. 环境配置与数据导入第一次打开MATLAB准备处理脑电数据时我对着空荡荡的命令窗口发呆了十分钟。后来才发现原来EEGLAB这个神器需要先安装才能用。现在每次看到新手同事犯同样的错误我都会想起自己当年的窘态。安装EEGLAB其实很简单就像给手机装个新APP。先去官网下载最新版本解压后把整个文件夹扔到MATLAB的toolbox目录下。然后在MATLAB命令行输入addpath(genpath(你的EEGLAB文件夹路径)); eeglab这个genpath命令特别重要它能自动添加所有子目录。我有次偷懒只加了主目录结果运行ICA时各种报错排查了半天才发现是漏了子文件夹。数据导入更是个技术活。不同实验室用的采集设备五花八门从NeuroScan到BrainProducts每种设备生成的文件格式都不同。EEGLAB支持超过20种格式我最常用的是这三种方式% 对于.edf格式 EEG pop_biosig(subj01.edf); % 对于BrainVision的.vhdr/.vmrk/.eeg三件套 EEG pop_loadbv(文件夹路径, subj01.vhdr); % 对于NeuroScan的.cnt文件 EEG pop_loadcnt(subj01.cnt);有个坑要特别注意采样率不一致问题。有次导入的数据显示采样率是1000Hz但实际只有500Hz导致后续滤波全乱了。现在我都会先用这个命令检查fprintf(采样率: %d Hz\n, EEG.srate);2. 预处理全流程详解2.1 电极定位与重参考第一次做电极定位时我把所有电极位置都标在了头顶上活像个刺猬。后来才知道10-20系统有标准坐标文件。现在我的标准操作流程是% 加载标准电极位置 EEG pop_chanedit(EEG, lookup,standard-10-5-cap385.elp); % 可视化检查 pop_chanloc(EEG, plotrad, 0.6);重参考的选择直接影响结果。有篇论文比较了8种参考方式发现乳突参考对听觉实验最友好而全脑平均参考适合视觉实验。我的常用代码% 双侧乳突参考 EEG pop_reref(EEG, [find(strcmp({EEG.chanlocs.labels},TP9))... find(strcmp({EEG.chanlocs.labels},TP10))]); % 全脑平均参考 EEG pop_reref(EEG, []);2.2 滤波与去噪实战滤波参数设置是个经验活。有次我设了0.1Hz高通滤波结果把慢波全滤没了被导师好一顿说。现在我的黄金法则是% 带通滤波 (0.5-40Hz) EEG pop_eegfiltnew(EEG, 0.5, 40, [], 0, [], 1); % 50Hz工频陷波 EEG pop_eegfiltnew(EEG, 49, 51, [], 1, [], 1);ICA去伪迹最考验眼力。我总结了个快速识别伪迹成分的口诀眨眼看前额肌电看颞区心电看后枕。操作步骤% 运行ICA EEG pop_runica(EEG, extended,1,interupt,on); % 可视化成分 pop_selectcomps(EEG, 1:20); % 剔除眨眼成分(假设成分1和3是眨眼) EEG pop_subcomp(EEG, [1 3], 0);3. 分段与特征提取3.1 事件相关电位分析分段时最容易犯的错误是基线校正时间窗不对。有次我设成[-1 0]秒结果刺激前1秒根本没数据。标准做法% 分段 (-0.2到1秒) EEG pop_epoch(EEG, {event_code}, [-0.2 1.0]); % 基线校正 (-0.2到0秒) EEG pop_rmbase(EEG, [-200 0]);时域平均时要注意剔除异常试次。我常用这个三标准差法则% 计算各试次均值 trial_means squeeze(mean(EEG.data,2)); % 找出异常试次 outliers find(abs(zscore(trial_means)) 3); % 剔除异常试次 EEG pop_rejepoch(EEG, outliers, 0);3.2 时频分析技巧小波变换的参数设置直接影响结果。经过多次调试我发现这些参数最适合ERP分析% 时频分析参数 cycles [3 0.5]; % 低频用较少周期高频用较多 freqs 2:1:40; % 频率范围 times -0.5:0.01:1.5; % 时间窗 % 计算ERSP [ersp,itc,powbase,times,freqs] ... pop_newtimef(EEG, 1, Cz, [-200 1500], ... cycles, cycles, freqs, freqs);时频结果的可视化也有讲究。用jet色谱虽然好看但不科学现在推荐用parulafigure; contourf(times, freqs, ersp, 40, linecolor,none); colormap(parula); colorbar; title(时频能量变化);4. 高级分析与机器学习4.1 功能连接分析相位锁定值(PLV)计算要注意volume conduction问题。我的解决方案是% 先计算相位 phase_data angle(hilbert(EEG.data)); % 计算加权PLV plv zeros(EEG.nbchan, EEG.nbchan); for i 1:EEG.nbchan for j i1:EEG.nbchan plv(i,j) abs(mean(exp(1i*(phase_data(i,:,:) - phase_data(j,:,:))),3)); end end % 可视化 imagesc(plv); title(通道间相位锁定值);4.2 机器学习分类用EEG数据做分类时特征工程比算法选择更重要。这是我处理运动想象任务的流程% 提取时频特征 features []; for ch [find(strcmp({EEG.chanlocs.labels},C3))... find(strcmp({EEG.chanlocs.labels},C4))] [ersp,~,~,~] pop_newtimef(EEG, ch, [-1000 2000],... freqs, [8 30], plotersp,off); features [features; squeeze(mean(ersp,3))]; end % 标签准备 labels [ones(1,EEG.trials/2) zeros(1,EEG.trials/2)]; % 交叉验证 cv cvpartition(labels, KFold, 5); acc zeros(cv.NumTestSets,1); for i 1:cv.NumTestSets trainIdx cv.training(i); testIdx cv.test(i); model fitcsvm(features(:,trainIdx), labels(trainIdx)); acc(i) sum(predict(model, features(:,testIdx)) labels(testIdx))/sum(testIdx); end disp([平均准确率: num2str(mean(acc)*100) %]);5. 实战经验与避坑指南数据质量检查这个步骤太容易忽略。有次我分析完整个数据集才发现有个被试的数据全是噪声。现在我的检查清单包括查看原始波形是否有大面积漂移检查阻抗值是否过高如果元数据中有验证事件标记是否与实验log一致批处理脚本的编写也有讲究。我习惯用这种结构subjects {subj01,subj02,subj03}; for s 1:length(subjects) try % 数据导入 EEG pop_loadbv(datafolder, [subjects{s} .vhdr]); % 预处理流程 EEG pop_eegfiltnew(EEG, 0.5, 40); % ...其他处理步骤 % 保存结果 save([subjects{s} _preprocessed.mat], EEG); catch ME fprintf(处理 %s 时出错: %s\n, subjects{s}, ME.message); end end可视化是分析的最后一步但很重要。我常用的组合是单试次波形平均ERP时频能量图拓扑图功能连接矩阵网络图比如绘制ERP图像时这个代码很实用figure; plot(EEG.times, squeeze(mean(EEG.data(strcmp({EEG.chanlocs.labels},Pz),:,:),3))); hold on; plot(EEG.times, squeeze(mean(EEG.data(strcmp({EEG.chanlocs.labels},Cz),:,:),3))); xlabel(时间 (ms)); ylabel(幅值 (μV)); legend(Pz,Cz); title(事件相关电位);