ERP 波形计算与地形图绘制

% 用模拟数据直观感受"平均"的威力
% 每个试次 = 固定形状的"ERP"(信号) + 随机噪声
fs = 500;                                    % 采样率 (Hz)
t  = -0.2 : 1/fs : 0.8;                      % -200 ~ 800 ms 时间轴
sig = 6 * exp(-((t - 0.3) / 0.05).^2);       % 模拟 300ms 的正波，峰值 6 uV
nTrial = 100;
data = zeros(length(t), nTrial);             % 时间点 x 试次
for k = 1:nTrial
    data(:, k) = sig + 8 * randn(length(t), 1);   % 噪声比信号还大——真实脑电的常态
end

figure;
subplot(3, 1, 1); plot(t, data(:, 1));        title('单个试次：信号被噪声淹没');
subplot(3, 1, 2); plot(t, mean(data(:, 1:10), 2)); title('平均 10 个试次：隐约成形');
subplot(3, 1, 3); plot(t, mean(data, 2));     title('平均 100 个试次：ERP 浮出水面');
xlabel('时间 (s)'); ylabel('幅度 (\muV)');

% 前提：EEG 是模块 2 结束时的分段数据（已基线校正）
ERP = mean(EEG.data, 3);                    % 沿第 3 维（试次维）平均 -> 通道 x 时间
[~, chIdx] = ismember({'Pz'}, {EEG.chanlocs.labels});   % 找到 Pz 的通道编号

figure; plot(EEG.times, ERP(chIdx, :), 'LineWidth', 1.2); hold on;
xline(0, '--');                             % 刺激 onset 参考线
xlabel('时间 (ms)'); ylabel('幅度 (\muV)'); title('单被试 Pz 的 ERP 波形');

% EEGLAB 自带的画 ERP 函数：GUI 中 Plot > Channel ERPs 的脚本等价形式
pop_erpplot(EEG, 1, [chIdx], 0);            % 通道列表；0 表示不叠加标准差

%% 组平均：先算每个人的 ERP，再跨被试平均
clear; clc; close all;
dataDir = 'D:\project\processed\';          % 预处理输出目录（英文路径）
subs    = {'sub01','sub02','sub03','sub04','sub05'};

erpAll = [];  times = [];
for s = 1:numel(subs)
    EEG = pop_loadset('filename', [subs{s} '_pp.set'], 'filepath', dataDir);
    erp = mean(EEG.data, 3);                % 这名被试的 ERP：通道 x 时间
    if isempty(erpAll)
        % 用第一名被试初始化尺寸与通道编号
        [~, chIdx] = ismember({'Pz'}, {EEG.chanlocs.labels});
        erpAll = zeros(numel(subs), numel(EEG.times));
        times  = EEG.times;                 % 预处理一致 -> 时间轴一致
    end
    erpAll(s, :) = erp(chIdx, :);           % 收集 Pz 波形：每行一个被试
    fprintf('%s：剩余 %d 个试次\n', subs{s}, EEG.trials);   % 顺手做质控
end
GA = mean(erpAll, 1);                       % 跨被试平均 -> 组平均

figure; hold on;
plot(times, erpAll, 'Color', [0.75 0.75 0.75]);   % 灰色：每个被试
plot(times, GA, 'r', 'LineWidth', 2);             % 红色：组平均
xlabel('时间 (ms)'); ylabel('幅度 (\muV)');
title('灰色为单被试，红色为 grand average');

%% 手工画 ERPimage：imagesc 就够了
X = squeeze(EEG.data(chIdx, :, :))';        % 原为 1 x 时间 x 试次，转置成 试次 x 时间

figure;
subplot(2, 1, 1);
imagesc(EEG.times, 1:EEG.trials, X);        % 每一行 = 一个试次的完整波形
set(gca, 'YDir', 'normal'); colorbar; axis tight;
xlabel('时间 (ms)'); ylabel('试次编号'); title('按采集顺序排列');

% 换一种排法：按 250-350ms 窗口内的平均幅值从小到大排序
subplot(2, 1, 2);
idxWin = EEG.times >= 250 & EEG.times <= 350;
[~, order] = sort(mean(X(:, idxWin), 2));   % 每个试次在窗口内的幅值
imagesc(EEG.times, 1:EEG.trials, X(order, :));
set(gca, 'YDir', 'normal'); colorbar; axis tight;
xlabel('时间 (ms)'); ylabel('试次（按幅值排序）'); title('同一份数据，按窗口幅值排序');

% 画 300ms 时刻的地形图：给每个电极一个值
[~, tIdx] = min(abs(EEG.times - 300));      % 离 300ms 最近的时间点下标
vals = ERP(:, tIdx);                        % 所有电极在该时刻的平均电位 (uV)

figure;
topoplot(vals, EEG.chanlocs, ...
    'maplimits', 'maxabs', ...              % 颜色范围对称：红=正，蓝=负
    'electrodes', 'labels', ...             % 在头皮图上标出电极名
    'style', 'both', ...                    % 填色 + 等高线
    'numcontours', 6);                      % 等高线条数
title('300 ms 的头皮电位分布');