使用Matlab进行基于时间序列数据的滑窗分析解决方案

心理学的实验常常会记录到时间序列数据，包括脑电数据、眼动数据以及近红外成像的血氧信号等。由于仪器设备的采样率较高，在描述时间序列时，常常选定一个固定的时间窗口，以特定的步长滑动窗口，将每个窗口中的数据进行平均得到一组新的数据。这样得到的新的数据数据量较少，也更平滑。

在实际的数据处理过程中碰到了这样的一个问题：在2×2的混合设计中，使用EEG记录了被试的脑电，并已经按照ERP的数据处理流程，根据组内变量的不同实验条件，将每个被试、每个电极点的数据分段并叠加平均好了。现在想看不同实验条件（组内变量）下的差异波在不同组别间（组间变量）的差异是如何变化的（类似于交互作用）。这是一个典型的，需要使用Matlab进行基于时间序列数据的滑窗分析的应用场景，以下是代码和解决过程。

先简单介绍实验背景：被试被分到两组，分别使用不同的学习方法学习。然后使用经典的新旧词识别范式检测被试的学习效果，同时记录被试的脑电。以刺激出现的时刻为0点，以[-200~1000]的区间进行分段并叠加平均，可以得到一个被试判断新词和判断旧词时的脑电反应，使用新词减旧词得到的差异波（这里也可以是旧词减新词）可以看作被试脑电上的新旧效应。现在想看不同组别间新旧效应的差异是如何变化的。

首先，第一步是将所有被试的数据导入Matlab。

这里我拿到的数据是从Analyzer软件中导出的已经分段、叠加、平均结束的，不同电极点的一般数据，直接用死办法，将代码文件和所有被试数据文件放在同一文件夹，使用importdata命令将原始数据全部导入。

data = importdata('data.name')

第二步，将原始数据中相同兴趣区的不同电极点进行平均。

这里被试间变量为组别，一共两个组，分别为Ping组和Zhiz组。每个组共24名被试，每名被试的数据根据电极点分为4个兴趣区，分别为左侧额叶lf、右侧额叶rf、左侧顶叶lc和右侧顶叶rc。使用for循环对每个被试的数据进行批处理。结果得到每个被试在四个兴趣区各有一行数据，这行数据是该被试在该兴趣区脑电振幅的时间序列数据。

%% 平均电极点
lf_Ping_Sum = zeros(24, 600);
rf_Ping_Sum = zeros(24, 600);
lc_Ping_Sum = zeros(24, 600);
rc_Ping_Sum = zeros(24, 600);
lf_Zhiz_Sum = zeros(24, 600);
rf_Zhiz_Sum = zeros(24, 600);
lc_Zhiz_Sum = zeros(24, 600);
rc_Zhiz_Sum = zeros(24, 600);
for i =1:24
    %组1
    Ping = eval(['Ping', num2str(i)]); %将字符串返回为变量名
    lf_Ping = [Ping.data(1, :); Ping.data(3, :); Ping.data(11, :); Ping.data(17, :); Ping.data(21, :); Ping.data(23, :)];
    rf_Ping = [Ping.data(2, :); Ping.data(4, :); Ping.data(12, :); Ping.data(18, :); Ping.data(22, :); Ping.data(24, :)];
    lc_Ping = [Ping.data(5, :); Ping.data(7, :); Ping.data(9, :); Ping.data(13, :); Ping.data(15, :); Ping.data(19, :)];
    rc_Ping = [Ping.data(6, :); Ping.data(8, :); Ping.data(10, :); Ping.data(14, :); Ping.data(16, :); Ping.data(20, :)];
    lf_Ping_Sum(i, :) = mean(lf_Ping);
    rf_Ping_Sum(i, :) = mean(rf_Ping);
    lc_Ping_Sum(i, :) = mean(lc_Ping);
    rc_Ping_Sum(i, :) = mean(rc_Ping);
    
    %组2
    Zhiz = eval(['Zhiz', num2str(i)]);
    lf_Zhiz = [Zhiz.data(1, :); Zhiz.data(3, :); Zhiz.data(11, :); Zhiz.data(17, :); Zhiz.data(21, :); Zhiz.data(23, :)];
    rf_Zhiz = [Zhiz.data(2, :); Zhiz.data(4, :); Zhiz.data(12, :); Zhiz.data(18, :); Zhiz.data(22, :); Zhiz.data(24, :)];
    lc_Zhiz = [Zhiz.data(5, :); Zhiz.data(7, :); Zhiz.data(9, :); Zhiz.data(13, :); Zhiz.data(15, :); Zhiz.data(19, :)];
    rc_Zhiz = [Zhiz.data(6, :); Zhiz.data(8, :); Zhiz.data(10, :); Zhiz.data(14, :); Zhiz.data(16, :); Zhiz.data(20, :)];
    lf_Zhiz_Sum(i, :) = mean(lf_Zhiz);
    rf_Zhiz_Sum(i, :) = mean(rf_Zhiz);
    lc_Zhiz_Sum(i, :) = mean(lc_Zhiz);
    rc_Zhiz_Sum(i, :) = mean(rc_Zhiz);
