FieldTrip 安装与时域 / 频域 / 时频分析

% 假设 data 是 ft_preprocessing 的输出，先整体看一眼
data                 % 直接打印结构体
data.trial           % 元胞数组：1 个元胞 = 1 个试次
data.trial{1}        % 第 1 个试次的矩阵：通道 x 时间点，例如 64 x 1000
data.label           % Nx1 元胞数组：通道名，如 Fp1、Cz、O2
data.time            % 1x1 元胞数组：时间轴（单位是秒），对齐后一般从 -0.5 开始
data.fsample         % 采样率，例如 500 (Hz)
data.trialinfo       % 每个试次一行的附加信息：条件编号、反应时等

% 三个自检动作（拿到任何新数据都先做）
size(data.trial{1})          % 通道数 x 每试次采样点数，与采集记录核对
numel(data.trial)            % 试次总数
unique(data.trialinfo(:,1))  % 出现过哪些条件编号，例如 [1; 2]

%% 第 1 步：定义试次（只算试次表，还没读数据）
cfg = [];
cfg.dataset  = 'D:\project\sub01.vhdr';   % 原始数据（英文路径！）
cfg.trialfun = 'ft_trialfun_general';     % 通用试次定义函数
cfg.trialdef.prestim    = 0.5;            % 事件前 0.5 秒
cfg.trialdef.poststim   = 1.5;            % 事件后 1.5 秒
cfg.trialdef.eventtype  = 'Stimulus';     % 事件类型名（取决于采集软件）
cfg.trialdef.eventvalue = [1 2];          % 触发值：1=条件A，2=条件B
cfg = ft_definetrial(cfg);               % 注意：输出还是 cfg，trl 已填好

cfg.trl(1:3, :)    % 每行三列：[起始采样点  结束采样点  事件偏移]

%% 第 2 步：读入数据并完成预处理（继续往同一个 cfg 上加）
cfg.channel   = {'EEG'};                 % 只读脑电通道，自动排除 EOG/ECG
cfg.hpfilter  = 'yes';  cfg.hpfreq  = 0.1;    % 高通 0.1 Hz：去慢漂移（理由同 2.3）
cfg.lpfilter  = 'yes';  cfg.lpfreq  = 45;    % 低通 45 Hz：挡肌电、留 ERP 主能量
cfg.dftfilter = 'yes';  cfg.dftfreq  = 50;   % 50 Hz 工频陷波
cfg.reref     = 'yes';  cfg.refchannel = {'all'};  % all = 平均参考
cfg.demean    = 'yes';  cfg.baselinewindow = [-0.5 0];  % 基线校正（单位秒）
cfg.padding   = 1;                       % 滤波时两侧各补 1 秒，压制边缘振铃

data = ft_preprocessing(cfg);            % 一口气：读数据 + 上面全部处理

% 面孔一：数值阈值（zvalue 法）——把数据先 z 变换，超过 4 个标准差算伪迹
cfg = [];
cfg.artfctdef.type = {'zvalue'};          % 用哪些检测方法，可叠加多种
cfg.artfctdef.zvalue.channel = {'EEG'};   % 在哪些通道上检测
cfg.artfctdef.zvalue.cutoff  = 4;         % 超过 4 个标准差视为伪迹
data_clean = ft_rejectartifact(cfg, data); % 默认整试次剔除（reject='complete'）

% 好习惯：把拒绝率打出来，写进实验日志
fprintf('剔除 %d / %d 个试次\n', ...
    numel(data.trial) - numel(data_clean.trial), numel(data.trial));

% 面孔二：视觉拒绝——交互浏览、手工点选（小数据集与抽查用）
cfg = [];
cfg.method  = 'summary';    % 汇总指标视图；也可选 'trial' / 'channel'
cfg.channel = 'EEG';
data_clean = ft_rejectvisual(cfg, data);  % 弹出浏览窗口，勾选要删的试次/通道

% 1) 时域分析：按条件分别平均
cfg = [];
cfg.trials = find(data.trialinfo(:,1) == 1);   % 条件 A 的试次编号
tlk_A = ft_timelockanalysis(cfg, data);

cfg = [];
cfg.trials = find(data.trialinfo(:,1) == 2);   % 条件 B
tlk_B = ft_timelockanalysis(cfg, data);

% 输出结构：tlk_A.avg（通道x时间的平均）、tlk_A.var、tlk_A.time、tlk_A.label

% 2) ERP 绘图三件套
cfg = [];
cfg.channel = 'Cz';            % 指定一个通道
cfg.xlim    = [-0.5 1.5];
ft_singleplotER(cfg, tlk_A);   % 单通道 ERP 曲线：论文主图的高频操作

cfg = [];
cfg.layout = 'elec1005';       % 标准 10-05 系统电极布局（EEG 通用）
ft_multiplotER(cfg, tlk_A);     % 全通道阵列图：一眼看全局形态

cfg = [];
cfg.layout   = 'elec1005';
cfg.xlim     = [0.3 0.5];      % 显示 300-500 ms 的时间窗（如 P3 窗口）
cfg.zlim     = 'maxabs';
cfg.colorbar = 'yes';
ft_topoplotER(cfg, tlk_A);     % 头皮地形图

