multimoto_ica_pca_process_program.m

用于脑电分析的各类源码

%% 四种想像动作的EEG信号处理 ——  （基于ICA分析）
   % 四种动作分别为左手，右手，脚，舌头
   
   % 数据文件共三个，指定本次处理的文件名及路径
   path='F:\matlab_work\process_for_BCICOMP\for_C2005_dataset_iii\for_IIIa\k3b\'
   
%% 格式转换，将每导数据存为一个文件

    % 将原始数据由数百兆的一个大文件转换成每导数据一个文件
    for i=1:60
        data=s(:,i)';
        filename=['s',num2str(i),'.mat'];
        eval(['save ',path,filename,' data;']);
        clear data;
    end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 原始数据文件过大（约数百兆），此处理为减轻每次计算的内存负担    %    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    
%% 确定数据参数及时频参数
    eval(['load ',path,'HDR.mat;']);
    trig=HDR.TRIG;
    label=HDR.Classlabel;
    artifact=HDR.ArtifactSelection;
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 经检测各subject数据中噪声过大和含有NaN的trails的数量，分别为： %
        % subject:noise,NaN____k3b:62,2;    k6b:65,10;  l1b:73,0;      % 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %将噪声过大的数据段去除
        label(find(artifact==1))=-1;
        for i=1:4
            trig_s{i}=trig(find(label==i));
        end
        %将含NaN数值的数据段去除
        eval(['load ',path,'NaN_label.mat;']);
        label(find(NaN_label==1))=-1;
        for i=1:4
            trig_s{i}=trig(find(label==i));
        end
        
    % 设定计算参数
    fs=250;trail_t=8;frq=[8:12,20:24];

%%  确定各动作的ERD/ERS发生频段及显著导联位置
        %计算所有导联各频段PSD曲线，并将结果存储为中间数据

for k=1:60    
    % 输入单导数据
    filename=['s',num2str(k),'.mat'];
    eval(['load ',path,filename,';']);
    
    % 对四种动作分别将其PSD曲线处理
    for i=1:4
        trig_temp=trig_s{i};
        psd_motor=zeros(length(frq),1751);   % 只提取frq定义频段的psd数据
        % processing one trail by one trail
        for j=1:length(trig_temp)
            temp_index=trig_temp(j)+[0:fs*trail_t-1];
            tempdata=data(temp_index);
            
            temp_psd=spectrogram(tempdata,kaiser(fs,2),fs-1,fs,fs);
            temp_psd=abs(temp_psd(frq+1,:));
            psd_motor=psd_motor+temp_psd;
        end
        psd_motor=psd_motor/length(trig_temp);
        mi_chan_psd{k,i}=psd_motor;i
    end
    k
end
eval(['save ',path,'psd_aver_chan_frq.mat mi_chan_psd;']);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 本步处理较为费时，故将结果作为中间数据进行存储，以备调用        % 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 绘制各动作的frq段功率谱变化ERD/ERS系数地形图
    % 使用上述数据直接绘制与使用ICA＋PCA分解后的特征分量绘制的地形图进行对比

    % 各1hz频段的psd曲线作ica＋pca分解后的ERD/ERS系数
    ref_t=[1:1.5*fs]; erds_t=[1:1.5*fs]+3*fs; 
for i=1:4   % single moto
    for j=1:60  %single chan
        temp=mi_chan_psd{j,i};
        for k=1:length(frq)     %single hz frequency   
            psd_temp=temp(k,:);
            psd_ref=psd_temp(ref_t);
            psd_erds=psd_temp(erds_t);
            erds(j,k)=(mean(psd_erds)-mean(psd_ref))/abs(mean(psd_ref));
        end
    end
    % 绘出每种动作各频率段的ERD/ERS系数地形图
    q_multi_topo(erds,loc_leads);
end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 此处未使用图形存储命令,可手动加注图形说明并存储此4幅图形        %
        % 11－12Hz频段处的左右手ERD现象与脚舌部的ERS现象较为明显         %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  以ERD/ERS区别最明显的11－12Hz段psd值绘制连续显示的地形图动画
for i=1:4   % single moto
    for j=1:60  %single chan
        temp=mi_chan_psd{j,i};
        psd_temp(j,:)=mean(temp([5,6],:));   
    end
    eval(['psd_',num2str(i),'=psd_temp;']);
end

