ICA方法用于脑电图分析的程序

%% 四种想像动作的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 % %%%%%%%%%%%%%%%%%%%