%% 面孔一：mtmfft 算功率谱（PSD）——整段试次一次变换
cfg = [];
cfg.method    = 'mtmfft';      % 多窗 FFT
cfg.output    = 'pow';         % 输出功率
cfg.taper     = 'dpss';        % 多窗（Slepian 序列），比单窗谱更平滑
cfg.foilim    = [1 45];        % 1-45 Hz
cfg.tapsmofrq = 2;             % 谱平滑带宽 2 Hz：平滑与分辨率的交换
cfg.trials    = find(data.trialinfo(:,1) == 1);
freqA = ft_freqanalysis(cfg, data);
% 输出：freqA.powspctrm（通道x频率）、freqA.freq

% 画功率谱：x 轴是频率（复用了 TFR 绘图函数）
cfg = [];
cfg.channel = 'Cz';
cfg.xlim    = [1 45];          % 这里 xlim 指频率范围
ft_singleplotTFR(cfg, freqA);

%% 面孔二：mtmconvol 滑动窗算时频（TFR）
cfg = [];
cfg.method     = 'mtmconvol';        % 多窗卷积时频
cfg.output     = 'pow';
cfg.foi        = 2:1:40;              % 感兴趣频率
cfg.t_ftimwin  = 5 ./ cfg.foi;        % 窗长随频率变：每个频率取 5 个周期
cfg.tapsmofrq = 0.4 .* cfg.foi;      % 平滑带宽随频率变
cfg.toi        = -0.5:0.05:1.5;      % 时间中心点，50 ms 一步
cfg.trials     = find(data.trialinfo(:,1) == 1);
tfrA = ft_freqanalysis(cfg, data);
% 输出：tfrA.powspctrm（通道x频率x时间）

% 画时频图：单通道的"热图"
cfg = [];
cfg.channel = 'Cz';
cfg.xlim    = [-0.5 1.5];      % 时间轴
cfg.ylim    = [2 40];          % 频率轴
cfg.zlim    = 'maxabs';
ft_singleplotTFR(cfg, tfrA);

% 画时频地形图：某个时间-频率窗的空间分布
cfg = [];
cfg.layout = 'elec1005';
cfg.xlim   = [0.3 0.8];        % 时间窗
cfg.ylim   = [8 13];           % 频率窗：alpha 频段
ft_topoplotTFR(cfg, tfrA);

%% ============ 参数区：只改这里 ============
subjects = {'sub01','sub02','sub03'};     % 被试列表
dataDir  = 'D:\project\raw\';              % 原始数据目录（英文路径！）
outDir   = 'D:\project\ft_result\';        % 输出目录
conds    = [1 2];                          % 两个条件的触发值

if ~exist(outDir, 'dir'); mkdir(outDir); end
tlkA = cell(1, numel(subjects));           % 每被试每条件一个 timelock 结构
tlkB = cell(1, numel(subjects));

%% ============ 循环区：所有被试同一流程 ============
for k = 1:numel(subjects)
    fprintf('==== [%d/%d] %s ====\n', k, numel(subjects), subjects{k});

    % --- 试次定义 + 预处理（每人重配一份新 cfg，防止残留字段捣乱）---
    cfg = [];
    cfg.dataset  = [dataDir subjects{k} '.vhdr'];
    cfg.trialfun = 'ft_trialfun_general';
    cfg.trialdef.prestim    = 0.5;
    cfg.trialdef.poststim   = 1.5;
    cfg.trialdef.eventvalue = conds;
    cfg = ft_definetrial(cfg);

    cfg.channel   = {'EEG'};
    cfg.hpfilter  = 'yes';  cfg.hpfreq = 0.1;
    cfg.lpfilter  = 'yes';  cfg.lpfreq = 45;
    cfg.reref     = 'yes';  cfg.refchannel = {'all'};
    cfg.demean    = 'yes';  cfg.baselinewindow = [-0.5 0];
    data = ft_preprocessing(cfg);

    % --- 伪迹拒绝（数值阈值，可复现）---
    cfg = [];
    cfg.artfctdef.type = {'zvalue'};
    cfg.artfctdef.zvalue.channel = {'EEG'};
    cfg.artfctdef.zvalue.cutoff  = 4;
    data = ft_rejectartifact(cfg, data);

    % --- 按条件分别算 ERP ---
    cfg = [];  cfg.trials = find(data.trialinfo(:,1) == conds(1));
    tlkA{k} = ft_timelockanalysis(cfg, data);
    cfg = [];  cfg.trials = find(data.trialinfo(:,1) == conds(2));
    tlkB{k} = ft_timelockanalysis(cfg, data);

    nA = sum(data.trialinfo(:,1) == conds(1));
    nB = sum(data.trialinfo(:,1) == conds(2));
    fprintf('    保留试次：A %d 个，B %d 个\n', nA, nB);
end

%% ============ 保存：5.2 统计的直接输入 ============
save(fullfile(outDir, 'timelock_all.mat'), 'subjects', 'tlkA', 'tlkB');