功率谱密度（PSD）计算

一、为什么要看频率

想象你在听一首歌：示波器看到的是一条不断抖动的空气压力曲线， 但你的耳朵听到的是"低音、中音、高音"。耳朵做的事情就是 频率分解——把混合在一起的波形拆成各个频率的成分。

脑电同理。头皮记录到的信号是无数神经元群体同步振荡的混合物： 有 10Hz 左右的 α 节律、有 20Hz 的 β 活动、也有慢悠悠的 δ 波。 ERP 平均会把不锁时于刺激的节律"平均掉"，而这些节律恰恰携带了 大脑状态的关键信息——所以我们需要频域分析。

二、傅里叶变换：拆解波形的"三和弦"

傅里叶变换的神奇之处在于：任何一段合理信号，都可以唯一地 表示成一组不同频率、不同幅度、不同相位的正弦波的叠加。就像一个和弦—— 你听到的是 C、E、G 三个音的混合，而傅里叶变换能精确告诉你 "这个和弦里有哪几个音、每个音多响"。

先用一段仿真信号建立直觉：

% 构造一个"三和弦"信号：5Hz + 10Hz + 40Hz 三个正弦波叠加
fs = 500;                    % 采样率 500 Hz
t  = 0 : 1/fs : 2 - 1/fs;    % 2 秒时间轴
x  = 3*sin(2*pi*5*t) ...     % 5Hz，幅度 3
   + 7*sin(2*pi*10*t) ...    % 10Hz，幅度 7（最强）
   + 2*sin(2*pi*40*t);       % 40Hz，幅度 2

figure; plot(t(1:250), x(1:250));   % 看前 0.5 秒
xlabel('时间 (s)'); title('时域：一团乱麻？');

时域上这就是一团复杂的波形，肉眼看不出任何结构。但做一次 FFT：

N  = length(x);
X  = fft(x) / N;            % 除以 N 得到归一化幅值
f  = (0:N-1) * fs / N;      % 每个点对应的频率

half = 1 : floor(N/2);      % 只取一半（对称镜像）
figure; stem(f(half), abs(X(half)), 'filled');
xlabel('频率 (Hz)'); ylabel('幅度');
title('频域：清清楚楚三个峰 —— 5Hz(3) / 10Hz(7) / 40Hz(2)');

频谱上只有三个峰，分别在 5、10、40 Hz，高度恰好是 3、7、2。 混乱消失了——这就是频率视角的力量。

三、从频谱到功率谱

幅度谱（|X(f)|）衡量每个频率"有多强"，但物理和统计上我们更常用 功率——也就是幅度的平方：

P = abs(X).^2;              % 功率谱 = 幅度谱的平方
% 平方有两个好处：
% 1) 与能量直接挂钩（物理意义）
% 2) 让小信号更小、大信号更大，对比更鲜明

而当你把信号长度增加一倍、把相邻频点的功率加起来归一化，就得到 功率谱密度（PSD, Power Spectral Density），单位是 功率/Hz（如 μV²/Hz）。它描述的是"每个频率上、每 1Hz 带宽内分布着多少功率"， 与数据长度无关，因此可以在不同研究间直接比较。

四、周期图与 Welch 方法：为什么要"分段平均"

最直接的 PSD 估计叫周期图（periodogram）：整段数据做一次 FFT 再平方。 它的问题在于方差极大——对同一段平稳信号反复截取不同片段算周期图， 得到的谱会剧烈抖动，像噪声一样不可靠。

Welch 方法是标准解法，三步走：

1. 分段：把长数据切成若干段（可重叠）；
2. 加窗：每段乘一个窗函数（默认 Hann 窗），压住截断造成的泄漏；
3. 平均：把各段的功率谱平均，方差被压低到约 1/K（K 为段数）。

五、Matlab 标准实现：pwelch

% pwelch 是科研中最常用的 PSD 估计函数
fs = 500;
[pxx, f] = pwelch(eeg_data, ...
    hann(500), ...          % 窗：500 点 Hann 窗（约 1 秒）
    250, ...                % 相邻段重叠 250 点（50%）
    2048, ...               % nfft：补零到 2048 点，插密频轴
    fs);                    % 采样率

figure; plot(f, 10*log10(pxx));    % 习惯画成对数刻度（dB）
xlabel('频率 (Hz)'); ylabel('功率谱密度 (dB, re 1 \muV^2/Hz)');
xlim([1 80]);                     % 只看有生理意义的频段

参数选择的门道：

六、计算 EEG 功率谱的五个注意事项

1. 数据长度决定频率分辨率

分辨率 = fs / N。1 秒数据只有 1Hz 分辨率；想精细看 α 峰的形态，至少 2-4 秒连续数据。静息态分析常用 2 秒窗。

2. Nyquist 限制：一半以上的频率是"镜像"

采样率 500Hz 的数据，有效频率上限是 250Hz。频谱上半部分只是下半部分的镜像回折，画图时通常只画到 fs/2。

3. 1/f 背景趋势

真实 EEG 功率谱不是平的：低频功率远大于高频（近似 1/f 规律），叠加慢漂移后更明显。直接看绝对功率会"低频独大"。对策：分段前去趋势（detrend）或高通滤波 0.5-1Hz，或分析时改用相对功率 / 对数刻度。

4. 典型频段与生理含义

5. epoch 数据 vs 连续数据

同样的分段窗（如 2 秒），从连续数据切窗和从预处理后的 epochs 切窗，结果口径不同——epochs 已经过基线校正和伪迹剔除，通常更干净；但若在 epochs 边缘切窗会有截断效应。报告方法时必须写清楚用的是哪种、窗长多少、重叠多少。

七、把一课串起来：单通道 EEG 的 PSD 完整流程

% 输入：eeg_cont (1 x Ntimes) 连续数据，fs = 500
fs = 500;
eeg_det = detrend(eeg_cont');          % 去趋势，压掉慢漂移

[pxx, f] = pwelch(eeg_det, hann(2*fs), fs, 4096, fs);

% 提取各经典频段的平均功率
bands = struct( ...
  'delta',  [1 4], 'theta', [4 8], ...
  'alpha',  [8 13], 'beta', [13 30], 'gamma', [30 45]);
bn = fieldnames(bands);
for k = 1:numel(bn)
    idx = f >= bands.(bn{k})(1) & f <= bands.(bn{k})(2);
    pw(k) = mean(pxx(idx));            % 该频段平均 PSD
end
disp(table(bn', pw', 'VariableNames', {'band','power'}));

figure; plot(f, 10*log10(pxx)); xlim([1 50]);
xlabel('频率 (Hz)'); ylabel('PSD (dB)');

这段代码就是静息态 EEG 频谱分析的骨架——把它包一层被试循环， 就是你在无数论文方法部分看到的"各频段功率"指标。