简介面向生物医学信号处理与睡眠医学研究者的完整技术文档系统阐述基于时频域分析的睡眠脑电EEG自动分期方法并附有可编辑的MATLAB实现代码。资源为1个PDF文件体积约195KB内容覆盖小波变换Haar、Morlet、Meyer与功率谱估计相关函数法、周期图法、AR模型核心原理以及MATLAB GUI可视化的完整源码。文档从脑电信号基础、睡眠分期标准NREM/REM、小波阀值去噪到现代功率谱估计算法逐层展开读者可据此复现基于Yule-Walker与Burg算法的功率谱散点图、零点穿越图、睡眠分期统计图并掌握GUI交互分析流程。已有258人浏览学习适合需要开展睡眠分期算法研究、信号处理课程设计或临床辅助诊断工具开发的用户参考使用可复制代码进行二次开发与扩展。1. 睡眠自动分期第一课先弄清30秒脑电片段和标签是怎么对齐的睡眠脑电自动分期这个方向最容易让人误判的不是模型选型而是数据准备。很多下载到一套号称“全部MATLAB原代码可编辑”的分期工程后第一反应是直接运行结果卡在EDF读取、通道命名、标签偏移这类看似琐碎的问题上。自动分期的输入是30秒一个epoch的脑电信号输出是W、N1、N2、N3、REM五类中的一类。医生肉眼判断靠的是alpha节律、theta波、delta波、睡眠纺锤波这些成分在时间轴上的此消彼长而这些成分的定义天生就是频段的所以“时频域”是绕不开的技术底座。这套东西适合生物医学工程研究生、睡眠中心数据处理人员和需要批量处理PSG数据的算法工程师。下面按一个可复现的MATLAB工程来拆。2. 时频特征怎么取STFT、小波包和HHT的选型与MATLAB实现2.1 三种时频方法的边界为什么主流分期代码都选STFT或小波包睡眠脑电是典型的非平稳信号30秒内可能同时存在低频慢波和高频纺锤波只用FFT看整段频谱会丢掉这些成分出现的时间信息只看时域波形又很难定量描述节律强度所以需要时频分析。常见的三种工具是短时傅里叶变换STFT、连续小波变换CWT和小波包分解WPD。STFT的思路是加窗切片对每一小段做FFT。它的核心矛盾是窗长窗越长频率分辨率越高但时间分辨率越低。睡眠分期的常用做法是30秒一个epoch特征分析时再用4秒窗、1秒步长在epoch内滑动这样能保留节律在30秒内的变化过程。MATLAB里用spectrogram函数即可优点是计算快、谱图直观缺点是频率分辨率在分析短片段时受限。CWT和WPD解决了STFT的固定分辨率问题低频用长窗、高频用短窗匹配睡眠脑电中delta慢波与sigma纺锤波共存的场景。其中WPD比CWT更适合做特征工程因为它把信号分解到一组固定频带的子带上正好对应临床分期规则里定义的delta、theta、alpha、sigma、beta频段后续按子带累加能量非常直观。HHT这类方法虽然对非平稳信号更“自适应”但EMD分解在批量处理时速度慢、模态混叠不稳定不同受试者的IMF数量还不一致特征矩阵难以对齐实际分期工程里很少作为首选。所以落到MATLAB代码上我会把WPD作为主特征提取路径STFT作为验证和可视化辅助。这样既保证特征有生理学解释又不会让代码在“参数调不出来”上耗时间。2.2 睡眠EEG中真正有效的4类时频特征睡眠分期规则依据的核心节律在频段上是固定的。delta波0.5到4Hz对应N3深睡theta波4到8Hz对应N1和REMalpha波8到12Hz出现在清醒闭眼和N1早期sigma波12到16Hz是N2期睡眠纺锤波的标志beta波16到30Hz在清醒和REM期活跃。基于这些频段特征设计要抓四类信息。第一类是各频段相对能量比。以30秒epoch为基本单元计算delta、theta、alpha、sigma、beta五个频段的能量分别占全频段总能量的比例。N3期的delta相对能量会显著升高W期alpha和beta占主导N2期sigma能量是亮点。第二类是谱质心也就是把整个epoch的功率谱按频率加权求平均得到“主轴频率”。它用单一数值概括当前脑电的主导节律清醒时偏高深睡时偏低对区分阶段很有帮助。第三类是sigma纺锤波能量。睡眠纺锤波是12到16Hz的短时爆发出现频率高是N2期的标志性特征。单看整个30秒的sigma平均能量会淹没纺锤波的爆发特性所以可以在4秒窗内找sigma频段的局部峰值取峰值强度作为特征。第四类是theta与beta的比值。theta/beta比值反映大脑从警觉转向困倦的程度在N1和N2期显著升高。做特征时加一个稳定项防止除以零。我一般会在这些基础上再补一两个衍生特征比如谱边缘频率SEF90即90%功率集中在其下的频率点。特征数量控制在10到20个以内维度太高对分类器反而是负担。2.3 用MATLAB实现小波包特征提取最小可跑代码与参数含义下面这段函数是整套分期代码里最核心的部分。输入x是单导联30秒脑电信号向量Fs是采样率输出是一个特征向量。function feat extractWPDfeatures(x, Fs, level) % x: 单导联30秒信号已去均值单位uV % Fs: 采样率常见200Hz或256Hz % level: 小波包分解层数推荐6 x x - mean(x); wpt wpdec(x, level, db4); nodes 2^level; e zeros(nodes, 1); for i 0:nodes-1 c wpcoef(wpt, i); e(i1) sum(c.^2) / length(c); end e e / (sum(e) 1e-12); eF e(wpfrqord(nodes) 1); bandWidth (Fs/2) / nodes; bands [0.5 4 4 8 8 12 12 16 16 30]; feat zeros(1, 5); for b 1:5 lo ceil(bands(2*b-1) / bandWidth) 1; hi floor(bands(2*b) / bandWidth) 1; lo max(lo, 1); hi min(hi, nodes); feat(b) sum(eF(lo:hi)); end freqAxis (0:nodes-1) * bandWidth; feat(6) sum(freqAxis .* eF) / (sum(eF) 1e-12); feat(7) feat(2) / (feat(1) 1e-8); end这段代码里有几个必须说清楚的点。首先是小波包分解后子带的排列顺序不是按频率从小到大MATLAB的wpdec返回的是Paley序直接按序号取子带会发现频率是乱序的。wpfrqord函数负责把Paley序转换成自然频率序这一步漏掉的话后面累加delta能量会算到别的频段上特征完全错乱。其次db4小波是Daubechies族里常用的基函数时域支撑短、计算快在睡眠分期文献里出现频率很高不必换成更长的小波。level取6时在Fs200Hz下每个子带带宽是1.5625Hzdelta频段大约覆盖三个子带频率分辨率足够支撑分期取5层带宽变成3.125Hzsigma频段只覆盖约1.3个子带能量计算会有毛刺。参数调优方面我的建议是level固定为6不要为了追求更细的频率分辨率盲目加到8或9。子带越多特征维度虽然可以由你手动聚合成5个频段但小波包分解的计算量成倍增加而且过细的子带对分类器没有额外信息量。另外wpcoef返回的是每个子带的系数用系数平方和近似子带能量的做法在小波基是正交基的前提下是成立的所以用小波包而不是连续小波变换。提示在跑这段代码前先确认MATLAB安装了Wavelet Toolbox。没有这个工具箱wpdec和wpfrqord都调用不了。30秒epoch内部如何应用这段函数我常用的做法是对每个epoch按4秒窗长、1秒步长滑窗共得到27个小窗对每个小窗调用extractWPDfeatures然后把27个特征向量按列做均值再补一列方差。这样做的好处是保留了节律在epoch内部的时变信息同时抑制了短时伪迹对单窗特征的冲击。代码上的代价是每个epoch要跑27次小波包分解批量处理时会稍慢但换来的是N1和N3的判别力明显提升。3. 分期主流程跑通从原始EDF到五分类结果的最小闭环3.1 数据切片与标签对齐一个epoch偏移就会让结果整体漂移原始PSG数据通常以EDF格式存储里面除了脑电还有眼电、肌电、心电等通道。自动分期最常见的做法是只取单导联脑电如C4-A1或F4-A1因为睡眠分期规则本身就是建立在单通道脑电上的。% 读取EDF文件老版本MATLAB的edfread返回hdr和record [hdr, rec] edfread(subj01.edf); Fs hdr.samples(1); % 确认通道顺序找到C4-A1或F4-A1所在列 chNames hdr.label; idx find(contains(chNames, C4)); eeg rec(:, idx(1)); % 物理单位校准EDF里没有统一单位常见uV或ADC if hdr.units{idx(1)} uV eeg eeg; else eeg eeg * hdr.calib(idx(1)); end eeg eeg - mean(eeg); % 按30秒一个epoch切分 epochLen Fs * 30; n floor(length(eeg) / epochLen); epochs reshape(eeg(1:n*epochLen), epochLen, n);这段代码有几点要特别注意。一是edfread在不同MATLAB版本里的返回值不一样R2020a及之前返回hdr和record两个变量之后版本可能返回timetable或对象格式拿到原代码后先确认它适配的是哪种版本。二是EDS里通道名没有统一标准有的写C4-A1有的写EEG C4必须用contains做模糊匹配找不到就直接报错列出所有通道名。三是物理单位EDF文件头里可能标着uV但实际数据是ADC原始值不做校准的话后面幅度阈值和能量计算全都不对。标签对齐是这里最隐蔽的坑。睡眠中心的标注软件通常每隔30秒输出一个分期标签但起始点可能不是0秒。例如有些软件把记录开始后的第1秒到第30秒标为第一个epoch有些则从第0秒开始。如果代码里直接按0到30秒切分而标签文件按1到31秒切分那么每个epoch的脑电内容就和标签错位了一个采样点。这种错位不会让训练报错只会让准确率神秘地停在70%上下且怎么调都上不去。解决方法是先读取标签文件的时间戳数组用diff检查相邻标签间隔是否为30秒再与信号时间轴做同步。3.2 随机森林还是SVM特征维度不高时先别上深度网络当特征是前面提到的7到20维时随机森林和SVM都是比深度学习更合适的选择。深度网络需要大量数据来拟合时序依赖而常见睡眠数据集一个受试者一夜也只有约1000个epoch样本量不足以支撑复杂模型。随机森林的优势在于对特征尺度不敏感不需要做标准化能直接输出特征重要度。SVM在特征维度低时训练极快但必须做z-score标准化且对类别不平衡更敏感。MATLAB里没有单独叫fitcrandomforest的函数随机森林是fitcensemble中用Bag方法加决策树学习器实现的。SVM多分类则用fitcecoc内部自动做一对一编码。我这里推荐先用随机森林打底。原因是睡眠分期特征之间有明显的非线性关系比如delta相对能量高时theta相对能量低这种关系用树模型能自然捕捉。而且树模型对离群点有容忍度脑电伪迹不可能完全剔除干净。超参数设置上树的棵数取100到200棵足够再增加只拖慢训练速度准确率几乎不动。MinLeafSize是防止过拟合最关键的参数默认值是1但睡眠分期任务里我一般取5到10。MinLeafSize太小树会深挖个别受试者的特殊模式跨受试者验证时掉点明显。类别不平衡方面N2通常占比最高N1最少可以在fitcensemble里设置Prior,uniform来弱化多数类主导。3.3 用MATLAB完成训练与交叉验证按受试者分组是底线交叉验证这里必须强调一个原则训练集和测试集不能包含同一个受试者的epoch。睡眠相邻epoch之间的特征高度相似如果同一受试者前半夜的epoch进了训练集、后半夜的进了测试集模型等于变相看到了测试数据评估结果会虚高。MATLAB里实现按受试者分组的交叉验证用cvpartition的group版本。rng(2025); cvp cvpartition(subjectID, KFold, 5); acc zeros(cvp.NumTestSets, 1); confMat zeros(5, 5); for k 1:cvp.NumTestSets tr cvp.training(k); te cvp.test(k); % 标准化只拟合训练集 mu mean(X(tr, :)); s std(X(tr, :)); Xtr (X(tr, :) - mu) ./ (s 1e-8); Xte (X(te, :) - mu) ./ (s 1e-8); mdl fitcensemble(Xtr, y(tr), ... Method, Bag, ... NumLearningCycles, 150, ... Learners, templateTree(MinLeafSize, 5), ... Prior, uniform); pred predict(mdl, Xte); acc(k) mean(pred y(te)); confMat confMat confusionmat(y(te), pred); end这段代码单独拎出来解释两个细节。第一标准化均值方差的计算只能发生在训练集内部测试集直接套用训练集的mu和s。很多新手把整个X先标准化再切分这会让每折验证都泄露测试集统计信息。第二cvpartition的group参数传入的是受试者编号向量每个epoch对应一行同编号的epoch永远不会跨入训练集和测试集。如果不这么做后面看混淆矩阵会得到一个过于乐观的准确率实际部署时立刻垮掉。训练结束后把每折的混淆矩阵累加计算每类的precision和recall。睡眠分期任务的评估不能只看整体准确率因为N2占比可能超过40%一个全部预测为N2的模型也能有40%以上的准确率看起来不太差但毫无用处。要额外盯住N1的召回率它是所有分类器最难啃的阶段。4. 自动分期最常见的5个翻车现场与排查方法4.1 现象N1的召回率只有二成大量N1被判断成N2或W原因N1本身是过渡期特征上既不像清醒期那样有明显alpha节律也不像N2那样有稳定的纺锤波专家标注N1的一致性本身就低。另外单个30秒epoch的特征只描述当前片段内部信息没有考虑前后阶段的上下文而N1的判定恰恰依赖它出现在W之后、N2之前的这个时序位置。解决给特征加上下文。把第k-1和第k1个epoch的特征向量直接拼接到当前epoch特征后面让分类器看到“之前是什么、之后是什么”。如果特征维度是20拼接后变成60随机森林依然能处理。另一个做法是预测完所有epoch后用滑动窗口对预测概率做平滑因为睡眠阶段转移不可能一秒内从W跳到N3。4.2 现象训练集准确率95%验证集只有73%感觉像换了套数据原因几乎可以确定是数据泄漏。第一种泄漏是全局标准化把全部epoch的均值和方差算完后才划分训练测试集测试集的统计信息提前进入了模型。第二种泄漏是同一受试者的epoch被随机分到了训练集和测试集相邻epoch高度相关模型等同于见过了测试样本的“邻居”。解决标准化改写成交叉验证循环内部只fit训练集交叉验证一律用cvpartition(subjectID, KFold, 5)按受试者分组。改完这两个地方验证集准确率会下降一点但那是真实水平。4.3 现象A受试者模型表现很好换到B受试者直接垮掉原因不同受试者的脑电幅度、节律分布、电极位置都存在差异。有的受试者alpha节律很强有的几乎看不到老年人慢波普遍偏少N3阶段特征与青年人有差异。如果模型只在少数几个受试者上训练泛化能力必然受限。解决先检查训练集是否覆盖了不同性别、年龄段的人群。特征层面可以按受试者做基准校正取该受试者入睡前的清醒期信号计算alpha能量基线之后所有特征减去这个基线消除个体alpha强度的系统性差异。如果数据量实在有限至少要在报告里写明模型适用的受试者范围。4.4 现象小波包特征在每个epoch边界处突变N3误报特别多原因小波分解对信号边界敏感默认的边界延拓模式是零填充信号两端会出现不连续导致分解系数在边界处产生虚假高能量。epoch边界正好是特征计算的起点和终点虚假能量会污染整个频段能量比。解决在调用wpdec前切换边界延拓模式或在特征提取时丢弃每个窗口两侧的系数。常用做法是执行dwtmode(symw)把默认延拓改为对称延拓然后对4秒窗提取特征时只保留中间3秒的窗输出丢弃前后各半秒。另外伪迹容易出现在epoch边界因为医生打标时不会把伪迹精确截掉。4.5 现象epoch内有大量伪迹能量特征被整体抬高W期被误判成REM原因头部运动、电极接触不良、快速眼动都会产生大幅值信号在频域上表现为全频段能量抬升。如果伪迹幅度太大小波包子带能量会整体失真任何依赖绝对能量的特征都会失效。解决在特征提取前做两层过滤。第一层是幅度阈值|x|超过500uV的采样点所在的epoch直接标记为伪迹第二层是差分阈值相邻采样点之差超过某一阈值说明存在高频突变。伪迹epoch在训练阶段直接剔除在测试阶段可以标记为无法判断而不是强行给一个分期结果。5. 把“全部MATLAB原代码可编辑”改造成真正能调的工程5.1 先确认PDF里的代码到底能不能编辑很多以PDF形式分发的MATLAB代码关键问题不是算法而是代码层是不是文本。如果PDF里每一行代码是图片那“可编辑”就是一句空话只能手动重敲。拿到资料后第一件事是用鼠标框选一段代码试试能否选中复制。能选中意味着代码层是文本后续才谈得上改造。其次要看代码里是否硬编码了数据路径例如直接写着load(C:\Users...\data.mat)这类代码换一台机器必然报错需要全局搜索.mat和load。代码完整性检查也有套路。查看主脚本末尾是否有disp或fprintf之类的输出提示如果有说明作者至少自己跑通过。查看是否有functions文件夹或明显的子函数没有的话所有代码挤在一个脚本里可维护性较差。查看注释密度好代码的注释会告诉你每个参数在生理上对应什么烂代码只有光秃秃的变量名。提示正文里的代码块可以直接复制到MATLAB的.m文件里。PDF如果是从网页转换来的要小心代码中英文引号被替换成全角引号运行时会报“未定义变量”或字符串错误。5.2 建立参数集中化文件以后改参数只动一个文件把散落在主脚本各处的参数集中到一个params.m文件里是所有后续改造的第一步。睡眠分期需要用到的参数就十来个集中管理后想换数据集、换电极、换小波基函数都只改一处。% params.m 集中管理所有可调参数 Fs 200; % 采样率HzEDF头文件可查 EPOCH_SEC 30; % 每个epoch时长AASM标准固定30秒 CHANNEL_KEYWORD C4; % 通道名模糊匹配关键词 FILTER_BAND [0.3 35]; % 带通滤波范围Hz去除直流和超高频 NOTCH_FREQ 50; % 工频陷波频率中国电网50Hz北美60Hz WPD_LEVEL 6; % 小波包分解层数 WAVELET db4; % 小波基函数 SLIDE_WIN 4; % 滑动窗口秒数 SLIDE_STEP 1; % 滑动步长秒数 TRAIN_METHOD bag; % 分类器类型bag或svm NUM_TREES 150; % 随机森林树棵数 MIN_LEAF 5; % 叶子节点最小样本数 PROB_SMOOTH_WIN 5; % 预测概率平滑窗口单位epoch给这段参数文件两个建议。一是每个参数都写注释说明取值范围和修改后果这会在一个月后你重新拿起代码时救你一命。二是陷波频率要注意电网差异国内EDF数据普遍是50Hz工频但有些公开数据集来自北美是60Hz用错了陷波频率会把脑电信号本身也吃掉一部分。5.3 把特征提取、训练、预测拆成独立函数原版代码如果是上千行堆在main.m里先别急着全盘重写而是按功能边界拆成四个函数。第一个是数据加载函数负责读EDF、切epoch、读标签返回一个standardData包。第二个是特征提取函数接收epoch和params返回特征矩阵。第三个是训练函数接收特征矩阵和标签返回训练好的模型。第四个是预测函数接收模型和待预测特征返回分期序列和概率矩阵。拆成函数后每个函数可以单独验证。比如特征提取函数写完后输入一段已知有纺锤波的N2期数据输出特征里sigma频段相对能量应当明显偏高如果不符合这个预期说明小波包子带映射写错了。训练函数写完后用一个只有50个epoch的小数据集跑通确认没有行列维度问题再上全量数据。这个习惯能节省大量调试时间。6. 不花算力的提分技巧前后文特征拼接与概率平滑睡眠分期和其他分类任务最大的区别是强时序依赖。医生看一个30秒epoch时脑子里装着前一个epoch是什么、后面会过渡到什么阶段。把这种时序常识注入模型的常见技巧是前后文特征拼接也就是把第k-1、k、k1三个epoch的特征向量拼接成第k个样本的特征。这个技巧不需要增加任何计算资源没有新数据也没有新特征只是把特征矩阵从X变成[X(k-1,:), X(k,:), X(k1,:)]随机森林就能捕捉到“从W进入N1再进入N2”的过渡模式。% 前后文特征拼接示例 n size(X, 1); Xc zeros(n, size(X, 2) * 3); for k 2:n-1 Xc(k, :) [X(k-1, :), X(k, :), X(k1, :)]; end % 首尾epoch无法拼接直接复制单epoch特征 Xc(1, :) [X(1, :), X(1, :), X(2, :)]; Xc(n, :) [X(n-1, :), X(n, :), X(n, :)];概率平滑比硬标签平滑更可靠。MATLAB的predict函数对fitcensemble返回的是最终分类标签要拿到每个epoch属于五类的概率调用predict(model, Xte, Learners, 1:NumTrained)这种方式比较麻烦更通用的做法是predict函数输出[labels, scores]scores就是每个类别的分数。把scores按epoch顺序做移动平均再取每行最大分数对应的类别作为最终分期。[~, scores] predict(mdl, Xte); % 对概率做移动平均窗口取5个epoch scoresS movmean(scores, 5, 1); [~, predSmooth] max(scoresS, [], 2);平滑窗口不是越大越好。取5是合理的相当于前后各看2个epoch取15以上会把短暂的N1完全抹掉反而降低总体准确率。验证这个技巧有没有效的方法很简单比较平滑前后混淆矩阵中N1的召回率以及hypnogram里是否出现不合理的瞬间切换比如W直接切N3再切回W。我现在的习惯是每跑完一批数据先看两样东西一个是混淆矩阵里N1那一行另一个是平滑前后hypnogram的差异。如果N1还是丑先别动模型而是回头查标签切分有没有偏移、边界延拓方式改没改对这两个基础问题修好之后比换任何分类器都涨得多。希望帮到你。本文还有配套的精品资源点击获取