台大 MEG 分析處理流程

這裡是選擇sub8exp2shortrun2_sss.fif資料進行示範

主要有四個步驟

  1. 根據R波位置切epoch,存成Epo.mat (其中要選比較乾淨的ECG channel)
  2. 去除jump和muscle artifact,存成Epo_noartifact.mat
  3. 做ICA,存成Epo_IC.mat
  4. 去除component,存成Epo_Rej.mat,內含4個檔案 (step 4做0~100hz的coherence,這個頻率範圍不要改)

(1) EOG與component的coherence metrix,存成eog_fdcomp.mat

(2) ECG與component的coherence metrix,存成ecg_fdcomp.mat

(3) 去除EOG與ECG的component,存成Epo_JMICA.mat

(4) 只去除EOG的component,存成Epo_JMICA_ECG.mat

第一步

流程說明

根據R波位置切epoch,存成Epo.mat (其中要選比較乾淨的ECG channel)

設定要選擇的實驗數據

subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%% short: 1,2   long: 1,2
ID = num2str(subjectID);
R = num2str(run);

載入實驗數據

這裡要注意一下

要把路徑位置改成自己本機的位置

%%% load fif file
try
%fname = sprintf(['/raid/smhsu/PA/MEG_data/sub' ID '/exp2' condition 'run' R '_sss.fif']);
if strcmp(condition,'nose') ==1 
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp1' condition '_sss.fif']);
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp1' condition '_sss.fif']);
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
end
end

預處理 MEG 訊號

這裡主要是將台大那邊給的訊號檔

利用ft_preprocessing進行處理的動作

%preprocessing
try
cfg = []; 
cfg.dataset    = fname;                                                              
data_orig = ft_preprocessing(cfg);

% ready triggers and redefine the onset and offset of the signal
event = ft_read_event(cfg.dataset); 
sample_onset = [event(find(strcmp('STI001', {event.type}))).sample]';
sample_offset = [event(find(strcmp('STI002', {event.type}))).sample]';

cfg = [];
cfg.toilim    = [(sample_onset+20000)/data_orig.fsample sample_offset/data_orig.fsample]; % discard the initial period, redefine a latency window in seconds
%%cfg.toilim    = [1/data_orig.fsample sample_offset/data_orig.fsample];% for sub2 longrun1 only
data_new = ft_redefinetrial(cfg, data_orig);

% %%%downsample to 500hz
% cfg = [];
% cfg.resamplefs = 500;
% cfg.detrend    = 'no';
% data_new_down = ft_resampledata(cfg, data_new);
end

載入 R 波索引位置

