基于 FieldTrip 的功能连接实践

%% ===== 头皮电极水平连接分析（单被试完整流程） =====
clear; clc; close all;
ft_defaults                               % 启动 FieldTrip（每个脚本开头跑一次）

subj = 'sub01';

% --- 1) 读入模块 2 预处理好的 EEGLAB .set（FieldTrip 能直接读 .set） ---
cfg = [];
cfg.dataset = ['D:\project\processed\' subj '_pp.set'];
cfg.channel = 'eeg';                      % 只保留脑电通道（自动剔除 EOG/ECG 等）
data = ft_preprocessing(cfg);

% （若数据是未分段的连续记录：先切等长小段充当"试次"，连接统计需要重复观测）
% cfg = []; cfg.length = 2;               % 每段 2 秒
% data = ft_redefinetrial(cfg, data);

% --- 2) 频域分解：要 Fourier 谱，不要功率谱 ---
% 连接靠相位；功率谱是"相位盲"的：output 必须设为 'fourier'
cfg = [];
cfg.method     = 'mtmfft';                % 多窗 FFT：逐试次的频谱估计
cfg.output     = 'fourier';               % 关键：输出复数频谱（保留相位）
cfg.taper      = 'dpss';
cfg.foilim     = [8 13];                  % 只算 alpha 频段：内存与计算量直线下降
cfg.tapsmofrq  = 2;                       % 2Hz 频率平滑：稳定相位估计
cfg.keeptrials = 'yes';                   % 保留试次维！互谱要跨试次求
freq = ft_freqanalysis(cfg, data);

% --- 3) 连接计算 ---
cfg = [];
cfg.method     = 'wpli';                  % 头皮水平抗零滞后首选（6.1 对比表）
cfg.channelcmb = {'F3' 'O1'; 'F4' 'O2'; 'Fz' 'Pz'};   % 只算这 3 对：内存闸门
conn = ft_connectivityanalysis(cfg, freq);
disp(conn)                                % 重点看两个字段：labelcmb 和 wplispctrm

% --- 4) 自带绘图：通道对 x 频率的折线图 ---
cfg = [];
cfg.parameter = 'wplispctrm';             % 必须与输出字段名完全一致
cfg.zlim      = [0 1];
ft_connectivityplot(cfg, conn);
% 注意：wpli 的输出可能带符号（含方向信息），看"强度"时对结果取 abs

function C = conn2matrix(conn, param, freqOfInterest)
% 把 ft_connectivityanalysis 的输出（labelcmb + *_spctrm）转成 N×N 对称矩阵
% 转完直接喂给 6.2 的矩阵图 / draw_circle_net 代码
% conn: 连接输出; param: 'cohspctrm' / 'wplispctrm' / 'plvspctrm'
% freqOfInterest: 感兴趣频率（Hz），取最近的频率 bin；不填则对整个频段取均值
labels = unique(conn.labelcmb(:));        % 由通道对反推通道名（按字母排序）
N = numel(labels);
vals = squeeze(conn.(param));             % squeeze 后：通道对 x 频率
if nargin < 3 || isempty(freqOfInterest)
    v = mean(vals, 2);                    % 不指定频率：整个频段取均值
else
    [~, kfr] = min(abs(conn.freq - freqOfInterest));
    v = vals(:, kfr);
end
C = nan(N, N);
for k = 1:size(conn.labelcmb, 1)
    i = find(strcmp(labels, conn.labelcmb{k,1}));
    j = find(strcmp(labels, conn.labelcmb{k,2}));
    C(i, j) = v(k);
    C(j, i) = v(k);                       % 对称指标：两格都填
end
end

% 用法：和 6.2 的画图代码无缝衔接
% C = conn2matrix(conn, 'wplispctrm', 10);   % 取最接近 10Hz 的结果
% draw_circle_net(C, labels, 0.3);           % 6.2 的圆形网络图

%% ===== 源空间连接：先源后连（骨架） =====
% 前置（模块 7）：headmodel = ft_prepare_headmodel(...)；
%                sourcemodel = ft_prepare_leadfield(...)（含每个源点的导联场）

% 1) 在数据上估计 LCMV 空域滤波器，并保留滤波器
cfg = [];
cfg.method = 'lcmv';
cfg.headmodel = headmodel;
cfg.sourcemodel = sourcemodel;
cfg.lcmv.keepfilter = 'yes';              % 保留每个源点的滤波器（下一步要用）
source = ft_sourceanalysis(cfg, dataAll);
% 注意：raw 输入时滤波器在 source.filter，timelock 输入时在 source.avg.filter

% 2) 用同一组滤波器提取各 ROI 的虚拟电极时间序列
cfg = [];
cfg.method = 'lcmv';
cfg.sourcemodel = sourcemodel;
cfg.sourcemodel.filter = source.avg.filter;   % 复用第 1 步估计的滤波器
roi = ft_sourceanalysis(cfg, dataTrial);      % roi.time{trial} = 各源点的时间序列

% 3) 之后与头皮水平完全一样：
%    对 ROI 时间序列做 ft_freqanalysis(output='fourier') -> ft_connectivityanalysis
%    区别只是：labelcmb 里现在是"ROI 对"而不是"电极对"

%% ===== 组级骨架：循环被试 -> 被试 x 边矩阵 -> 6.2 的统计 =====
subjects = {'sub01','sub02','sub03','sub04','sub05'};
edges    = {'F3','O1'; 'F4','O2'; 'Fz','Pz'; 'C3','C4'; 'T7','T8'};  % 预选的边
nSub  = numel(subjects);
nEdge = size(edges, 1);
W = nan(nSub, nEdge);                     % 被试 x 边 的 wPLI 矩阵（核心结果）

for s = 1:nSub
    try
        freq = run_subject_connectivity(subjects{s});  % 第三节第 1-2 步打包的函数
        cfg = [];
        cfg.method = 'wpli';
        cfg.channelcmb = edges;
        conn = ft_connectivityanalysis(cfg, freq);
        v = squeeze(conn.wplispctrm);
        if size(v, 2) > 1, v = mean(v, 2); end    % 频段内取均值
        W(s, :) = abs(v)';                        % abs 取强度（wpli 输出可带符号）
        fprintf('%s done\n', subjects{s});
    catch err
        fprintf(2, '%s FAILED: %s\n', subjects{s}, err.message);  % 单个被试失败不中断
        W(s, :) = NaN;                            % 标记缺失，统计时剔除
    end
end

% --- 两条件配对比较（示意）：逐边配对 t 检验 + 6.2 的 BH-FDR ---
% WA、WB：两条件各自的"被试 x 边"矩阵（同上循环各跑一遍）
p = nan(nEdge, 1);
for e = 1:nEdge
    [~, p(e)] = ttest(WA(:, e), WB(:, e));  % ttest 需 Statistics Toolbox；也可换 signflip 置换
end
sig = fdr_bh(p, 0.05);                     % 6.2 写的 BH-FDR，直接复用
fprintf('通过 FDR 的边：%d / %d\n', nnz(sig), nEdge);

save('group_wpli.mat', 'subjects', 'edges', 'W', 'p', 'sig');  % 存档：可复现