end

第三步，进行滑窗的拆分。

这里一定要注意，最好将基本参数空出来，这样后面修改代码时就会更简单，直接修改参数即可。一共有5个基本参数：数据起始时间点sample_start、数据结束时间点sample_end、设备采样率samplerate、滑窗时间长度slidewind_length和滑窗时间步长slidestep_length（变量名可以根据自己的喜好定义）。这一步的目的时确定每一个窗口的起始点和窗口的个数。

%% 滑窗拆分（参数可修改）
%采样率500Hz，完整数据1200ms（data points 600），时间窗50ms（data points 25），窗口步长10ms（data points 5）
sample_start = -200; %起始点 -200ms
sample_end = 1000; %结束点 1000ms
samplerate = 500; %单位 Hz
samplerate_t = 1000*1/samplerate; %单位 ms
epoch_length = sample_end - sample_start; %单位 ms
slidewind_length = 50; %单位 ms
slidestep_length = 10; %单位 ms
datapoints_total = epoch_length/samplerate_t;
datapoints_slidewind = slidewind_length/samplerate_t;
datapoints_slidestep = slidestep_length/samplerate_t;
sw_start = [1: datapoints_slidestep: (datapoints_total-datapoints_slidewind+1)]; %确定滑窗起始点
Nsw =  length(sw_start); %确定滑窗个数

第四步，进行滑窗的运算。

这里使用for循环语句，进行滑窗的运算，对每个窗口下的平均幅值进行独立样本t检验（注意：Matlab中独立样本t检验的函数是ttest2）

[h,p,ci,t] = ttest2(X,Y);
t_value = t.tstat

这里只展示左侧额叶的处理代码，其他的兴趣区只需重复执行此代码，改下变量名称即可。这里也可以用for循环语句，但是只有四个兴趣区，我就没有写循环结构。

%% 独立样本t检验
%左侧额叶
lf_P = zeros(1, Nsw); lf_T = zeros(1, Nsw);
for i = 1:Nsw
    lf_ping = mean(lf_Ping_Sum(:, sw_start(i): sw_start(i)+datapoints_slidewind-1), 2); %向量的引用以及基于列向量的平均
    lf_zhi = mean(lf_Zhiz_Sum(:, sw_start(i): sw_start(i)+datapoints_slidewind-1), 2);
    [lf_a, lf_b, lf_c, lf_d] = ttest2(lf_ping, lf_zhi);
    lf_P(i) = lf_b; lf_T(i) = lf_d.tstat;
end
% lf_P = mafdr(lf_P, 'BHFDR', true); %多重比较校正

最后，对结果进行可视化处理。

这里主要使用Matlab中的plot函数对结果进行可视化处理。注意不同画布中的横坐标尺度（时间）要保持一致才有可比性，同样以左侧额叶兴趣区为例。这里plot函数的使用参考了网友的总结，非常全面：【MATLAB】基本绘图函数（涵盖所有基本绘图指令）_Do My Love的博客-CSDN博客_matlab函数绘图

%% 可视化
%左侧额叶
figure; %建立新图窗
subplot(3, 1, 1); %将图窗分割成三行一列的三个子图窗
plot([sample_start: samplerate_t: sample_end - samplerate_t], mean(lf_Ping_Sum), 'b'); %绘制平板组原始值
hold on; %在同一张画布上继续绘图
plot([sample_start: samplerate_t: sample_end - samplerate_t], mean(lf_Zhiz_Sum), 'r'); %绘制纸质组原始值
hold off;
axis([sample_start, sample_end, -3, 1]); % 设置坐标轴尺度
legend({'平板组', '纸质组'},'Location', 'southwest') %设置图例
xlabel('time') %设置x轴名称
ylabel('new-old waves') %设置y轴名称
title('左侧额叶结果') %设置图的名称
subplot(3, 1, 2); 
plot([sample_start: slidestep_length: sample_end - slidewind_length], lf_T, 'k');  %绘制t值
axis([sample_start, sample_end, -2.5, 2.5]);
xlabel('time')
ylabel('t-value')
subplot(3, 1, 3); 
plot([sample_start: slidestep_length: sample_end - slidewind_length], lf_P, 'k');  %绘制p值
axis([sample_start, sample_end, 0, 1]);
xlabel('time')
ylabel('p-value')

最终能够看到两个组的被试的新旧效应（差异波）的差异随着时间的动态变化，作为寻找规律与科学解释的参考。

欢迎交流与分享，可从联系页面联系到我。