%%%%%%%%load ECG R wave location
try
fname2 = ['sub' ID '_' condition '_heart_rate'];
ftpath   = ['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '']; 
load(fullfile(ftpath, fname2));
if run ==1
    ECG_wave = ECG_wave1;
elseif run ==2
    ECG_wave = ECG_wave2;
end
end

重新轉換 ECG 訊號至原本的時間

% transform ECG_wave back to the original timeframe
ECG_wave_new = (sample_onset+20000)+ECG_wave.*2;

檢查 ECG 並選擇通道

ECG 通道主要位於底下 兩個

因為 MEG 處理出來的訊號

總共有 255 個

  1. ECG1=data_new.trial{1,1}(3,:);
  2. ECG2=data_new.trial{1,1}(4,:);
%%%check and plot, choose ECG channel
t=[180000:184000];   %%%%time (ms)
try
    a=find(ECG_wave_new<=t(end) & ECG_wave_new>=t(1));
    ECG1=data_new.trial{1,1}(3,:);
    ECG2=data_new.trial{1,1}(4,:);
    subplot(2,1,1); plot(t,ECG1(t-(sample_onset+20000))); hold on;  plot(ECG_wave_new(a),ECG1(a),'rs'); xlabel('ECG003');
    subplot(2,1,2); plot(t,ECG2(t-(sample_onset+20000))); hold on;  plot(ECG_wave_new(a),ECG2(a),'rs'); xlabel('ECG004');
end

Check_ECG_Channel

利用 ECG 的 R Peak 位置重新定義 MEG 訊號

這裡會將 MEG 訊號

透過 ECG 的 R Peak 進行切割

會分割成每一個 R Peak 就一個 MEG Data

這裡我們會定義副程式

來進行分割的動作

分割副程式

將該副程式檔名命名成trialfun_ECG.m

function [trl, event] = trialfun_ECG(cfg)

% read the header information and the events from the data
hdr   = ft_read_header(cfg.dataset);
event = ft_read_event(cfg.dataset);               %%%%%add

% %%%%%%%%%%%load
% if cfg.exp == 1 && cfg.condition == 1 && cfg.run ==0
%     load(['sub',num2str(cfg.ID),'_nose_idRespi.mat']);
% elseif cfg.exp == 1 && cfg.condition == 2 && cfg.run ==0
%     load(['sub',num2str(cfg.ID),'_mouth_idRespi.mat']);
% elseif cfg.exp == 2 && cfg.condition == 1 && cfg.run ==1
%     load(['sub',num2str(cfg.ID),'_short1_idRespi.mat']);
% elseif cfg.exp == 2 && cfg.condition == 1 && cfg.run ==2
%     load(['sub',num2str(cfg.ID),'_short2_idRespi.mat']);
% elseif cfg.exp == 2 && cfg.condition == 2 && cfg.run ==1
%     load(['sub',num2str(cfg.ID),'_long1_idRespi.mat']);
% elseif cfg.exp == 2 && cfg.condition == 2&& cfg.run ==2
%     load(['sub',num2str(cfg.ID),'_long2_idRespi.mat']);
% end

% determine the number of samples before and after the period of interest
pretrig  = -round(cfg.trialdef.prestim  * hdr.Fs); 
posttrig =  round(cfg.trialdef.poststim * hdr.Fs);

% trial epoching between two minimal points (one breath cycle)
trl = [];
for j = 1:(length(cfg.ECG_wave))
    if (j+1<=length(cfg.ECG_wave))       %%%%%%%%%%%must not exceed dimension
      trlbegin = cfg.ECG_wave(j+1) + pretrig; % j+1 to prevent no sufficient prestimulus period      
      trlend   = cfg.ECG_wave(j+1) + posttrig;       
      offset   = pretrig;
      newtrl   = [trlbegin trlend offset];
      trl      = [trl; newtrl];
    end
end

透過ft_definetrial來重新定義訊號

副程式寫完後

cfg.trialfun = 'trialfun_ECG';

並透過ft_definetrial來重新定義訊號

也就是進行訊號分割的動作

%% define trial (0.4s+0.7s)
try
cfg = [];
cfg.dataset     = fname;
%cfg.channel ='meg';           %%%%%%%%%
% cfg.channel ='all';   
% cfg.datafile = 'data_new_down';
% cfg.headerfile   = hdr;
cfg.trialfun = 'trialfun_ECG';
cfg.trialdef.prestim        = 0.4;% in seconds
cfg.trialdef.poststim       = 0.7; % in seconds
% cfg.ID=subjectID;               %%%%%%%%%%%%%%%%%%%%%%%%%%%add ID to cfg to trialfun
% cfg.exp=exp;
% cfg.condition=condition;
% cfg.run=run;
cfg.ECG_wave=ECG_wave_new;

cfg1 = ft_definetrial(cfg);
Epo=ft_preprocessing(cfg1);
end

檢查分割後訊號

這裡選擇 ECG 通道進行檢查的動作

%%%check and plot trial (check 0.4s+R wave +0.7s)
trial=145;    %%%%the trial to plot
cha=2;     %%%%1 or 2  (cha=1 -> ECG003  , cha=2  -> ECG004)
try
if cha==1
    plot(Epo.trial{1,trial}(3,:));
elseif cha==2
    plot(Epo.trial{1,trial}(4,:));
end 
end

Check_Cut_MEG_Data

儲存檔案為 Epo 檔

%%%save
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
end

%%save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');

第一步完整程式碼

%%ECG_step1
subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%% short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);
%%%load fif file
try
%fname = sprintf(['/raid/smhsu/PA/MEG_data/sub' ID '/exp2' condition 'run' R '_sss.fif']);
if strcmp(condition,'nose') ==1 
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp1' condition '_sss.fif']);
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp1' condition '_sss.fif']);
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '/exp2' condition 'run' R '_sss.fif']);
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
end
end

%% check ECG
%preprocessing
try
cfg = []; 
cfg.dataset    = fname;                                                              
data_orig = ft_preprocessing(cfg);

