基于时频域睡眠脑电信号自动分期与可编辑MATLAB源码实现

简介：这份资源是围绕时频域睡眠脑电信号自动分期撰写的完整技术论文PDF，面向生物医学工程、信号处理方向的学生与研究人员，帮助理解如何用小波变换与功率谱估计实现睡眠阶段自动识别。全文系统梳理了脑电信号形成机理、睡眠分期标准，并重点讲解Haar、Morlet、Meyer三种经典小波，以及相关函数法、周期图法、AR模型等功率谱估计方法，最后借助MATLAB的GUI界面展示睡眠分期统计、功率谱散点图与零点穿越图等可视化结果。资源包内仅含1个PDF文件，约195KB，篇幅紧凑，目录涵盖绪论、脑电与睡眠基础知识、小波变换与功率谱估计、MATLAB分析实现及总结等章节，结构完整，便于按模块查阅。目前已有258人学习下载，适合作为睡眠脑电自动分期入门与算法复现的参考资料。

1. 从一份能改的 MATLAB 睡眠分期源码说起

如果你正在做睡眠脑电信号自动分期，大概率遇到过两种资源：一种是论文里只有几张流程图和准确率表格，代码影子都看不到；另一种是丢给你一堆 .m 文件，打开全是加密的 .p 文件或者函数名被混淆过，想改一个滤波参数都无从下手。这份「基于时频域睡眠脑电信号自动分期和全部 MATLAB 原代码可编辑」的资源，解决的就是后一种痛点——它给的是可读、可改、可复现的 MATLAB 源码，而不是一个只能跑的黑匣子。

睡眠分期本身是脑电分析里落地性很强的任务：把整夜 PSG 信号切成 30 秒一帧，判定为 W、N1、N2、N3、REM 五类。难点在于脑电是非平稳信号，单看时域波形区分度低，所以主流做法是转到时频域提取特征，再送分类器。这份资源的技术栈就落在这条主线上，适合做生物医学信号处理的研究生、准备睡眠算法落地的工程师，以及需要一份能改的 MATLAB 基线代码来对比自己方法的人。下面我按「它怎么组织、怎么跑起来、参数怎么调、哪里会翻车」的顺序拆一遍。

2. 时频域特征提取：从原始 EEG 到分类器输入

2.1 为什么睡眠分期绕不开时频域

先讲清楚选型理由，不然后面调参全是玄学。睡眠脑电的核心节律——delta（0.5–4 Hz）、theta（4–8 Hz）、alpha（8–13 Hz）、spindle（12–14 Hz）、K 复合波——它们的区分信息几乎都藏在频率成分里。你在时域上看 N2 和 N3 的波形，都是慢波为主，肉眼很难稳定区分；但一转到功率谱，N3 的 delta 功率占比明显压过 N2。这就是为什么几乎所有自动分期方法都要做时频变换。

常见做法有三条路：短时傅里叶变换（STFT）、小波变换（WT）、以及基于滤波器组的频带能量提取。这份源码走的是可编辑的 MATLAB 实现，意味着你能看到每一步的窗长、重叠、频带边界是怎么设的，而不是被封装成一个 extractFeatures() 调用完事。对做研究的人来说，这个透明度比准确率多两个点更值钱。

2.2 预处理与分帧的可执行步骤

拿到整夜 EDF 或 MAT 格式的脑电后，第一步不是急着做变换，而是把信号规整成统一采样率和统一帧长。睡眠分期的行业惯例是 30 秒一帧，采样率常见 100 Hz 或 200 Hz。下面这段是我一般会先跑的预处理骨架，参数按你数据实际情况改：

% 读取单通道 EEG，假设已加载为变量 eeg，fs 为采样率
fs = 200;                      % 采样率，按你的数据改
epochLen = 30;                 % 每帧 30 秒，睡眠分期标准
samplesPerEpoch = fs * epochLen;

% 带通滤波 0.3-35 Hz，去掉基线漂移和工频以上噪声
[b, a] = butter(4, [0.3 35]/(fs/2), 'bandpass');
eegFilt = filtfilt(b, a, eeg);  % 零相位滤波，避免帧间相位偏移

% 按 30 秒切帧，丢弃不足一帧的尾部
numEpochs = floor(length(eegFilt) / samplesPerEpoch);
epochs = reshape(eegFilt(1:numEpochs*samplesPerEpoch), ...
                 samplesPerEpoch, numEpochs).';

逻辑说明： filtfilt 而不是 filter ，是因为睡眠脑电做特征时对相位失真敏感，零相位滤波能避免慢波被平移导致帧边界错位。 butter 的阶数选 4 是经验值，阶数太高在 0.3 Hz 附近容易数值不稳定。参数说明： epochLen 不要随便改，改成 20 秒或 60 秒会让你的结果没法和其他文献对比； fs 必须和滤波截止频率联动，截止频率超过 fs/2 会直接报错。

