超详细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)');

超详细TMS-EEG数据处理教程（下）相关推荐

超详细的gnuplot使用教程【2】
超详细的gnuplot使用教程 1.gnuplot参数介绍及演示 1.1首先来解释一下会用到的各类参数以及其解释 1.2 画图实际测试 1.3 其它参数介绍(约定范围.坐标轴设定) 1.3.1 约束画 ...
超详细的 Wireshark 使用教程
原文超详细的 Wireshark 使用教程一.wireshark是什么? wireshark是非常流行的网络封包分析软件,简称小鲨鱼,功能十分强大.可以截取各种网络封包,显示网络封包的详细信息. ...
【02】2022.11最新超详细Vuforia图片识别教程
[02]2022.11最新超详细Vuforia图片识别教程文章目录 [02]2022.11最新超详细Vuforia图片识别教程 1.Vuforia环境搭建 2.License Key获取及注册 3. ...
Anaconda超详细下载安装配置教程(Windows)
Anaconda最新超详细下载安装配置教程(Windows) 命令总结写在最前面 1.查看conda版本: conda --version 2.进入python交互模式: python 3.退出pyt ...
超详细的MySQL入门教程（四）
MySQL:简单的增删改查查询数据基本语法介绍打印任意值查询表中全部数据查询表中部分字段限定条件查询例1:查询编号值小于指定值的记录例2:查询地址不等于某值的记录例3:查询一级地址等 ...
【数据的存储】浮点数在内存中的存储详解【超详细的保姆级别教程，让面试官心服口服】手撕浮点数存储使用方式
[数据的存储]浮点数在内存中的存储详解[超详细的保姆级别教程,让面试官对你心服口服]手撕浮点数存储使用方式作者: @小小Programmer 这是我的主页:@小小Programmer 在食用这篇博客 ...
华为C8500S 超详细线刷刷机教程
华为C8500S 超详细线刷刷机教程,菜鸟也可以刷机手机, 华为C8500S, 教程华为C8500S详细线刷刷机教程,开始: 鉴于很多网友是新手,都在询问刷机过程,怎样刷机,下面我就做个详细教程, ...
超详细的cmake入门教程【转载】
这篇文章主要介绍了超详细的cmake入门教程,需要的朋友可以参考下源出处超详细的cmake入门教程什么是cmake 在 linux 平台下使用 CMake 生成 Makefile 并编译的流程 ...
【排序】什么都能排的C语言qsort排序详解【超详细的宝藏级别教程】深度理解qsort排序
[排序]什么都能排的C语言qsort排序详解[超详细的宝藏级别教程]深度理解qsort排序作者: @小小Programmer 这是我的主页:@小小Programmer 在食用这篇博客之前,博主在这里 ...
DBeaver安装与使用教程（超详细安装与使用教程）
文章预览: DBeaver安装与使用教程(超详细安装与使用教程) 一.DBeaver安装教程 ①下载地址 ②图文安装教程二.DBeaver使用教程 ①mysql数据库为例 1>填写数据库信息 ...

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

超详细TMS-EEG数据处理教程（下）相关推荐

最新文章

热门文章