脑电日记 之 MIEEG数据matlab预处理代码实践

已于 2025-04-01 18:59:12 修改
阅读量923 收藏 10
点赞数 9
分类专栏： 脑电学习 文章标签： matlab
于 2025-04-01 17:57:53 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_46234417/article/details/146915464
版权

常导PPT

先上代码

%% 
clear all; 
close all; clc;
eegroot = 'E:\test_data';                                % data path for EEG raw data
cd([eegroot]);                                    

[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;                
%
allSubs = dir([eegroot filesep 'S*']);                % need to check all Subject in the raw data file

isub = 1;          % 1:length(allSubs)，每次对单个被试一个一个的处理                  

allSubs_data = dir([eegroot filesep allSubs(isub).name  filesep 'S*']);


double EEG_tmp;  % Initialize as empty array

for sess = 1:length(allSubs_data)              % 1:length(allSubs_data); % all the data files for each subject
    EEG_tmp(sess) = pop_importNeuracle( {'data.bdf','evt.bdf'}, [allSubs_data(sess).folder filesep allSubs_data(sess).name]);  % load EEG data one by one
    
    str_results = strsplit(allSubs_data(sess).name,'_');         % to rename the event,here need to split the name files
    
    for evt =1:length([EEG_tmp(sess).event])
        EEG_tmp(sess).event(evt).subject = str_results{1};        % subject informaton
        if length(str_results) == 3                             % check the file name
            EEG_tmp(sess).event(evt).task_pre =  [str_results{2} '_' str_results{3}];
        else
            EEG_tmp(sess).event(evt).task_pre = str_results{2};
        end 
        EEG_tmp(sess).event(evt).task = EEG_tmp(sess).event(evt).task_pre(1:end-1);  % task, session, task_type is depend on subject and task
        EEG_tmp(sess).event(evt).session = EEG_tmp(sess).event(evt).task_pre(end);
        
        EEG_tmp(sess).event(evt).task_type = [EEG_tmp(sess).event(evt).task EEG_tmp(sess).event(evt).type];
    end    
end

EEG = pop_mergeset(EEG_tmp,[1:length(allSubs_data)]);     % 这个地方将这个被试各个任务数据进行合并
EEG = eeg_checkset(EEG);
clear EEG_tmp;
% eeglab redraw     

%% Downsample data: 
if EEG.srate>512
    NewSamplingRate = 512;
    EEG = pop_resample(EEG, NewSamplingRate);
end
% channel location
EEG.chanlocs = pop_chanedit(EEG.chanlocs, 'lookup','D:\MyApp\MatLab Working Path\eeglab2024.2\plugins\dipfit\standard_BESA\standard-10-5-cap385.elp'); %check the path of the EEGLAb

% 滤波 对数据的滤波是在0.1-100，同时采用了陷波滤波器49-51将交流电频率从中滤掉
% 带通滤波 0.1-100Hz,当然可以使用之前那种ERPLAB中的函数也是可以的
EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',100);
% 陷波滤波：49-51，去除50Hz的干扰
EEG = pop_eegfiltnew(EEG, 'locutoff',49,'hicutoff',51,'revfilt',1);
% cleanline，然后使用cleanline和dtrend 函数去除实验过程中漂移严重的数据段（这个注意选做项）
EEG = pop_cleanline(EEG, 'bandwidth',2,'chanlist',[1:length(EEG.chanlocs)] ,'computepower',0,'linefreqs',[50 100 150 200 250] ,'normSpectrum',0,'p',0.01,'pad',2,'plotfigures',0,'scanforlines',1,'sigtype','Channels','tau',100,'verb',1,'winsize',4,'winstep',4);
EEG.data = detrend(EEG.data, 0);

% 删除多余的数据通道，数据显示，观察是否有坏电极
EEG = pop_select( EEG, 'nochannel',{'ECG','HEOR','HEOL','VEOU','VEOL'});      % 删掉这些不同的电极
EEG = eeg_checkset( EEG );
eeglab redraw   
pop_eegplot(EEG, 1, 1, 1);                            % 视觉检查，reject掉很糟糕的数据，并且标注坏电极

originalEEG = EEG;   % 对处理坏电极之前对原始数据通道等信息进行保存
%% （选做）-1 标注坏电极，并删除坏电极（根据数据质量情况选择），在ICA之前如果是有坏电极，则要把坏电极删掉
badchannel{1,2}={'Fpz'};
badchannel{1,3}={'Fp2'};

EEG = pop_select( EEG,'nochannel',badchannel{1,2});
EEG = pop_select( EEG,'nochannel',badchannel{1,3});

EEG = eeg_checkset( EEG );
pop_eegplot(EEG, 1, 1, 1);   
%% 截断
EEG = pop_epoch(EEG, {'1','2','3'}, [0 4], 'epochinfo', 'no');         % 
EEG = pop_rmbase( EEG, [],[1:length(EEG.times)]);    % which is to say using the whole epoch to do the baseline correction=pop_rmbase( EEG, [],[1:length(EEG.times)]) = pop_rmbase( EEG, 1000*[EEG.xmin,EEG.xmax]); 
EEG = eeg_checkset(EEG, 'eventconsistency');
% 截断之后visual inspection,检查一遍是否有伪迹和噪声较多的epoch，如果有，将其删除。
pop_eegplot(EEG, 1, 1, 1); 
EEG = pop_saveset(EEG);            % ICA之前数据保存  S***_preICA.set
%% ICA
[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;                
originalEEG = pop_loadset('S220_preICA.set');

EEG=originalEEG;

EEG.data = reshape(EEG.data,size(EEG.data,1),size(EEG.data,2)*size(EEG.data,3));
[COEFF,SCORE,latent,tsquare] = pca(EEG.data');      % princomp change to pca
num_components_to_keep = find(cumsum(latent) ./ sum(latent)>0.9999,1)+1;   %change the value from 0.98-0.99-0.995-0.998-1, if sometimes the ICA is not done
% if use 0.99 to get the ICA result and still there are lots of eye blick signal(out of -100 100), then use 0.98/0.95 to get less component and reject eye component by visual
% sometimes the ICA cannot reject the eye artifact automatic, you need to judje and reject manually
% here run ICA using function runica, and binica is more faster than this one, need to check and install
%
EEG = pop_runica(EEG, 'icatype','runica','pca', num_components_to_keep-2, 'chanind', 1:length(EEG.chanlocs));      %, 'pca', length(EEG.chanlocs));,'pca', 8,
%%
[EEG, com] = SASICA(EEG);  
EEG.rm_ICs = find(EEG.reject.gcompreject);
ncomps = sum(EEG.reject.gcompreject);
[EEG,com] = pop_subcomp(EEG,find(EEG.reject.gcompreject),1);
clear ncomps
%% 平均重参考

EEG = pop_reref( EEG, [], 'refstate',-1);   %here is averaged referene, also can referenced to one certain channel
EEG = eeg_checkset(EEG);

%% 坏电极插值
EEG = pop_interp(EEG,originalEEG.chanlocs,'spherical');   % 之后对坏电极进行插值处理
EEG = eeg_checkset( EEG );
pop_eegplot(EEG, 1, 1, 1);

EEG = pop_saveset(EEG);                       % 最后处理完之后数据保存：S***_postICA.set

拆解分析

· 数据结构介绍

（从整体实验室采集数据中随机选择了5个被试数据的4组运动执行数据进行测试，旨在初学习代码预处理过程）

· 初始化与数据导入

clear all; 
close all; clc;
eegroot = 'E:\test_data';                                % data path for EEG raw data
cd([eegroot]);                                    

[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;                
%
allSubs = dir([eegroot filesep 'S*']);                % need to check all Subject in the raw data file

isub = 1;          % 1:length(allSubs)，每次对单个被试一个一个的处理                  

allSubs_data = dir([eegroot filesep allSubs(isub).name  filesep 'S*']);

  • 功能：初始化环境，设置 EEG 数据路径，并加载 EEG 数据。

  • 解释：

    • eegroot：EEG 数据的根目录。

    • cd([eegroot])：切换到数据目录。

    • dir([eegroot filesep 'S*'])：列出所有以 S 开头的被试文件夹。

    • isub：当前处理的被试索引。

· 加载EEG数据

double EEG_tmp;  % Initialize as empty array

for sess = 1:length(allSubs_data)              % 1:length(allSubs_data); % all the data files for each subject
    EEG_tmp(sess) = pop_importNeuracle( {'data.bdf','evt.bdf'}, [allSubs_data(sess).folder filesep allSubs_data(sess).name]);  % load EEG data one by one
    
    str_results = strsplit(allSubs_data(sess).name,'_');         % to rename the event,here need to split the name files
    
    for evt =1:length([EEG_tmp(sess).event])
        EEG_tmp(sess).event(evt).subject = str_results{1};        % subject informaton
        if length(str_results) == 3                             % check the file name
            EEG_tmp(sess).event(evt).task_pre =  [str_results{2} '_' str_results{3}];
        else
            EEG_tmp(sess).event(evt).task_pre = str_results{2};
        end 
        EEG_tmp(sess).event(evt).task = EEG_tmp(sess).event(evt).task_pre(1:end-1);  % task, session, task_type is depend on subject and task
        EEG_tmp(sess).event(evt).session = EEG_tmp(sess).event(evt).task_pre(end);
        
        EEG_tmp(sess).event(evt).task_type = [EEG_tmp(sess).event(evt).task EEG_tmp(sess).event(evt).type];
    end    
end

  • 功能：加载 EEG 数据，并对事件信息进行重命名和整理。

  • 解释：

    • pop_importNeuracle：加载 EEG 数据。

    • strsplit：将文件名按 _ 分割，提取被试和任务信息。

    • EEG_tmp(sess).event(evt).subject：设置事件的被试信息。

    • EEG_tmp(sess).event(evt).task 和 EEG_tmp(sess).event(evt).session：设置任务和会话信息。

· 合并数据集

EEG = pop_mergeset(EEG_tmp,[1:length(allSubs_data)]);     % 这个地方将这个被试各个任务数据进行合并
EEG = eeg_checkset(EEG);
clear EEG_tmp;

  • 功能：将同一被试的不同任务数据合并。

  • 解释：

    • pop_mergeset：合并多个 EEG 数据集。

    • eeg_checkset：检查 EEG 数据集的完整性。

· 数据预处理

%% Downsample data: 
if EEG.srate>512
    NewSamplingRate = 512;
    EEG = pop_resample(EEG, NewSamplingRate);
end
% channel location
EEG.chanlocs = pop_chanedit(EEG.chanlocs, 'lookup','D:\MyApp\MatLab Working Path\eeglab2024.2\plugins\dipfit\standard_BESA\standard-10-5-cap385.elp'); %check the path of the EEGLAb

% 滤波 对数据的滤波是在0.1-100，同时采用了陷波滤波器49-51将交流电频率从中滤掉
% 带通滤波 0.1-100Hz,当然可以使用之前那种ERPLAB中的函数也是可以的
EEG = pop_eegfiltnew(EEG, 'locutoff',0.1,'hicutoff',100);
% 陷波滤波：49-51，去除50Hz的干扰
EEG = pop_eegfiltnew(EEG, 'locutoff',49,'hicutoff',51,'revfilt',1);
% cleanline，然后使用cleanline和dtrend 函数去除实验过程中漂移严重的数据段（这个注意选做项）
EEG = pop_cleanline(EEG, 'bandwidth',2,'chanlist',[1:length(EEG.chanlocs)] ,'computepower',0,'linefreqs',[50 100 150 200 250] ,'normSpectrum',0,'p',0.01,'pad',2,'plotfigures',0,'scanforlines',1,'sigtype','Channels','tau',100,'verb',1,'winsize',4,'winstep',4);
EEG.data = detrend(EEG.data, 0);

  • 功能：对 EEG 数据进行重采样、滤波、去噪和去趋势处理。

  • 解释：

    • pop_resample：重采样到 512 Hz。

    • pop_chanedit：设置通道位置。

    • pop_eegfiltnew：带通滤波（0.1-100 Hz）和陷波滤波（49-51 Hz）。

    • pop_cleanline：去除交流电干扰。

    • detrend：去除数据中的趋势。

· 删除坏电极

% 删除多余的数据通道，数据显示，观察是否有坏电极
EEG = pop_select( EEG, 'nochannel',{'ECG','HEOR','HEOL','VEOU','VEOL'});      % 删掉这些不同的电极
EEG = eeg_checkset( EEG );
eeglab redraw   
pop_eegplot(EEG, 1, 1, 1);                            % 视觉检查，reject掉很糟糕的数据，并且标注坏电极

  • 功能：删除不需要的通道，并检查数据质量。

  • 解释：

    • pop_select：删除指定的坏电极。

    • pop_eegplot：可视化 EEG 数据。

· 标注和删除坏电极（选做）

%% （选做）-1 标注坏电极，并删除坏电极（根据数据质量情况选择），在ICA之前如果是有坏电极，则要把坏电极删掉
badchannel{1,2}={'Fpz'};
badchannel{1,3}={'Fp2'};

EEG = pop_select( EEG,'nochannel',badchannel{1,2});
EEG = pop_select( EEG,'nochannel',badchannel{1,3});

EEG = eeg_checkset( EEG );
pop_eegplot(EEG, 1, 1, 1);

  • 功能：标注并删除坏电极。

  • 解释：

    • badchannel：标注坏电极。

    • pop_select：删除坏电极。

    • badchannel{1,2} 和 badchannel{1,3} 是一个元胞数组，分别存储了坏电极的名称。

    • 这里标注了两个坏电极：'Fpz' 和 'Fp2'。

· 截段和基线校正

%% 截断
EEG = pop_epoch(EEG, {'1','2','3'}, [0 4], 'epochinfo', 'no');         % 
EEG = pop_rmbase( EEG, [],[1:length(EEG.times)]);    % which is to say using the whole epoch to do the baseline correction=pop_rmbase( EEG, [],[1:length(EEG.times)]) = pop_rmbase( EEG, 1000*[EEG.xmin,EEG.xmax]); 
EEG = eeg_checkset(EEG, 'eventconsistency');
% 截断之后visual inspection,检查一遍是否有伪迹和噪声较多的epoch，如果有，将其删除。
pop_eegplot(EEG, 1, 1, 1); 
EEG = pop_saveset(EEG);            % ICA之前数据保存  S***_preICA.set

  • 功能：截断数据并进行基线校正。

  • 解释：

    • pop_epoch：截断数据。

    • pop_rmbase：基线校正。

    • pop_saveset：保存处理后的数据。

· ICA检测和校正（！）

%% ICA
[ALLEEG EEG CURRENTSET ALLCOM] = eeglab;                
originalEEG = pop_loadset('S220_preICA.set');

EEG=originalEEG;

EEG.data = reshape(EEG.data,size(EEG.data,1),size(EEG.data,2)*size(EEG.data,3));
[COEFF,SCORE,latent,tsquare] = pca(EEG.data');      % princomp change to pca
num_components_to_keep = find(cumsum(latent) ./ sum(latent)>0.9999,1)+1;   %change the value from 0.98-0.99-0.995-0.998-1, if sometimes the ICA is not done
% if use 0.99 to get the ICA result and still there are lots of eye blick signal(out of -100 100), then use 0.98/0.95 to get less component and reject eye component by visual
% sometimes the ICA cannot reject the eye artifact automatic, you need to judje and reject manually
% here run ICA using function runica, and binica is more faster than this one, need to check and install
%
EEG = pop_runica(EEG, 'icatype','runica','pca', num_components_to_keep-2, 'chanind', 1:length(EEG.chanlocs));      %, 'pca', length(EEG.chanlocs));,'pca', 8,
%%
[EEG, com] = SASICA(EEG);  
EEG.rm_ICs = find(EEG.reject.gcompreject);
ncomps = sum(EEG.reject.gcompreject);
[EEG,com] = pop_subcomp(EEG,find(EEG.reject.gcompreject),1);
clear ncomps

  • 功能：运行 ICA，检测并移除伪迹成分。

  • 解释：

    • pca：主成分分析。

    • pop_runica：运行 ICA。

    • SASICA：自动检测伪迹成分。

    • pop_subcomp：移除伪迹成分。

ICA详细学习

运行操作如下

这个窗口是用于选择独立成分（ICA components）的选项设置界面，通常用于 EEG 数据处理中的伪迹检测和校正。以下是每个选项的详细解释：

1. Autocorrelation（自相关性）

  • Enable：启用自相关性检测。

  • Threshold (r)：设置自相关系数的阈值，用于检测低自相关性成分。

  • Lag (ms)：设置自相关性检测的时间滞后（毫秒）。

2. Focal components（焦点成分）

  • Enable：启用焦点成分检测。

  • Threshold (z)：设置焦点成分的阈值，用于检测高度局部化的成分。

3. Focal trial activity（焦点试验活动）

  • Enable：启用焦点试验活动检测。

  • Threshold (z)：设置焦点试验活动的阈值，用于检测试验中高度局部化的活动。

4. Signal to noise ratio（信噪比）

  • Enable：启用信噪比检测。

  • POI (min max ms)：设置感兴趣的信号时间范围（毫秒）。

  • BL (min max ms)：设置基线时间范围（毫秒）。

  • Threshold ratio：设置信噪比的阈值。

5. Dipole fit residual variance（偶极子拟合残差方差）

  • Enable：启用偶极子拟合残差方差检测。

  • Threshold (%)：设置残差方差的百分比阈值。

6. Correlation with EOG（与EOG的关联）

  • Enable：启用与EOG（眼电图）的关联检测。

  • Vertical EOG threshold (r)：设置垂直EOG的关联阈值。

  • Vertical EOG channel：选择垂直EOG通道。

  • Horizontal EOG threshold (r)：设置水平EOG的关联阈值。

  • Horizontal EOG channel：选择水平EOG通道。

7. Correlation with other channel(s)（与其他通道的关联）

  • Enable：启用与其他通道的关联检测。

  • Other channel(s) threshold (r)：设置与其他通道的关联阈值。

  • Other channel(s)：选择其他通道。

8. ADJUST Selection（ADJUST选择）

  • Enable：启用ADJUST伪迹检测方法。

9. FASTER selection（FASTER选择）

  • Enable：启用FASTER伪迹检测方法。

  • Blink channels：选择眨眼通道。

底部按钮

  • Just replot：仅重新绘制图形。

  • No plots：不绘制图形。

  • Help：显示帮助信息。

  • Close：关闭窗口。

  • Compute：执行计算。

Compute...

OK...删除坏成分

然而...我的疑问 --> 这个一个个彩色的脑子图究竟是代表什莫意思呢？？？为什莫这些脑图的数量恰好等于EEG数据的通道数呢？？？难道是每一个图代表一个通道吗？？？

问AI后 得知这每一个彩色脑图代表一个独立成分（IC）在头皮上的电活动分布，红色表示高电活动，蓝色表示低电活动，通过观察这些头皮分布图可以识别出与伪迹相关的成分

那什么又是独立成分（IC）呢还是不解？

再问我的好搭档AI ... ... " 独立成分（Independent Component, IC）是 EEG 数据中的一种信号表示，通过独立成分分析（ICA）算法从 EEG 数据中分离出来的独立信号源。每个独立成分代表一种潜在的信号源，可能是脑电活动，也可能是伪迹（如眼动、肌电等）"

暂且理解独立成分就是一种通过一种强大算法 从包罗众多混合信息的脑电数据中 分离出来的一个个独立的信号源吧，这些信号源呢 对于我们想要分析的脑活动有的是有用的，有的是无用的（比如眼动和肌电等它们可能会影响我们的后续分析，因此需要把它们除掉），通过我们matlab中强大的工具箱SASICA可以快速的识别出哪些是有用成分，那些又是无用成分，进而将于我们无用的成分除掉

下一个疑问 为什莫得到的这些脑图的数量恰好等于我们数据的通道数呢，到底和通道有啥联系呢？？？

再问AI ... ...

为什莫这样巧合与通道数一模一样呢？？？

Soga！balabala...总之和我冥冥之中的方阵猜测契合啦！真相只有一个那就是ICA算法原理的某些数学特性需要的是一个个方阵，因此独立成分IC的数量就恰好与通道数一致啦！

· 平均重参考和坏电极插值

%% 平均重参考

EEG = pop_reref( EEG, [], 'refstate',-1);   %here is averaged referene, also can referenced to one certain channel
EEG = eeg_checkset(EEG);

%% 坏电极插值
EEG = pop_interp(EEG,originalEEG.chanlocs,'spherical');   % 之后对坏电极进行插值处理
EEG = eeg_checkset( EEG );
pop_eegplot(EEG, 1, 1, 1);

EEG = pop_saveset(EEG);                       % 最后处理完之后数据保存：S***_postICA.set

  • 功能：进行平均重参考和坏电极插值。

  • 解释：

    • pop_reref：平均重参考。

    • pop_interp：坏电极插值。

    • pop_saveset：保存最终处理后的数据。

数据对比

删除电极之后...

截段之后...

ICA之后...

评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值