2.3 时频特征怎么算、怎么存

分帧之后，每一帧做时频变换并提取频带能量。STFT 的窗长是个关键参数：窗太长频率分辨率高但时间分辨率差，窗太短反过来。睡眠慢波周期在 1 秒量级，我一般用 2 秒窗、50% 重叠：

winLen = 2 * fs;               % 2 秒窗
overlap = winLen / 2;          % 50% 重叠
nfft = 2^nextpow2(winLen);     % FFT 点数取 2 的幂，加速
bands = [0.5 4; 4 8; 8 13; 13 30];  % delta/theta/alpha/beta

featMat = zeros(numEpochs, size(bands,1) + 2);
for k = 1:numEpochs
    [S, f, t] = spectrogram(epochs(k,:), hamming(winLen), overlap, nfft, fs);
    psd = mean(abs(S).^2, 2);          % 时间维平均得功率谱
    for b = 1:size(bands,1)
        idx = f >= bands(b,1) & f < bands(b,2);
        featMat(k, b) = bandpower(psd(idx));   % 各频带绝对功率
    end
    totalPower = sum(psd);
    featMat(k, end-1) = featMat(k,1) / totalPower;  % delta 相对功率
    featMat(k, end)   = featMat(k,3) / totalPower;  % alpha 相对功率
end

逻辑说明：先算功率谱再按频带积分，比直接对时域做带通滤波再算方差更稳，因为 spectrogram 已经把窗效应处理掉了。相对功率（频带功率除以总功率）比绝对功率更抗个体差异和电极阻抗差异，这是睡眠分期里很实用的一招。参数说明： bands 的边界可以按你的数据微调，但 delta 上界别超过 4 Hz，否则会把 theta 混进来； nfft 取 2 的幂纯粹是为了 FFT 速度，不影响物理意义。

特征矩阵 featMat 每行对应一帧，列是各频带功率加相对功率。这份源码可编辑的价值在这里体现得最明显——你可以直接加一列 spindle 频带能量，或者把 STFT 换成小波，而不用去猜封装函数里干了什么。

3. 分类器搭建与整夜分期流程

3.1 分类器选型：为什么先用传统方法打底

时频特征提出来之后，分类器选择上有个常见误区：一上来就上深度学习。但睡眠分期的公开数据集标注量有限，而且个体差异极大，直接上 LSTM 或 Transformer 很容易过拟合到某几个受试者。这份源码走的是传统机器学习路线，常见做法是 SVM 或随机森林打底，先跑通整条链路，再考虑换模型。

选 SVM 的理由很实际：特征维度不高（十几个频带特征），样本量几千帧，SVM 在这个规模下训练快、调参少、可解释。随机森林则对特征尺度不敏感，不用做归一化。我一般两个都跑一遍对比，谁在留出受试者上表现稳就用谁。

3.2 训练与交叉验证的代码骨架

睡眠分期最忌讳的是随机划分训练测试集——同一受试者的相邻帧高度相关，随机划分会让准确率虚高十几个点。正确做法是按受试者划分，或者用留一受试者交叉验证：

% featMat 为特征，labels 为 1-5 的睡眠期标签
subjID = ...;                  % 每帧对应的受试者编号
uniqSubj = unique(subjID);
accList = [];

for s = 1:length(uniqSubj)
    testMask = subjID == uniqSubj(s);
    Xtrain = featMat(~testMask, :);
    Ytrain = labels(~testMask);
    Xtest  = featMat(testMask, :);
    Ytest  = labels(testMask);

    % 标准化：用训练集均值方差，避免信息泄漏
    mu = mean(Xtrain); sg = std(Xtrain);
    Xtrain = (Xtrain - mu) ./ sg;
    Xtest  = (Xtest  - mu) ./ sg;

    mdl = fitcecoc(Xtrain, Ytrain, 'Learners', 'svm');
    pred = predict(mdl, Xtest);
    accList(end+1) = mean(pred == Ytest);
end
fprintf('留一受试者平均准确率: %.2f%%\n', mean(accList)*100);

逻辑说明： fitcecoc 是 MATLAB 里做多分类 SVM 的标准入口，五类睡眠期用它比二分类 SVM 手动拼要省事。标准化必须用训练集的 mu 和 sg 去变换测试集，如果对全体数据先标准化再划分，测试集信息就泄漏进训练了，这是血泪经验里最常见的一类错误。参数说明： 'Learners','svm' 可以换成 'tree' 试随机森林； accList 按受试者记录，能看出模型是不是对某类人特别差。

3.3 整夜分期的后处理

逐帧分类出来的标签往往有毛刺——比如 N2 中间突然蹦出一帧 W，这在生理上几乎不可能。加一层平滑能明显提升整夜一致性：

