文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。
上一期的文章TMS-EEG数据处理教程（上）中详细地介绍了TMS伪影类型和预处理步骤。这期主要讲了完成数据预处理后，再进行一些（后）处理步骤，如过滤、去趋势、去均值和降采样。但要注意的是，一些分析步骤可能需要对数据进行不同的处理。例如，当查看经颅磁刺激诱发电位（TEPs）时，你可能想要滤除数据中的高频噪声，但在执行时频分析时（滤除高频噪声）是不必要的；你可能也希望对数据进行去趋势操作，但这同样不建议用于分析TEPs。
幸运的是，这些函数（例如ft_timelockanalysis和ft_freqanalysis）允许对输入数据进行与ft_preprocessing相同的预处理步骤。这适用于对每步分析进行单独的处理，而不必创建额外的数据结构。
这里先只对数据进行降采样，接下来会应用到不同的（后）处理步骤。降采样需要使用ft_resampledata。
cfg = [];
cfg.resamplefs = 1000;cfg.detrend = 'no';
cfg.demean = 'yes';
% Prior to downsampling a lowpass filter is applied, demeaning avoids artifacts at the edges of your trialdata_tms_clean = ft_resampledata(cfg, data_tms_clean);
接下来先需要将数据进行保存。
save('data_tms_clean','data_tms_clean','-v7.3');
分析
我们最初的问题是预收缩是否会影响经颅磁刺激诱发电位。为了解决这个问题，接下来将比较TEP的振幅，检查TMS的响应频率，并查看全局平均场功率。
锁时平均
当前数据集中有两个条件要比较：' relax ' & ' contract '。条件在数据集中通过EEG中的标记显示出来。读取数据时，基于这些标记的试次，可以在数据结构trialinfo字段中找到每个试次属于哪个条件的信息。在这里，“放松”条件用数字1表示，“收缩”条件用数字3表示。
>> data_tms_clean.trialinfo
ans =
1
1
1
1
.
.
.
3
3
3
.
.
.
我们可以使用该字段中的信息分别对这两个条件执行锁时分析。trialinfo字段中的每一行都对应于数据集中的一个试次。
在计算锁时平均值时，应用基线校正（TMS发生前50ms到1ms），进行35Hz的低通滤波。
cfg = [];
cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = [-0.05 -0.001];cfg.preproc.lpfilter = 'yes';
cfg.preproc.lpfreq = 35;
% Find all trials corresponding to the relax conditioncfg.trials = find(data_tms_clean.trialinfo==1);relax_avg = ft_timelockanalysis(cfg, data_tms_clean);
% Find all trials corresponding to the contract conditioncfg.trials = find(data_tms_clean.trialinfo==3);
contract_avg = ft_timelockanalysis(cfg, data_tms_clean);
计算两个条件的平均值之间的差，使用函数ft_math，对一个或多个数据结构执行一定数量的数学操作。
% calculate difference
cfg = [];
cfg.operation ='subtract'; % Operation to apply
cfg.parameter ='avg'; % The fieldinthe data structure towhichto apply the operation
difference_avg = ft_math(cfg, contract_avg, relax_avg);
使用ft_singleplotER绘制条件及其差异，该函数非常适合绘制和绘制比较条件。
% Plot TEPsofboth conditions
cfg = [];
cfg.layout ='easycapM10'; % Specifyingthisallows you to produce topographical plotsofyour datacfg.channel ='17';
cfg.xlim = [-0.10.6];
ft_singleplotER(cfg, relax_avg, contract_avg, difference_avg);
ylabel('Amplitude (uV)');
xlabel('time (s)');
title('Relax vs Contract');
legend({'relax''contract''contract-relax'});
ft_singleplotER一个很好的特点是，如果指定一个布局，你可以在绘图窗口中选择一个时间范围，然后单击它来生成该时间点的振幅会在地形图上显示出来。你也可以使用ft_topoplotER函数来完成此步操作。
%% Plotting topographies
figure;
cfg = [];cfg.layout = 'easycapM10';
cfg.xlim = 0:0.05:0.55; % Here we've specified a vector between 0 and 0.55 seconds in steps of 0.05 seconds. A topoplot will be created for each time point specified here.
cfg.zlim = [-2 2]; % Here you can specify the limit of thevaluescorrespondingtothe colors.Ifyoudonotspecify this the limits will be estimated automaticallyforeachplot making it difficulttocompare subsequent plots.
ft_topoplotER(cfg, difference_avg);
全局平均场功率
全局平均场功率（GMFP）是Lehmann和Skandries(1979)首次提出的一种测量方法，例如，Esser等人(2006)使用该方法来表征全局EEG活动。GMFP可以用以下公式计算（来自Esser et al. (2006)）。
其中t为时间，V为channel i处的电压，K为channel数。在Esser等人(2006)中，GMFP是根据所有被试的平均水平计算的。因为这里只有一个被试的数据，所以只计算这个被试的GMFP。但如果有多个被试，那么可以在进行总平均时用相同的方法进行计算。基本上，GMFP是channel上的标准差。
FieldTrip有一个内置的函数来计算GMFP；ft_globalmeanfield，这个函数需要输入锁时数据。这里将使用与Esser等人(2006)类似的预处理。
% Create time-locked average
cfg = [];cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = [-0.1 -.001];cfg.preproc.bpfilter = 'yes';
cfg.preproc.bpfreq = [5 100];cfg.trials = find(data_tms_clean.trialinfo==1);
% 'relax' trialsrelax_avg = ft_timelockanalysis(cfg, data_tms_clean);
cfg.trials = find(data_tms_clean.trialinfo==3);
% 'contract' trialscontract_avg = ft_timelockanalysis(cfg, data_tms_clean);
% GMFP calculation
cfg = [];
cfg.method = 'amplitude';
relax_gmfp = ft_globalmeanfield(cfg, relax_avg);
contract_gmfp = ft_globalmeanfield(cfg, contract_avg);
现在可以画出两个条件的GMFP。
%Plot GMFP
figure;
plot(relax_gmfp.time, relax_gmfp.avg,'b');holdon;
plot(contract_gmfp.time, contract_gmfp.avg,'r');
xlabel('time (s)');
ylabel('GMFP (uv^2)');
legend({'Relax''Contract'});
xlim([-0.10.6]);
ylim([03]);
时频分析
把信号分解成频率，然后看这些频率的功率平均值。首先使用ft_freqanalysis将信号分解成不同的频率。在做光谱分析时，重要的是在分解成不同频率之前去趋势和去均值，以避免功率谱看起来很奇怪。因此，可以使用preproc选项对数据进行去趋势和去均值。
% Calculate Induced TFRs fpor both conditions
cfg = [];
cfg.polyremoval =1; % Removes mean and linear trend
cfg.output ='pow'; % Output the powerspectrum
cfg.method ='mtmconvol';cfg.taper ='hanning';
cfg.foi =1:50; % Our frequenciesofinterest. Now:1to50,instepsof1.
cfg.t_ftimwin =0.3.*ones(1,numel(cfg.foi));
cfg.toi =-0.5:0.05:1.5;
% Calculate TFRforrelax trials
cfg.trials = find(data_tms_clean.trialinfo==1);
relax_freq = ft_freqanalysis(cfg, data_tms_clean);
% Calculate TFRforcontracttrials
cfg.trials = find(data_tms_clean.trialinfo==3);
contract_freq = ft_freqanalysis(cfg, data_tms_clean);
计算条件之间的差异。通常在绘制TFRs时，指定一个基线窗口。由于这里需计算条件之间的差异，并且对基线校正后的两个条件之间的差异感兴趣，所以这里需要先从条件中删除基线。
% Remove baseline
cfg = [];
cfg.baselinetype = 'relchange'; % Calculate the change relative to the baseline ((data-baseline) / baseline). You can also use 'absolute', 'relative', or 'db'.
cfg.baseline = [-0.5 -0.3];
relax_freq_bc = ft_freqbaseline(cfg, relax_freq);
contract_freq_bc = ft_freqbaseline(cfg, contract_freq);
% Calculate the difference between both conditions
cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'powspctrm';
difference_freq = ft_math(cfg, contract_freq_bc, relax_freq_bc);
现在已经计算了两种条件的TFR及其差异，我们可以用不同的方法绘制结果。使用ft_multiplotTFR脚本以头部2D形式绘制所有TFR。
cfg = [];
cfg.xlim = [-0.1 1.0];
cfg.zlim = [-1.5 1.5];
cfg.layout = 'easycapM10';
figure;
ft_multiplotTFR(cfg, difference_freq);
此图是完全交互式的，单击并拖动以选择一个或多个channel，单击以查看所选channel的均值。还可以使用ft_singleplotTFR在单个视图中绘制单个（或多个）channel。
cfg = [];
cfg.channel ='17'; % Specify the channel to plot
cfg.xlim = [-0.11.0]; % Specify thetimerange to plot
cfg.zlim = [-33];
cfg.layout ='easycapM10';
figure;
subplot(1,3,1); % Use MATLAB's subplot function to divide plots into one figure
ft_singleplotTFR(cfg, relax_freq_bc);
ylabel('Frequency (Hz)');
xlabel('time(s)');
title('Relax');
subplot(1,3,2);
ft_singleplotTFR(cfg, contract_freq_bc);
title('Contract');
ylabel('Frequency (Hz)');
xlabel('time(s)');
subplot(1,3,3);
cfg.zlim = [-1.5 1.5];
ft_singleplotTFR(cfg, difference_freq);
title('Contract - Relax');
ylabel('Frequency (Hz)');
xlabel('time(s)');