cmax=max(max([psd_1;psd_2;psd_3;psd_4]));
cmin=min(min([psd_1;psd_2;psd_3;psd_4]));
title_char={'左手';'右手';'脚部';'舌头'};
for i=1:175
    for j=1:4
        subplot(2,2,j)
        title(title_char{j});
        eval(['temp=psd_',num2str(j),';']);
        plotdata=temp(:,1+10*(i-1));
        [plotdata,loc_lead_temp]=q_nan_del(plotdata,loc_leads);
        q_topoplot(plotdata,loc_lead_temp,'label','char');
        caxis([cmin cmax]);colorbar;
    end    
    F(i)=getframe(gcf);
    i
end
save f_invers.mat F;
movie(F,1,1)
movie2avi(F,'4to1.avi','fps',25);

%%  绘制各动作的frq段功率谱变化ICA+PCA分解后的ERD/ERS分量地形图

% pinciple component numbers
pca_num=3;

% creat pic saving directory
pic_dic_name=['topo_pic\'];
mkdir([path,pic_dic_name(1:length(pic_dic_name)-1)]);

% processing
for i=1:4   % single moto
    for j=1:length(frq)  % single hz frequency
        for k=1:60     % single chan   
            temp=mi_chan_psd{k,i};
            psd_temp(k,:)=temp(j,:);  
        end

        % run ica+pca process and plot components curves and topograph
        [ww,ss]=runica(psd_temp,'pca',pca_num,'extended',-1);
        for m=1:pca_num
            psd(m,:)=q_norm(ww(m,:)*ss*psd_temp)+3-m;
        end
        
        % plotting
        ha=figure;
        subplot(2,1,1);plot(0:1/250:7,psd);
        for m=1:pca_num
            subplot(2,pca_num,pca_num+m);
            q_topoplot(ww(m,:),loc_leads,'label','char');
            colorbar;
        end
        
        % save pic
        pic_name=['moto_',num2str(i),'_psd_',num2str(frq(j)),'hz_ica_topo.jpg'];
        eval(['saveas(ha,''',path,pic_dic_name,pic_name,''');']);
        delete(ha);
    end
end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 人工选取具有明显erd及ers现象的频段                             %
        % 对比以上两种地形图的相似与区别                                 % 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  以ERD/ERS区别最明显的11－12Hz段psd曲线进行ICA+PCA处理，

% pinciple component numbers
pca_num=3;

% ica+pca processing
for i=1:4   % single moto
    for j=1:60  %single chan
        temp=mi_chan_psd{j,i};
        psd_temp(j,:)=mean(temp([5,6],:));   
    end
    clear psd;
    [ww,ss]=runica(psd_temp,'pca',pca_num,'extended',-1);
    for m=1:pca_num
        psd(m,:)=q_norm(ww(m,:)*ss*psd_temp)+3-m;
        coeff(i,m)=q_rattio_erds(psd(m,:),fs,[0 2],[3.5 6]);
    end
    ww_model{i,1}=ww;    ww_model{i,2}=ss;

    % plotting
    ha=figure;
    subplot(2,1,1);plot(0:1/250:7,psd);
    for m=1:pca_num
        subplot(2,pca_num,pca_num+m);
        q_topoplot(ww(m,:),loc_leads,'label','char');
        colorbar;
    end
end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 各动作明显的响应分量为：                                      % 
        %左手：第一分量ERD；右手：第二分量ERD；                          %
        %脚：第三分量ERD；舌头：第一分量ERS                              % 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  使用frq段所有trail的平均psd曲线进行ica＋pca处理
clear psd_temp;
for i=1:4   % single moto
    for j=1:60  %single chan
        temp=mi_chan_psd{j,i};
        psd_temp(j,:)=mean(temp);
    end
    % run ica+pca process and plot components curves and topograph
    pca_num=3;clear psd;
    [ww,ss]=runica(psd_temp,'pca',pca_num,'extended',-1);
    for m=1:pca_num
        psd(m,:)=q_norm(ww(m,:)*ss*psd_temp)+3-m;
    end
    % plotting
    ha=figure;
    subplot(2,1,1);plot(0:1/250:7,psd);
    for m=1:pca_num
        subplot(2,pca_num,pca_num+m);
        q_topoplot(ww(m,:),loc_leads,'label','char');
        colorbar;
    end

end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 上两个频段比较，11－12Hz段效果明显好于frq段效果，故取11－12Hz  %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 将各动作对应的分量权值取出，作为标准权重模板
component_no=[1,2,3,1];
for i=1:4
    ww_temp=ww_model{i,1};
    w_std(i,:)=ww_temp(component_no(i),:);
end

%% 使用ICA滤波处理单trail数据,并将滤波后数据存储以备下一步进行识别之用

data_ica_dic=['data_ica_filted\'];
mkdir([path,data_ica_dic(1:length(data_ica_dic)-1)]);

test_index=find(isnan(label)==1);
for i=1:length(test_index)
    trig_temp=trig(test_index(i));
    data_index=trig_temp+[0:fs*trail_t-1];
    clear test_data;
    for k=1:60    
        % 输入单导数据
        filename=['s',num2str(k),'.mat'];
        eval(['load ',path,filename,';']);
        test_data(k,:)=data(data_index);
    end
    i
    % filter data using ica
    kept_c_num=20; % kept components number
    test_data=q_ica_filter(test_data,trail_t,fs,frq,kept_c_num);
    
    % save data filted
    eval(['save ',path,data_ica_dic,'test_data_ica_filted_trail_',num2str(i),'.mat test_data;']);
end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % ica滤波参数可以调整，中间数据按照单个trail一个文件的格式存储    %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 对test数据进行识别，直接使用11－12hz频段psd曲线计算ERD/ERS系数作为特征值

last_frq=[11,12];fs=250;

for i=1:149
    eval(['load ',path,data_ica_dic,'test_data_ica_filted_trail_',num2str(i),'.mat;']);
    for j=1:60
        psd_temp(j,:)=q_psfd_cal(test_data(j,:),last_frq,fs);
        coef(i,j)=q_rattio_erds(psd_temp(j,:),fs,[0 2],[3.5 5.5]);
    end    
    i
end

% test数据的标准值
load true_label_k3.txt  % 将具体路径加入
test_index=find(isnan(label)==1);
true_label=true_label_k3(test_index);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 此处未考虑噪声过大的test数据                                  %
        % 直接使用各导处的11－12hz的erds系数作为参数                    %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 进行识别
for i=1:149
    for j=1:4
        temp1=q_norm(w_std(j,:));
        temp2=q_norm(coef(i,:)');
        temp3(j)=abs(temp1*temp2);
    end
    results(i)=find(temp3==max(temp3));
end
correct_rattio=sum(results'==true_label)/149
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 此处未考虑噪声过大的test数据                                  %
        % 该识别方法的具体参数定义和选择需进一步调整                     %
        % 最好使用交叉验证的方式来确定最优参数                           %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 简单的距离识别方法以训练数组的11－12Hz的ERD/ERS系数为标准样板集对test数据进行识别

% 训练样本集的ERD/ERS系数
    
last_frq=[11,12];
for k=1:60    
    % 输入单导数据
    filename=['s',num2str(k),'.mat'];
    eval(['load ',path,filename,';']);
    
    % 对四种动作分别将其PSD曲线处理
    for i=1:4
        trig_temp=trig_s{i};
        clear coef_temp;
        % processing one trail by one trail
        for j=1:length(trig_temp)
            temp_index=trig_temp(j)+[0:fs*trail_t-1];
            tempdata=data(temp_index);
            psd_temp=q_psfd_cal(tempdata,last_frq,fs);
            coef_temp(j)=q_rattio_erds(psd_temp,fs,[0 2],[3.5 5.5]);
        end
        coef_train{k,i}=coef_temp;
    end
    k
end
    % ERD/ERS系数:coef_tr及标准动作标志:std_group
coef_tr=[];std_group=[];
for i=1:4
    clear temp2;
    for j=1:60
        temp1=coef_train{j,i};
        temp2(j,:)=temp1;
    end
    coef_tr=[coef_tr;temp2'];
    std_group=[std_group;ones(size(temp2,2),1)*i];
end

% 分类识别
[class,err,posterior,logp] = classify(coef,coef_tr,std_group,'diagquadratic') ;
correct_rattio=sum(class==true_label)/149
sum(class==true_label)

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 注 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 本节作为简单的距离识别算法示例，可作为前面识别算法的对照组      %
        % classify 函数来自matlab自带的statistic toolbox                %
        % 使用何种距离定义可以在classify函数给定的类型中选择             %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%