% ready triggers and redefine the onset and offset of the signal
event = ft_read_event(cfg.dataset); 
sample_onset = [event(find(strcmp('STI001', {event.type}))).sample]';
sample_offset = [event(find(strcmp('STI002', {event.type}))).sample]';

cfg = [];
cfg.toilim    = [(sample_onset+20000)/data_orig.fsample sample_offset/data_orig.fsample]; % discard the initial period, redefine a latency window in seconds
%%cfg.toilim    = [1/data_orig.fsample sample_offset/data_orig.fsample];% for sub2 longrun1 only
data_new = ft_redefinetrial(cfg, data_orig);

% %%%downsample to 500hz
% cfg = [];
% cfg.resamplefs = 500;
% cfg.detrend    = 'no';
% data_new_down = ft_resampledata(cfg, data_new);
end

%%%%%%%%load ECG R wave location
try
fname2 = ['sub' ID '_' condition '_heart_rate'];
ftpath   = ['C:/Users/king0/Documents/Matlab Code/MEG_Preprocess/sub' ID '']; 
load(fullfile(ftpath, fname2));
if run ==1
    ECG_wave = ECG_wave1;
elseif run ==2
    ECG_wave = ECG_wave2;
end
end

% transform ECG_wave back to the original timeframe
ECG_wave_new = (sample_onset+20000)+ECG_wave.*2;

%%%check and plot, choose ECG channel
t=[180000:184000];   %%%%time (ms)
try
    a=find(ECG_wave_new<=t(end) & ECG_wave_new>=t(1));
    ECG1=data_new.trial{1,1}(3,:);
    ECG2=data_new.trial{1,1}(4,:);
    subplot(2,1,1); plot(t,ECG1(t-(sample_onset+20000))); hold on;  plot(ECG_wave_new(a),ECG1(a),'rs'); xlabel('ECG003');
    subplot(2,1,2); plot(t,ECG2(t-(sample_onset+20000))); hold on;  plot(ECG_wave_new(a),ECG2(a),'rs'); xlabel('ECG004');
end

%% define trial (0.4s+0.7s)
try
cfg = [];
cfg.dataset     = fname;
%cfg.channel ='meg';           %%%%%%%%%
% cfg.channel ='all';   
% cfg.datafile = 'data_new_down';
% cfg.headerfile   = hdr;
cfg.trialfun = 'trialfun_ECG';
cfg.trialdef.prestim        = 0.4;% in seconds
cfg.trialdef.poststim       = 0.7; % in seconds
% cfg.ID=subjectID;               %%%%%%%%%%%%%%%%%%%%%%%%%%%add ID to cfg to trialfun
% cfg.exp=exp;
% cfg.condition=condition;
% cfg.run=run;
cfg.ECG_wave=ECG_wave_new;

cfg1 = ft_definetrial(cfg);
Epo=ft_preprocessing(cfg1);
end

%%%check and plot trial (check 0.4s+R wave +0.7s)
trial=145;    %%%%the trial to plot
cha=2;     %%%%1 or 2  (cha=1 -> ECG003  , cha=2  -> ECG004)
try
if cha==1
    plot(Epo.trial{1,trial}(3,:));
elseif cha==2
    plot(Epo.trial{1,trial}(4,:));
end 
end

%%%save
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
end

%%save(['sub' ID '_' condition 'run' R '_Epo'],'Epo');

第二步

流程說明

去除jump和muscle artifact,存成Epo_noartifact.mat

也就是進行去雜訊的動作

設定要選擇的實驗數據

%%ECG_step2
subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%%short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);

載入分割完的資料

這裡載入的是剛剛第一步分割完畢的資料

