超详细TMS-EEG数据处理教程（下）

文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。

上一期的文章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 differencecfg = [];cfg.operation = 'subtract'; % Operation to applycfg.parameter = 'avg'; % The field in the data structure to which to apply the operationdifference_avg = ft_math(cfg, contract_avg, relax_avg);

使用ft_singleplotER绘制条件及其差异，该函数非常适合绘制和绘制比较条件。

% Plot TEPs of both conditionscfg = [];cfg.layout = 'easycapM10'; % Specifying this allows you to produce topographical plots of your datacfg.channel = '17';cfg.xlim = [-0.1 0.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 topographiesfigure;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 the values corresponding to the colors. If you do not specify this the limits will be estimated automatically for each plot making it difficult to compare 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 averagecfg = [];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 calculationcfg = [];cfg.method = 'amplitude';relax_gmfp = ft_globalmeanfield(cfg, relax_avg);contract_gmfp = ft_globalmeanfield(cfg, contract_avg);

现在可以画出两个条件的GMFP。

%Plot GMFPfigure;plot(relax_gmfp.time, relax_gmfp.avg,'b');hold on;plot(contract_gmfp.time, contract_gmfp.avg,'r');xlabel('time (s)');ylabel('GMFP (uv^2)');legend({'Relax' 'Contract'});xlim([-0.1 0.6]);ylim([0 3]);

时频分析

把信号分解成频率，然后看这些频率的功率平均值。首先使用ft_freqanalysis将信号分解成不同的频率。在做光谱分析时，重要的是在分解成不同频率之前去趋势和去均值，以避免功率谱看起来很奇怪。因此，可以使用preproc选项对数据进行去趋势和去均值。

% Calculate Induced TFRs fpor both conditionscfg = [];cfg.polyremoval     = 1; % Removes mean and linear trendcfg.output          = 'pow'; % Output the powerspectrumcfg.method          = 'mtmconvol';cfg.taper           = 'hanning';cfg.foi             = 1:50; % Our frequencies of interest. Now: 1 to 50, in steps of 1.cfg.t_ftimwin       = 0.3.*ones(1,numel(cfg.foi));cfg.toi             = -0.5:0.05:1.5;
% Calculate TFR for relax trialscfg.trials         = find(data_tms_clean.trialinfo==1);relax_freq         = ft_freqanalysis(cfg, data_tms_clean);
% Calculate TFR for contract trialscfg.trials         = find(data_tms_clean.trialinfo==3);contract_freq      = ft_freqanalysis(cfg, data_tms_clean);

计算条件之间的差异。通常在绘制TFRs时，指定一个基线窗口。由于这里需计算条件之间的差异，并且对基线校正后的两个条件之间的差异感兴趣，所以这里需要先从条件中删除基线。

% Remove baselinecfg = [];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 conditionscfg = [];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 plotcfg.xlim = [-0.1 1.0]; % Specify the time range to plotcfg.zlim = [-3 3];cfg.layout = 'easycapM10';
figure;subplot(1,3,1); % Use MATLAB's subplot function to divide plots into one figureft_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)');