% 中值滤波平滑标签序列，窗长 5 帧
labelsSmoothed = medfilt1(double(pred), 5);
% 睡眠期转换规则：W 不会孤立出现在 N3 中间
for k = 2:length(labelsSmoothed)-1
    if labelsSmoothed(k) == 1 && ...
       labelsSmoothed(k-1) == 4 && labelsSmoothed(k+1) == 4
        labelsSmoothed(k) = 4;   % 孤立 W 帧改判为 N3
    end
end

逻辑说明： medfilt1 对标签序列做中值滤波，能压掉单帧跳变。后面那段规则是睡眠分期里常用的生理约束——睡眠是连续过程，不会在深睡里插一帧清醒。参数说明：中值窗长别超过 5，太长会把真实的期转换也抹平；规则可以按 AASM 标准继续加，比如 REM 后通常接 N2 而不是直接跳 N3。

4. 避坑与排查：跑不通时先看这几条

4.1 现象：滤波后信号幅度异常或出现振荡

原因： butter 阶数过高加上低截止频率（0.3 Hz）在低采样率下数值不稳定，或者用了 filter 导致相位偏移叠加。解决：把阶数降到 2–4，确认 fs 至少是最高截止频率的 2.5 倍以上，滤波一律用 filtfilt 。如果还振荡，检查信号里有没有 NaN， filtfilt 遇到 NaN 会整段污染。

4.2 现象：准确率高达 95% 以上，但换个人就崩

原因：随机划分训练测试集导致同受试者相邻帧泄漏，或者标准化用了全体数据。解决：改成按受试者划分，标准化只用训练集统计量。如果准确率还是虚高，检查标签是不是严重不平衡——N2 通常占一半，全预测 N2 也能有 50% 准确率，这时候要看每类召回率而不是总准确率。

4.3 现象： spectrogram 报错或输出维度对不上

原因： nfft 小于窗长，或者 overlap 大于等于 winLen 。解决： nfft 取 2^nextpow2(winLen) 保证不小于窗长； overlap 严格小于 winLen ，50% 重叠就写 winLen/2 。另外注意 spectrogram 返回的 f 是列向量，做频带索引时方向要对齐。

4.4 现象：整夜跑完内存爆掉

原因：把整夜多通道信号一次性读进内存做 STFT，或者 featMat 预分配维度算错。解决：按帧循环处理，每帧算完就存特征，不要保留中间 S 矩阵； featMat 用 zeros(numEpochs, nFeat) 预分配，别在循环里动态增长数组，MATLAB 动态增长是性能杀手。

4.5 现象：换自己的数据后标签全乱

原因：标签编码和源码假设不一致，比如源码假设 1=W、2=N1…，你的数据是 0=W 或字符串标签。解决：训练前先 unique(labels) 看一眼实际取值，用映射表统一到源码的编码。这一步看着简单，但翻车的人特别多，因为报错往往出现在 fitcecoc 而不是数据加载处，容易误判成模型问题。

5. 进阶：把小波变换接进同一条流水线

传统 STFT 的窗长固定，对睡眠纺锤波这种短时高频事件和慢波这种长时低频事件没法同时兼顾。想再往上走一步，可以把特征提取那一段换成小波变换，其余流程不动。MATLAB 里 cwt 或 modwt 都能用，我一般用 modwt 因为它是非抽取的，帧长不会变，接后面的分类器最省事：

wname = 'db4';                 % 小波基，db4 在脑电里常用
level = 5;                     % 分解层数，对应到约 0.5-30 Hz
for k = 1:numEpochs
    w = modwt(epochs(k,:), wname, level);
    % 各层小波系数能量作为特征
    for L = 1:level
        featMat(k, nBase + L) = sum(w(L,:).^2);
    end
end

逻辑说明： modwt 不做下采样，每层系数长度和原信号一致，所以能量特征直接可算，不用对齐长度。 db4 的消失矩适中，对脑电这种非平稳信号够用。参数说明： level 和采样率联动，200 Hz 下 5 层大致覆盖到 3 Hz 左右，想覆盖 delta 要加到 6 层；小波基换成 sym8 或 coif3 可以对比，但别频繁换，换一次整套特征分布就变了。

验证方法上，我习惯做两件事：一是把新特征和原 STFT 特征拼在一起，看留一受试者准确率有没有稳定提升，而不是只看单次划分；二是画混淆矩阵，重点看 N1 的召回率——N1 是五类里最难分的，如果 N1 召回率没动，说明新特征没抓到关键信息，提升可能来自过拟合。

还有个具体技巧：整夜分期结果最好和原始信号叠在一起画出来，肉眼扫一遍期转换是否连续。我踩过的坑是模型指标好看，但画出来发现 REM 段被切得七零八落，回头查是平滑窗长设太大把 REM 的快速眼动特征抹掉了。从那以后我每次改完特征或后处理，都强制把预测标签和原始 EEG 叠图走一遍，指标只是参考，生理合理性才是底线。希望帮到你。

本文还有配套的精品资源，点击获取