%% 四種想像動作的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函數給定的類型中選擇             %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%