静息态 EEG 微状态分析（上）

% ---- 第 1 步：GFP 曲线与峰值提取 ----
% 前提：数据已完成模块 2 的预处理，并使用平均参考（微状态对参考极其敏感）
fs   = EEG.srate;
data = EEG.data - mean(EEG.data, 1);   % 再做一次空间中心化，保证拓扑"相对平均"
gfp  = sqrt(mean(data.^2, 1));         % GFP(t)：每个时刻的空间标准差

figure; plot(EEG.times/1000, gfp);     % 目检：GFP 应是快速起伏的曲线
xlabel('时间 (s)'); ylabel('GFP (\muV)'); title('全局场功率');

% 只保留 GFP 局部极大值处的拓扑：间隔至少 10ms、高度超过均值
[pk, idx] = findpeaks(gfp, 'MinPeakDistance', round(0.01*fs), ...
                      'MinPeakHeight', mean(gfp));
maps = data(:, idx);                   % 通道 x 峰值时刻：候选拓扑集合
fprintf('采样拓扑数 %d（占全部时间点 %.1f%%）\n', ...
        size(maps,2), 100*size(maps,2)/size(data,2));
% findpeaks 来自 Signal Processing Toolbox（后面的 filtfilt 也是）

% ---- 第 2 步：极性无关 k-means 聚类（教学简化版）----
K = 4;  nInit = 50;                     % 经典微状态研究常取 K=4
nS = size(maps, 2);
mapsN = maps ./ vecnorm(maps);          % 样本归一化：只比较"形状"，不比较强弱
gfp2  = sum(maps.^2, 1);                % 每个样本的 GFP^2（后面当权重）
bestGEV = -inf;
for init = 1:nInit                      % 随机重启：对抗 k-means 的初值敏感性
    C = maps(:, randperm(nS, K));       % 初始中心：随机挑 K 张真实拓扑
    for it = 1:50                       % 教学版固定迭代；工程版加收敛判断
        sim = abs((C ./ vecnorm(C))' * mapsN);   % 分配步：K x nS 的|相关|矩阵
        [~, seg] = max(sim, [], 1);
        for k = 1:K                     % 更新步：类内空间第一主成分当新中心
            Xk = mapsN(:, seg == k);
            if size(Xk, 2) < 2, continue; end
            [u, ~, ~] = svd(Xk * Xk', 'econ');
            C(:, k) = u(:, 1);
        end
    end
    % 用全局解释方差 GEV 给本次初始化打分（这 K 张地图解释了多少电位方差）
    gev = 0;
    for k = 1:K
        sel = (seg == k);
        cn = C(:, k) / norm(C(:, k));
        gev = gev + sum((cn' * mapsN(:, sel)).^2 .* gfp2(sel));
    end
    gev = gev / sum(gfp2);
    if gev > bestGEV
        bestGEV = gev;
        Cbest = C ./ vecnorm(C);        % 归一化的 K 张原型地图
    end
end
C = Cbest;
fprintf('K = %d, 全局解释方差 GEV = %.1f%%\n', K, 100*bestGEV);
% 直觉：GEV 是"这 K 张地图演了全剧的多大比例"；K=4 时静息态典型值约 70-80%

% ---- 第 3 步：反向拟合——给每个时间点贴标签 ----
dataN = data ./ max(vecnorm(data), eps);   % 逐点归一化：不让高 GFP 时刻主导
sim   = abs(C' * dataN);                   % K x 时间点 的|空间相关|
[~, labels] = max(sim, [], 1);             % 每个时刻最像哪张地图
labels(gfp < 0.5*mean(gfp)) = 0;           % GFP 太低的时刻不归类（0 = 噪声主导）

% 时间平滑 + 最短持续时间：真实微状态持续约 60ms 以上
labels = medfilt1(double(labels), 3);      % 3 点中值滤波去掉孤立标签（奇数窗不出小数）
% labels = enforce_min_duration(labels, round(0.03*fs));  % <30ms 并入邻段：课后练习 3

% ---- 第 4 步：指标计算（游程编码思路：先切段，再按类汇总）----
d      = diff(labels);
edges  = [1, find(d ~= 0) + 1, numel(labels) + 1];   % 每一段的起止边界
durPts = diff(edges);                                % 每段持续的采样点数
segLab = labels(edges(1:end-1));                     % 每段的微状态编号
valid  = segLab > 0;                                 % 排除未归类(0)的段
recSec = numel(labels) / fs;

occurrence = accumarray(segLab(valid)', 1, [K 1]) / recSec;              % 次/秒
nSeg       = accumarray(segLab(valid)', 1, [K 1]);
coverage   = accumarray(segLab(valid)', durPts(valid)') / sum(durPts);   % 时间占比
duration   = accumarray(segLab(valid)', durPts(valid)', [K 1]) ./ max(nSeg, 1);
duration_ms = duration / fs * 1000;                  % 每段平均持续毫秒数

% 转移概率矩阵：行 = 当前状态，列 = 下一状态（行归一化的条件概率）
T = zeros(K);
seq = segLab(valid);
for s = 1:numel(seq) - 1
    if seq(s) ~= seq(s+1)
        T(seq(s), seq(s+1)) = T(seq(s), seq(s+1)) + 1;
    end
end
T = T ./ max(sum(T, 2), 1);
disp('转移概率矩阵：'); disp(round(T, 3));