基于 STFT 的时频分析（上）：原理与实现

一、两种"看不见"：为什么需要时频分析

到目前为止我们有两大工具，但它们各有盲区：

• FFT / PSD 的盲区：整段信号压成一张频谱清单，完全丢失时间信息。你只知道"这段 10 秒里有 α 节律"，不知道它是前 3 秒有还是后 3 秒有、是持续还是阵发。
• ERP 的盲区：平均只保留锁时锁相于刺激的成分。而大脑另一大类活动——比如刺激诱发但相位随机的节律增强或抑制（ERD/ERS，统称"诱导响应"）——在平均中被正负抵消，彻底消失。

解决思路朴素而优雅：给频谱装上时间轴——用一个短窗在信号上滑动， 每个位置做一次 FFT，把每次的频谱按时间排开。

二、STFT 原理：滑动窗 + 逐段 FFT

短时傅里叶变换（STFT）的完整流程只有三步：

1. 取一个固定长度的窗（如 500ms 的 Hann 窗）；
2. 把它乘在信号当前位置，对窗内数据做 FFT，得到这一小段的频谱；
3. 窗向右滑动一段距离（如 50ms），重复，直到滑过整段信号。

把每次得到的频谱按时间排成一排，就得到一个频率 × 时间的功率矩阵。 频率是纵轴、时间是横轴、功率用颜色表示——这就是时频图（spectrogram）。

三、Matlab 实现：spectrogram 函数

先用一段"频率随时间变化"的仿真信号练手：

% 构造信号：前 1 秒 10Hz，后 1 秒突然变成 25Hz
fs = 500;
t  = 0 : 1/fs : 2 - 1/fs;
x  = [sin(2*pi*10*t(1:500)), sin(2*pi*25*t(501:end))];

figure;
subplot(2,1,1); plot(t, x); xlim([0.9 1.1]);
title('时域：1 秒处频率从 10Hz 跳到 25Hz');

subplot(2,1,2);
spectrogram(x, hann(250), 125, 512, fs, 'yaxis');
title('时频图：两条水平能量带，跳变时刻一目了然');
ylim([0 50]);

时域图上你只能看到波形突然变密；而时频图上，10Hz 和 25Hz 各自是一条水平的能量带，切换发生在 1 秒处，清楚得像红绿灯。

spectrogram 也可以用输出参数形式调用，拿到矩阵自己画：

[s, f, tt] = spectrogram(x, hann(250), 125, 512, fs);
% s  : 频率 x 时间 的复数谱；功率 = abs(s).^2
% f  : 纵轴频率向量
% tt : 横轴各窗中心时刻

P = abs(s).^2;
figure; imagesc(tt, f, 10*log10(P)); axis xy;
xlabel('时间 (s)'); ylabel('频率 (Hz)');
colorbar; ctitle_str = '功率 (dB)'; ylim([0 50]);

四、三个关键参数

逐个说透：

窗长：分辨率的天平

窗内必须容纳至少一个完整周期才能"看见"该频率。500Hz 采样下：250ms 窗的频率分辨率是 4Hz——α 和 β 的细节全糊在一起；500ms 窗是 2Hz。但窗越长，时间上越"迟钝"：快变的事件会被抹平。你不可能同时要高时间分辨率和高频率分辨率，这就是不确定性原理，4.3 展开。

重叠：让时间轴连续

若窗每次只滑动"一个窗长"，相邻两窗不重叠，时间轴像珠串一样有缝。50% 重叠是惯例；追踪快速动态时可用 90%。代价只是计算量。

nfft：只是插值，不是魔法

nfft 大于窗长的部分是补零，让频轴点更密、图画得更平滑。真实频率分辨率不变——与 4.1 的结论一致。

五、怎么读一张时频图

拿到任何时频图，按固定清单过一遍：

1. 横轴范围与零点：时间零点通常对齐刺激/事件；零之前是基线期；
2. 纵轴范围：是否只画到关心的频段（如 2-45Hz），γ 之上多半是肌电；
3. 颜色刻度：线性功率还是 dB？dB 刻度把动态范围压进可见区间，是论文标准；
4. 有没有做基线校正：未校正的图里 1/f 背景会"下亮上暗"，像瀑布（4.3 教你校正）；
5. 能量块的走向：水平条带=节律性持续活动；短促竖条=瞬态；斜向条纹多为伪迹。

六、从仿真走向 EEG：epoch 的时频分析骨架

% 输入：epochs (ntrials x nchans x ntimes)，fs = 500
% 对每个试次做 STFT，再跨试次平均（抓住诱导响应的关键！）
fs = 500; ntr = size(epochs, 1);
win = hann(round(0.4*fs));         % 400ms 窗
nov = round(0.9 * length(win));     % 90% 重叠
nfft = 512;

for tr = 1:ntr
    [s, f, tt] = spectrogram(squeeze(epochs(tr,1,:)), win, nov, nfft, fs);
    P(:, :, tr) = abs(s).^2;        % 每试次的功率
end
Pavg = mean(P, 3);                  % 跨试次平均 —— 诱导响应在这里显形

figure; imagesc(tt, f, 10*log10(Pavg)); axis xy;
xlabel('时间 (s)'); ylabel('频率 (Hz)'); colorbar;
ylim([2 45]); title('试次平均时频图');