matlab ismember_基于MATLAB源码逐行超详解ERP Decoding

前面曾经写过一个基于Python进行脑电Decoding的教程：

教你进行认知神经科学领域中的EEG Decoding研究(基于Python实现)

由于现在MATLAB仍然在Neuroscience领域属于一个使用率更高的编程语言，今天用一段开源代码超超超详解基于MATLAB的ERP Decoding的实现。

所谓超超超详解就是路同学会逐行中文注释解说每一行代码，希望自己完成的这个艰巨的工作能对大家起到帮助！

源码来自Bae&Luck于2018年发表在Journal of Neuroscience上的一篇文章《Dissociable Decoding of Spatial Attention and Working Memory from EEG Oscillations and Sustained Potentials》，关于实验方面的问题可以参考这篇文章，代码用于实现文中Experiment 2的ERP Decoding部分。

% 该源码来自于Bae & Luck 2018 Journal of Neuroscience的% Experiment 2用于实现ERP Decoding的部分% 在其基础上进行了中文详细注解function SVM_ECOC_ERP_Decoding(subs)% 关闭并行计算delete(gcp)parpoolif nargin ==0    % 被试id    subs = [201 202 203 204 205 206 207 208 209 210 212 213 215 216 217 218];       end% 获取被试数量nSubs = length(subs);%% 设置变量% 分类的条件数svmECOC.nChans = 16; svmECOC.nBins = svmECOC.nChans; % 迭代次数svmECOC.nIter = 10;% 交叉验证用的Blocks数svmECOC.nBlocks = 3; % 低通滤波参数svmECOC.frequencies = [0 6];% 用于分析的时间段% -500ms到1496ms% 间隔20ms一次采样svmECOC.time = -500:20:1496; % 1个时间点代表1+4个采样点(前后各2)的时间窗svmECOC.window = 4;% 采样率svmECOC.Fs = 250;% 16个用于分析的导联对应数据中的编号ReleventChan = sort([2,3,4,18,19, 5,6,20, 7,8,21, 9,10,11,12,13,14, 22,23,24,25,26,27, 15,16,1,17]);% 相关导联的数量svmECOC.nElectrodes = length(ReleventChan);%% 简化变量nChans = svmECOC.nChans;nBins = svmECOC.nBins;nIter = svmECOC.nIter;nBlocks = svmECOC.nBlocks;freqs = svmECOC.frequencies;times = svmECOC.time;nElectrodes = svmECOC.nElectrodes;nSamps = length(svmECOC.time);Fs = svmECOC.Fs;%% 遍历两种条件for cond = 1:2    % 1代表位置, 2代表朝向    %% 遍历所有被试for s = 1:nSubs    % 提取第s个被试的编号    sn = subs(s);    fprintf('Subject:\t%d\n',sn)    %% 读取数据    currentSub = num2str(sn); % 被试编号的类型变成字符串    dataLocation = pwd; % 设置数据存储的路径    % 完整的该被试的数据地址    loadThis = strcat(dataLocation,'/Decoding_OL_',currentSub,'.mat');    % 加载数据    load(loadThis)        saveLocation = pwd; % 设置存储路径        %% 设置每一个试次的位置/朝向的信息        if cond ==1    channel = data.targetLocBin; % 获取数据的位置信息    else    channel = data.targetOriBin; % 获取数据的朝向信息    end        % 对channel转置并存为posBin作为特征信息    svmECOC.posBin = channel';    posBin = svmECOC.posBin;        % 提取用于分析的16个条件的eeg数据    eegs = data.eeg(:,ReleventChan,:);         %% 设置时间点        % 将数据记录的时间点与需要用于分析的时间点做匹配    % 数据时间点序列上有匹配的时间点为1，无匹配的时间点为0    tois = ismember(data.time.pre:4:data.time.post,svmECOC.time);     % 数据时间点的数量    nTimes = length(tois);        % 试次数    svmECOC.nTrials = length(posBin);    nTrials = svmECOC.nTrials;     %% 预设矩阵    % 用来存储SVM的预测结果的矩阵svm_predict    % [迭代次数, 样本数, block数, 条件数]    svm_predict = nan(nIter,nSamps,nBlocks,nChans);    % 用来存储真实的目标值的矩阵tst_target    tst_target = nan(nIter,nSamps,nBlocks,nChans);  % a matrix to save true target values    % 存储block的分配信息    svmECOC.blocks = nan(nTrials,nIter);  % create svmECOC.block to save block assignments    %% 低通滤波        % 初始化矩阵用来存储滤波后的数据    filtData = nan(nTrials,nElectrodes,nTimes);    % 并行对每一个导联进行滤波    parfor c = 1:nElectrodes                      filtData(:,c,:) = eegfilt(squeeze(eegs(:,c,:)),Fs,freqs(1,1),freqs(1,2)); % low pass filter    end    %% 遍历迭代        tic % 开始计时    for iter = 1:nIter        % 预设置两个数组，长度为试次数        blocks = nan(size(posBin));        shuffBlocks = nan(size(posBin));        clear binCnt        % 确认每一种条件的试次数        for bin = 1:nBins            binCnt(bin) = sum(posBin == bin);         end        % 最小的条件试次数        minCnt = min(binCnt);        % 每一个block的最小的条件试次数(向下取整)        nPerBin = floor(minCnt/nBlocks);        % 打乱试次        shuffInd = randperm(nTrials)'; % 随机排列试次id        shuffBin = posBin(shuffInd); % 随后试次对应的特征信息        %% 获取每一种条件的试次数据        %% 生成每种条件试次一致的label                for bin = 1:nBins               % 获取条件bin对应的试次位置            idx = find(shuffBin == bin);             % 剔除多余的试次            % 只保留nPerBin*nBlocks这么多的试次            % 这是为了保证每一种条件的试次数一致            idx = idx(1:nPerBin*nBlocks);            % 生成[1,nBlocks]个相同条件的labels            x = repmat((1:nBlocks)',nPerBin,1);            % 把label传入对应的试次            shuffBlocks(idx) = x;        end        % 将未打乱的信息存于blocks矩阵        blocks(shuffInd) = shuffBlocks;        svmECOC.blocks(:,iter) = blocks; % block assignment        svmECOC.nTrialsPerBlock = length(blocks(blocks == 1));        % Average data for each position bin across blocks                posBins = 1:nBins;        % 定义一个矩阵用来存储对滤波后数据进行50hz重采样后的数据        blockDat_filtData = nan(nBins*nBlocks,nElectrodes,nSamps);        % 定义一个矩阵用于存储labels        labels = nan(nBins*nBlocks,1);        % 定义一个矩阵用于存储block标号        blockNum = nan(nBins*nBlocks,1);                                % block numbers for averaged & filtered EEG data        bCnt = 1;        % 向blockDat_filtData矩阵传入数据        for ii = 1:nBins            for iii = 1:nBlocks                blockDat_filtData(bCnt,:,:) = squeeze(mean(filtData(posBin==posBins(ii) & blocks==iii,:,tois),1));                labels(bCnt) = ii;                blockNum(bCnt) = iii;                bCnt = bCnt+1;            end        end        %% 对每个时间点进行SVM分类                parfor t = 1:nSamps            % 提取时间点t的数据            % 这里的t是按照采样点而非实际时间取            % 提取的是一个时间窗的数据(5个时间点)            toi = ismember(times,times(t)-svmECOC.window/2:times(t)+svmECOC.window/2);                        % 平均时间点的数据            dataAtTimeT = squeeze(mean(blockDat_filtData(:,:,toi),3));              % 对每一个block迭代            % 每一个block都作为一次测试集            for i=1:nBlocks                % 训练labels                trnl = labels(blockNum~=i);                % 测试labels                tstl = labels(blockNum==i);                % 训练数据                trnD = dataAtTimeT(blockNum~=i,:);                % 测试数据                tstD = dataAtTimeT(blockNum==i,:);    % test data                % 训练SVM                mdl = fitcecoc(trnD,trnl, 'Coding','onevsall','Learners','SVM' );   %train support vector mahcine                % 用训练好的分类器进行预测(测试)                LabelPredicted = predict(mdl, tstD);                % 保存预测得到的labels                svm_predict(iter,t,i,:) = LabelPredicted;                % 保存实际数据的labels                tst_target(iter,t,i,:) = tstl;            end        end    end    % 结束计时    toc        % 存储对位置decoding的结果    if cond ==1    OutputfName = strcat(saveLocation,'/Location_Results_ERPbased_',currentSub,'.mat');    % 存储对朝向decoding的结果    elseif cond ==2    OutputfName = strcat(saveLocation,'/Orientation_Results_ERPbased_',currentSub,'.mat');    end    svmECOC.targets = tst_target;    svmECOC.modelPredict = svm_predict;     svmECOC.nBlocks = nBlocks;    save(OutputfName,'svmECOC','-v7.3');endend

原版源码及数据可见：

https://osf.io/bpexa/

路同学修改并添加超详细的中文注释的版本可见GItHub：

https://github.com/ZitongLu1996/ERP_Decoding_withChineseAnnotations

记得点关注、点在看哦！也欢迎给作者打赏！

/END/

zitonglu1996.github.io

爱胡思乱想

渴望成为优秀学生的

路同学愿和你谈谈心。

路同学