基于 STFT 的时频分析（下）：加窗的权衡

一、窗函数：从"为什么默认就有窗"说起

你可能没注意：spectrogram 和 pwelch 默认都给数据乘了一个窗。 为什么？因为截断本身就是加窗——从连续信号上剪下一段， 等价于乘上一个"中间是 1、两端瞬间归零"的矩形窗。

瞬间归零是灾难：信号在窗两端被"硬切"，产生大量高频毛刺， 频谱上表现为能量从真实频率泼洒到邻近频点——这就是频谱泄漏。 矩形窗主瓣最窄（分辨率最好），但旁瓣最高（泄漏最凶）。

解决办法：用两端平滑过渡到零的窗，把"硬切"变"柔切"。 Hann 窗是最常用的折中——旁瓣比矩形窗低约 31dB（能量泄漏压低三个数量级）， 主瓣只宽一倍。

% 直观对比：对 10.3Hz（非整周期截断）信号做两种窗的周期图
fs = 500; t = 0:1/fs:1-1/fs;
x  = sin(2*pi*10.3*t);              % 1 秒装不下整数个周期

N = length(x);
f = (0:N/2-1) * fs / N;
P_rect = abs(fft(x.*rectwin(N))/N).^2;
P_hann = abs(fft(x.*hann(N))/N).^2;

figure;
subplot(2,1,1); plot(f, 10*log10(P_rect(1:N/2)));
xlim([5 16]); title('矩形窗：10.3Hz 附近的能量泼得到处都是');
subplot(2,1,2); plot(f, 10*log10(P_hann(1:N/2)));
xlim([5 16]); title('Hann 窗：能量收拢在 10.3Hz 峰下');

二、不确定性原理：一对不可兼得的分辨率

时频分析有一条绕不过去的物理定律：时间分辨率 × 频率分辨率 ≳ 常数。 窗越短，你越清楚"什么时候"，越糊涂"什么频率"；窗越长，反之。

用具体数字把账算明白（500Hz 采样）：

三、经验法则：高频短窗、低频长窗

既然单一窗长无法两头讨好，标准策略是按频率定制窗长—— 这正是小波变换自动完成的事情（4.4 详述），但在 STFT 框架内， 可以对低频段和高频段分别跑两次不同窗长的分析，各取所长。

• 低频（δ/θ/α，2-13Hz）：周期长（100-500ms），需要 ≥2 个周期的长窗（400-1000ms）；
• 高频（β/γ，20Hz 以上）：周期短（<50ms），150-250ms 窗已含多个周期，时间精度优先。

四、基线校正：让 1/f 不再抢戏

原始时频功率图有一个顽疾：1/f 背景。低频功率天然比高频大一到两个数量级， 画出来下面一片红、上面一片黑，实验效应完全看不见。必须做基线校正。

两种标准变换：

1. dB 变换（最常用）：每个频率上，用基线期（如刺激前 -500 到 -100ms）的平均功率作参考： P_db = 10*log10(P / P_baseline)。0 表示与基线持平，正值为增强（ERS），负值为抑制（ERD）。
2. 百分比变换：P_pct = 100 * (P - P_baseline) / P_baseline。直观（"相对基线变化百分之几"），但高频段基线极小，百分比容易爆炸。

五、完整实战：从 epoch 到可发表时频图

% 输入：epochs (ntrials x nchans x ntimes)，500Hz，epoch = [-1 1]s
fs   = 500;
time = linspace(-1, 1, size(epochs, 3));
win  = hann(round(0.4*fs));            % 400ms Hann 窗
nov  = round(0.9*length(win));
nfft = 512;

% 1) 逐试次时频变换并平均功率
ntr = size(epochs, 1);
for tr = 1:ntr
    [s, f, tt] = spectrogram(squeeze(epochs(tr, chanIdx, :)), ...
                             win, nov, nfft, fs);
    P(:, :, tr) = abs(s).^2;
end
Pavg = mean(P, 3);                     % 频率 x 时间

% 2) 逐频率 dB 基线校正（-500 ~ -100ms）
bl = tt >= -0.5 & tt <= -0.1;
Pdb = 10 * log10( Pavg ./ mean(Pavg(:, bl), 2) );

% 3) 成图
figure;
imagesc(tt, f, Pdb); axis xy;
colormap(parula); colorbar;
caxis([-2 2]);                         % 固定色标范围，跨被试可比
ylim([2 45]);
xlabel('时间 (s)'); ylabel('频率 (Hz)');
title(sprintf('通道 %d · dB 基线校正（-500\~-100ms）', chanIdx));

% 4) 保存为 300dpi 图片（论文用）
print('-dpng', '-r300', 'tfr_figure.png');

这张图已具备论文时频图的所有要素：dB 刻度、逐频率基线、固定色标、 合理频段范围。接下来 8.7 课的 ITPC 会复用同一套逐试次框架—— 只是把"平均功率"换成"平均相位向量"。