%%%load Epo
try
%load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
if strcmp(condition,'nose') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo']);
    load(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo']);
    load(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
end
end

移除 Jump 訊號

%%%%%%%%%%%% jump rejection   
try
cfg                    = [];
%cfg.trl = trl;                                                 
% cfg.continuous = 'yes';
% channel selection, cutoff and padding
cfg.artfctdef.zvalue.channel    = 'meg';
cfg.artfctdef.zvalue.cutoff     = 50;
cfg.artfctdef.zvalue.trlpadding = 0;
% cfg.artfctdef.zvalue.trlpadding = -0.2;
cfg.artfctdef.zvalue.artpadding = 0;
cfg.artfctdef.zvalue.fltpadding = 0;
 
% algorithmic parameters
cfg.artfctdef.zvalue.cumulative    = 'yes';
cfg.artfctdef.zvalue.medianfilter  = 'yes';
cfg.artfctdef.zvalue.medianfiltord = 9;
cfg.artfctdef.zvalue.absdiff       = 'yes';
 
% make the process interactive
cfg.artfctdef.zvalue.interactive = 'yes';
 
[cfg, Epo_J] = ft_artifact_zvalue(cfg, Epo);
end

該程式碼執行過後

會出現以圖形介面的方式

讓你進行移除訊號的動作

選擇欲移除的訊號片段

點擊reject full後 關閉該頁面

即可移除該訊號

remove_jump_trail

移除肌肉跳動的訊號

%%%%% muscle rejection
try
  % muscle
  cfg            = [];
  %cfg.trl        = trl;                      
%   cfg.continuous = 'yes'; 
 
  % channel selection, cutoff and padding
  cfg.artfctdef.zvalue.channel = 'meg';
  cfg.artfctdef.zvalue.cutoff      = 50;
  cfg.artfctdef.zvalue.trlpadding = 0;
%   cfg.artfctdef.zvalue.trlpadding  = -0.2;
  cfg.artfctdef.zvalue.fltpadding  = 0;
%   cfg.artfctdef.zvalue.artpadding  = 0.1;
 cfg.artfctdef.zvalue.artpadding  = 0;
 
  % algorithmic parameters
  cfg.artfctdef.zvalue.bpfilter    = 'yes';
  cfg.artfctdef.zvalue.bpfreq      = [110 140];
  cfg.artfctdef.zvalue.bpfiltord   = 9;
  cfg.artfctdef.zvalue.bpfilttype  = 'but';
  cfg.artfctdef.zvalue.hilbert     = 'yes';
  cfg.artfctdef.zvalue.boxcar      = 0.2;
 
  % make the process interactive
  cfg.artfctdef.zvalue.interactive = 'yes';
 
  [cfg, Epo_M] = ft_artifact_zvalue(cfg, Epo);
end

這一步跟上一步一樣

選擇你覺得像雜訊的波型

進行移除就行了

remove_muscle_trail

移除 jump 和 muscle

因為前兩步 是將雜訊位置儲存成變數

所以才是正式將原訊號進行移除雜訊的動作

透過artfctdef進行去除的動作

更詳細的方法可以參考官網

artfctdefhttp://www.fieldtriptoolbox.org/reference/ft_rejectartifact/

%%%%%reject jump and muscle
try
cfg=[]; 
% cfg.artfctdef.reject = 'partial';
cfg.artfctdef.jump.artifact = Epo_J;
cfg.artfctdef.muscle.artifact = Epo_M;
Epo_no_artifacts = ft_rejectartifact(cfg,Epo); 
end
%

當用圖形介面去除 jump 和 muscle 失敗時的方法

%%%%%%when jump and muscle can't work,  %%%%%%%%%%%%(very important)
Epo_no_artifacts=Epo;

%%%%%only keep ECG channel
try
cfg = [];                             
cfg.channel    = {'ECG003','ECG004'};  %MISC(respiration),STI001(exp onset), STI002(exp offset)                                    
Epo_temp = ft_preprocessing(cfg,Epo_no_artifacts);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

檢查訊號全部分割訊號

%%%%%%%%%%%%%%%%%%%check trial
cha=2;    %%%the same as step 1
%%%press > or < button to ECG003 or ECG004 channel to check every trial
try
%  cfg=[];
%  cfg.method      = 'trial';
%  if cha==1
%     cfg.channel = {'ECG003'};
% end
% if cha==2
%     cfg.channel = {'ECG004'};
% end
  cfg=[];
  cfg.method      = 'channel';
 ft_rejectvisual(cfg, Epo_no_artifacts)
end
%

這裡可以讓你概觀的檢查所有的分割訊號

如果你有不滿意的訊號

你可以紀錄 該訊號的編號

後面會寫個程式 來將其進行刪除的動作

Check_All_ECG_Data

刪除訊號片段

delete_trial裡面的陣列

填寫欲刪除的編號

%%%delete trial
delete_trial=[127];
try
trialsToKeep = 1:length(Epo_no_artifacts.trial);
trialsToKeep(delete_trial) = []; 
cfg = [];
cfg.trials = trialsToKeep;
Epo_no_artifacts = ft_preprocessing(cfg, Epo_no_artifacts);
end
%

儲存檔案

%%%save
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
end

完整程式碼

%%ECG_step2
subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%%short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);
%%%load Epo
try
%load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
if strcmp(condition,'nose') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo']);
    load(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo']);
    load(['sub' ID '_' condition '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo']);
    load(['sub' ID '_' condition 'run' R '_Epo'],'Epo');
end
end



%%%%%%%%%%%% jump rejection   
try
cfg                    = [];
%cfg.trl = trl;                                                 
% cfg.continuous = 'yes';
% channel selection, cutoff and padding
cfg.artfctdef.zvalue.channel    = 'meg';
cfg.artfctdef.zvalue.cutoff     = 50;
cfg.artfctdef.zvalue.trlpadding = 0;
% cfg.artfctdef.zvalue.trlpadding = -0.2;
cfg.artfctdef.zvalue.artpadding = 0;
cfg.artfctdef.zvalue.fltpadding = 0;
 
% algorithmic parameters
cfg.artfctdef.zvalue.cumulative    = 'yes';
cfg.artfctdef.zvalue.medianfilter  = 'yes';
cfg.artfctdef.zvalue.medianfiltord = 9;
cfg.artfctdef.zvalue.absdiff       = 'yes';
 
% make the process interactive
cfg.artfctdef.zvalue.interactive = 'yes';
 
[cfg, Epo_J] = ft_artifact_zvalue(cfg, Epo);
end

%%%%% muscle rejection
try
  % muscle
  cfg            = [];
  %cfg.trl        = trl;                      
%   cfg.continuous = 'yes'; 
 
  % channel selection, cutoff and padding
  cfg.artfctdef.zvalue.channel = 'meg';
  cfg.artfctdef.zvalue.cutoff      = 50;
  cfg.artfctdef.zvalue.trlpadding = 0;
%   cfg.artfctdef.zvalue.trlpadding  = -0.2;
  cfg.artfctdef.zvalue.fltpadding  = 0;
%   cfg.artfctdef.zvalue.artpadding  = 0.1;
 cfg.artfctdef.zvalue.artpadding  = 0;
 
  % algorithmic parameters
  cfg.artfctdef.zvalue.bpfilter    = 'yes';
  cfg.artfctdef.zvalue.bpfreq      = [110 140];
  cfg.artfctdef.zvalue.bpfiltord   = 9;
  cfg.artfctdef.zvalue.bpfilttype  = 'but';
  cfg.artfctdef.zvalue.hilbert     = 'yes';
  cfg.artfctdef.zvalue.boxcar      = 0.2;
 
  % make the process interactive
  cfg.artfctdef.zvalue.interactive = 'yes';
 
  [cfg, Epo_M] = ft_artifact_zvalue(cfg, Epo);
end
%

%%%%%reject jump and muscle
try
cfg=[]; 
% cfg.artfctdef.reject = 'partial';
cfg.artfctdef.jump.artifact = Epo_J;
 cfg.artfctdef.muscle.artifact = Epo_M;
Epo_no_artifacts = ft_rejectartifact(cfg,Epo); 
end
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%when jump and muscle can't work,  %%%%%%%%%%%%(very important)
Epo_no_artifacts=Epo;

%%%%%only keep ECG channel
try
cfg = [];                             
cfg.channel    = {'ECG003','ECG004'};  %MISC(respiration),STI001(exp onset), STI002(exp offset)                                    
Epo_temp = ft_preprocessing(cfg,Epo_no_artifacts);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%check trial
cha=2;    %%%the same as step 1
%%%press > or < button to ECG003 or ECG004 channel to check every trial
try
%  cfg=[];
%  cfg.method      = 'trial';
%  if cha==1
%     cfg.channel = {'ECG003'};
% end
% if cha==2
%     cfg.channel = {'ECG004'};
% end
  cfg=[];
  cfg.method      = 'channel';
 ft_rejectvisual(cfg, Epo_no_artifacts)
end
%

%%%delete trial
delete_trial=[127];
try
trialsToKeep = 1:length(Epo_no_artifacts.trial);
trialsToKeep(delete_trial) = []; 
cfg = [];
cfg.trials = trialsToKeep;
Epo_no_artifacts = ft_preprocessing(cfg, Epo_no_artifacts);
end
%


%%%save
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');
end

%%%save(['sub' ID '_' condition 'run' R '_Epo_noartifact'],'Epo_no_artifacts');

第三步

選擇資料

%%ECG_step2
subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%%short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);

載入資料

%%%load Epo_noartifact
try
%load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
if strcmp(condition,'nose') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
end
end

ICA 獨立成分分析

這裡主要將訊號重新取樣至 500

而我們在進行ICA前 需先設定 n_comp 參數

n_comp 參數為該訊號的個數或維度

%%%%%%%%%%%%% ICA
%%% take meg to compute n_comp,downsample
try
%%%% take meg to compute n_comp
cfg=[];
cfg.channel    = 'meg';                                   
Epo_meg = ft_preprocessing(cfg,Epo_no_artifacts);
n_comp = rank(squeeze(Epo_meg.trial{1}) * squeeze(Epo_meg.trial{1})');  
%%%%%%% downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
Epo_no_artifacts_down = ft_resampledata(cfg, Epo_no_artifacts);
end
%

獨立成分分析

% perform the independent component analysis (i.e., decompose the data)
try
cfg        = [];
cfg.channel = 'meg';
cfg.method = 'runica'; % this is the default and uses the implementation from EEGLAB
%n_comp = rank(squeeze(Epo_no_artifacts_down.trial{1}) * squeeze(Epo_no_artifacts_down.trial{1})');  
n_comp = rank(squeeze(Epo_meg.trial{1}) * squeeze(Epo_meg.trial{1})');
cfg.runica.pca = n_comp;
cfg.runica.stop = 1e-7;
comp = ft_componentanalysis(cfg, Epo_no_artifacts_down);
end
%

儲存檔案

%%%save comp
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo_IC'],'comp');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo_IC'],'comp');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
end

%%save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');

完整程式碼

%%ECG_step2
subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%%short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);
%%%load Epo_noartifact
try
%load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
if strcmp(condition,'nose') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
end
end


%%%%%%%%%%%%% ICA
%%%take meg to compute n_comp,downsample
try
%%%%take meg to compute n_comp
cfg=[];
cfg.channel    = 'meg';                                   
Epo_meg = ft_preprocessing(cfg,Epo_no_artifacts);
n_comp = rank(squeeze(Epo_meg.trial{1}) * squeeze(Epo_meg.trial{1})');  
%%%%%%%downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
Epo_no_artifacts_down = ft_resampledata(cfg, Epo_no_artifacts);
end
%

% perform the independent component analysis (i.e., decompose the data)
try
cfg        = [];
cfg.channel = 'meg';
cfg.method = 'runica'; % this is the default and uses the implementation from EEGLAB
%n_comp = rank(squeeze(Epo_no_artifacts_down.trial{1}) * squeeze(Epo_no_artifacts_down.trial{1})');  
n_comp = rank(squeeze(Epo_meg.trial{1}) * squeeze(Epo_meg.trial{1})');
cfg.runica.pca = n_comp;
cfg.runica.stop = 1e-7;
comp = ft_componentanalysis(cfg, Epo_no_artifacts_down);
end
%

%%%save comp
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo_IC'],'comp');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo_IC'],'comp');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');
end

%%save(['sub' ID '_' condition 'run' R '_Epo_IC'],'comp');

第四步

選擇資料

subjectID=8;
condition='short';   %%%% nose , mouth , short , long
run=2;              %%%%short: 1,2   long: 1,2

ID = num2str(subjectID);
R = num2str(run);

載入資料

這裡主要是載入

去除雜訊、做完ICA、原始訊號檔

%%%load Epo_noarttifact  and Epo_IC and fif file
try
%load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
%load(['sub' ID '_' condition 'run' R '_Epo_IC']);
if strcmp(condition,'nose') ==1 
    fname = sprintf(['D:/work/pa/sub' ID '/exp1' condition '_sss.fif']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_IC']);
elseif strcmp(condition,'mouth') ==1 
    fname = sprintf(['D:/work/pa/sub' ID '/exp1' condition '_sss.fif']);
    load(['sub' ID '_' condition '_Epo_noartifact']);
    load(['sub' ID '_' condition '_Epo_IC']);
elseif strcmp(condition,'short') ==1 && run ==1
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_IC']);
elseif strcmp(condition,'short') ==1 && run ==2
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_IC']);
elseif  strcmp(condition,'long') ==1 && run ==1
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_IC']);
elseif  strcmp(condition,'long') ==1 && run ==2
    fname = sprintf(['D:/work/pa/sub' ID '/exp2' condition 'run' R '_sss.fif']);
    load(['sub' ID '_' condition 'run' R '_Epo_noartifact']);
    load(['sub' ID '_' condition 'run' R '_Epo_IC']);
end
end

畫出腦磁圖

%%%%%%plot the components for visual inspection (modify or default) (picture)
figure;
cfg = [];
cfg.component = [1:30];     %(modify or default) specify the component(s) that should be plotted
cfg.layout    = 'neuromag306all.lay'; % specify the layout file that should be used for plotting
cfg.comment   = 'no';
ft_topoplotIC(cfg, comp)

plot_meg_img

重新處理訊號

%%%data
try
cfg = []; 
cfg.dataset    = fname;                               
cfg.channel    = {'MISC001','STI*','ecg','eog','emg'};%MISC(respiration),STI001(exp onset), STI002(exp offset)                                    
data_orig = ft_preprocessing(cfg);

% ready triggers and redefine the onset and offset of the signal
event = ft_read_event(cfg.dataset); 
sample_onset = [event(find(strcmp('STI001', {event.type}))).sample]';
sample_offset = [event(find(strcmp('STI002', {event.type}))).sample]';

cfg = [];
cfg.toilim    = [(sample_onset+20000)/data_orig.fsample sample_offset/data_orig.fsample]; % discard the initial period, redefine a latency window in seconds
data_new = ft_redefinetrial(cfg, data_orig);
end
%

畫出磁強計

%%%%%%visualization of magnetometers (auto) (picture)
try
cfg = [];
cfg.layout = 'neuromag306mag.lay'; % specify the layout file that should be used for plotting
cfg.viewmode = 'component';
ft_databrowser(cfg, comp)
end
%

magnetometers

檢查訊號

這裡可以自行選擇片段來進行檢查

可以檢查 EOG EOG ECG ECG ECG

%%%choose trial to check
trial=14;   %%%(modify)
try
figure; 
a=Epo_no_artifacts.sampleinfo(trial,1)-data_new.sampleinfo(1)+1;
b=Epo_no_artifacts.sampleinfo(trial,2)-data_new.sampleinfo(1)+1;

subplot(4,1,1), plot(data_new.trial{1,1}(1,a:b)); title('EOG'); %%%%%EOG EOG ECG ECG ECG
subplot(4,1,2), plot(data_new.trial{1,1}(2,a:b)); title('EOG');
subplot(4,1,3), plot(data_new.trial{1,1}(3,a:b)); title('ECG');
subplot(4,1,4), plot(data_new.trial{1,1}(4,a:b)); title('ECG');
end
%

Check_4_channel_data

畫出 gradiometers

%%%%%%visualization of gradiometers (auto) (picture)
try
cfg = [];
cfg.layout = 'neuromag306planar.lay'; % specify the layout file that should be used for plotting
cfg.viewmode = 'component';
ft_databrowser(cfg, comp)
end
%

gradiometers

compute coherence between all components and the EOG (around 0~100hz)

%%%%compute coherence between all components and the EOG (around 0~100hz)
cha=2;    %%%%1 or 2  (cha=1 -> EOG001  , cha=2  -> EOG002)
try
cfg=[];
if cha==1
    cfg.channel    = {'EOG001'};
elseif cha==2
    cfg.channel    = {'EOG002'};
end
eog            = ft_preprocessing(cfg,Epo_no_artifacts);

%%%%downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
eog_down = ft_resampledata(cfg, eog);

% append the eog channel to the data structure
comp_eog      = ft_appenddata([], eog_down, comp);

% compute a frequency decomposition of all components and the EOG
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.foilim     = [0 100];   %%%%%%%%%hz
cfg.taper      = 'hanning';
cfg.pad        = 'maxperlen';
freq           = ft_freqanalysis(cfg, comp_eog);

% compute coherence between all components and the EOG
cfg            = [];
cfg.channelcmb = {'all' 'eog'};
cfg.jackknife  = 'no';
cfg.method     = 'coh';
eog_fdcomp         = ft_connectivityanalysis(cfg, freq);
end
%

look at the coherence spectrum between all components and the EOG

figure;
subplot(2,1,1); plot(eog_fdcomp.freq, abs(eog_fdcomp.cohspctrm)); xlabel('Hz'); ylabel('coherence');
subplot(2,1,2); imagesc(eog_fdcomp.freq,[1:length(eog_fdcomp.labelcmb)],abs(eog_fdcomp.cohspctrm)); xlabel('Hz'); ylabel('component');
eog_comp_avg=mean(eog_fdcomp.cohspctrm(:,7),2);   %%%0~5.44hz avg
[value,probabaly_eog_comp]=max(eog_comp_avg)    %%%find component of max value

EOG_coherence

compute coherence between all components and the ECG (around 0~100hz)

%%%%compute coherence between all components and the ECG (around 0~100hz)
cha=2;    %%%the same as step 1
try
cfg=[];
if cha==1
    cfg.channel    = {'ECG003'};
elseif cha==2
    cfg.channel    = {'ECG004'};
end
ecg            = ft_preprocessing(cfg,Epo_no_artifacts);

%%%%downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
ecg_down = ft_resampledata(cfg, ecg);

% append the ecg channel to the data structure
comp_ecg      = ft_appenddata([], ecg_down, comp);

% compute a frequency decomposition of all components and the ECG
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.foilim     = [0 100];    %%%%hz
cfg.taper      = 'hanning';
cfg.pad        = 'maxperlen';
freq           = ft_freqanalysis(cfg, comp_ecg);

% compute coherence between all components and the ECG
cfg            = [];
cfg.channelcmb = {'all' 'ecg'};
cfg.jackknife  = 'no';
cfg.method     = 'coh';
ecg_fdcomp         = ft_connectivityanalysis(cfg, freq);
end
%

look at the coherence spectrum between all components and the ECG

figure;
subplot(2,1,1); plot(ecg_fdcomp.freq, abs(ecg_fdcomp.cohspctrm)); xlabel('Hz'); ylabel('coherence');
subplot(2,1,2); imagesc(ecg_fdcomp.freq,[1:length(ecg_fdcomp.labelcmb)],abs(ecg_fdcomp.cohspctrm)); xlabel('Hz'); ylabel('component');
ecg_comp_avg=mean(ecg_fdcomp.cohspctrm(:,35),2);  %%%0~30.85hz avg
[value,probabaly_ecg_comp]=max(ecg_comp_avg)   %%%find component of max value

ECG_coherence

刪除 EOG 和 ECG component

%%%%%%delete EOG and ECG component (modify) (very important)%%%%%%%%%
comp_delete=[1 37];  %(modify)modify which component to delete  (small to big)
try
%%%%%%%downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
Epo_no_artifacts_down = ft_resampledata(cfg, Epo_no_artifacts);
%%%delete
cfg = [];
cfg.component = comp_delete; %%
Epo_JMICA = ft_rejectcomponent(cfg, comp, Epo_no_artifacts_down);
end
%

%%%%%%only delete EOG component  (modify) (very important)%%%%%%%%%
comp_delete=[1];  %(modify)modify which component to delete  (small to big)
try
%%%%%%%downsample
cfg = [];
cfg.resamplefs = 500;
cfg.detrend    = 'no';
Epo_no_artifacts_down = ft_resampledata(cfg, Epo_no_artifacts);
%%%delete
cfg = [];
cfg.component = comp_delete; %%
Epo_JMICA_ECG = ft_rejectcomponent(cfg, comp, Epo_no_artifacts_down);
end
%

儲存Epo、EOG coherence、 ECG coherence

%%%%save Epo, EOG coherence, ECG coherence
if strcmp(condition,'nose') ==1 
    save(['sub' ID '_' condition '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
elseif strcmp(condition,'mouth') ==1 
    save(['sub' ID '_' condition '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
elseif strcmp(condition,'short') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
elseif strcmp(condition,'short') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
elseif  strcmp(condition,'long') ==1 && run ==1
    save(['sub' ID '_' condition 'run' R '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
elseif  strcmp(condition,'long') ==1 && run ==2
    save(['sub' ID '_' condition 'run' R '_Epo_Rej'],'Epo_JMICA','Epo_JMICA_ECG','eog_fdcomp','ecg_fdcomp');
end
Last modification:March 13, 2020
If you think my article is useful to you, please feel free to appreciate