企业官网建设流程全解析

简介：面向脑电信号处理与睡眠医学研究者的MATLAB源码工具，通过计算脑电信号的功率能量分布辅助判断睡眠阶段，适用于睡眠分期算法验证、课程设计或科研预研场景。压缩包共5个文件，含3个M脚本与2个ASV自动存档文件，整体仅3KB，代码精简，按功能拆分为信号读取、功率谱计算、小波变换等模块，便于快速阅读与二次开发。目前已有288人学习并下载。源码覆盖从脑电信号读取、功率谱计算到趋势分析与特征判别的完整环节，其中小波变换模块可供对比不同频带能量变化，帮助使用者理解REM、浅睡、深睡等阶段在频谱上的差异，同时可作为算法改进的起点。对于刚接触MATLAB信号处理或睡眠监测的读者，是一份轻量、直观的入门素材。

1. 脑电波睡眠监测的功率能量判断法：为什么先算谱而不是先看波形

一段整夜脑电按 100 Hz 采样，六小时就是 216 万个点。直接看波形去数纺锤波和慢波，人眼盯不住，程序也跑不动。睡眠分期的本质不是看波形长什么样，而是看不同频段的能量往哪集中：清醒闭眼时 alpha 带抬起来，N2 阶段纺锤波让 sigma 带冒尖，深睡 N3 则是慢波占主场。把脑电切成一帧帧 30 秒的 epoch，对每帧做功率谱估计，再用频带能量占比去贴 AASM 分期标签，这就是脑电波睡眠监测根据功率能量判断睡眠阶段的工程落地路径。适合做睡眠监测算法的工程师，也适合想把脑电数据快速跑成分期曲线的研究生：不需要商用分析软件，MATLAB 加信号处理工具箱就够。

2. 睡眠分期标准与五频带能量签名：AASM 规则下的功率分布

2.1 先定裁判：AASM 五类标签与 30 秒 epoch

做睡眠分期第一步不是写代码，而是定“裁判标准”。临床睡眠分期的主流依据是 AASM 规则：整夜记录被切成连续的 30 秒片段，每段打一个标签，标签有五类：W（清醒）、N1（入睡期）、N2（浅睡）、N3（深睡）、REM（快速眼动期）。AASM 判读依赖多导信号——至少一个中央区脑电通道（C3 或 C4）、两个眼电通道、一个下颌肌电通道。脑电功率谱反映的只是“皮层状态”，而 REM 的快速眼动和 W 的高肌电张力要借助 EOG 和 EMG 判断。

这直接划定了纯脑电功率方法的边界：单通道 C3 的频带能量，能稳妥区分 N3、N2 和大部分 N1；W 与 REM 在频谱上高度相似，单通道几乎无解。后面所有实现都围绕这个边界设计。先把这条边界说清楚，看到 W/REM 混淆分数难看时，就不会怀疑是代码写错，而是明白这是信息量不足的必然结果。

2.2 五个频带的睡眠分期能量签名

脑电功率谱在 0.5–30 Hz 内的分布形态是分期的核心依据。临床惯用的频带划分如下：

注意 sigma 带与 alpha、beta 的边界有重叠，这是刻意为之：睡眠纺锤主频落在 12–16 Hz，用这个窗口能把它从 alpha 和 beta 里干净切出来。分期在谱域的特征很明确：W 时 alpha 与 beta 占大头；进入 N1 后 alpha 占比下降、theta 抬升；进入 N2，纺锤让 sigma 出现局部峰值；N3 时 delta 相对功率常超过 0.5，整条谱线左移。把每帧 PSD 在五个频带上积分，就得到一帧的能量特征向量，这是后面所有分类工作的输入。

2.3 相对功率与比值特征：为什么不用绝对能量

绝对功率受头皮阻抗、放大器增益、电极位置影响很大，同一受试者隔天测量都可能差一倍。工程上优先用相对功率，即每个频带功率除以总功率，乘性增益会被消掉；再配合比值特征，如 (alpha+beta)/(delta+theta)，放大阶段间的对比差异。相对功率看“比例”，比值特征看“对比”，两者对个体差异都不敏感，这也是跨受试者部署时少数能直接沿用的特征。标题里的“功率能量”落实到代码，就是 PSD 在频带上的积分：输入单位是 µV 时结果单位是 µV²，30 秒帧长固定，积分结果与“能量”成正比，除以时长就是平均功率——两者只差一个常数，做分期时可以互相替代。

2.4 功率谱估计中的窗长与分辨率矛盾

30 秒帧、100 Hz 采样得到 3000 个点。用 pwelch 做谱估计时，频率分辨率 = fs / 窗长。取 4 秒汉宁窗，分辨率是 0.25 Hz；取 2 秒窗，分辨率降到 0.5 Hz，0.5–4 Hz 的 delta 带只剩 8 个频点，积分误差偏大。窗越长谱越平滑但方差越大，折中方案是 4–6 秒窗加 50% 重叠。先用一帧试看谱形，再批量提取：

fs = 100; x = epochs(:, 1); % 只看第一帧 win = hann(4*fs); % 4 秒窗，分辨率 0.25 Hz [pxx, f] = pwelch(x, win, round(2*fs), 1024, fs); plot(f, 10*log10(pxx), 'LineWidth', 1); xlim([0 35]); xlabel('Hz'); ylabel('dB/Hz'); grid on;

pwelch 四个核心参数的含义：win 是窗函数，决定频率分辨率和谱泄漏；noverlap 取窗长一半，兼顾平滑与计算量；nfft 补零到 1024 只是细化频点显示，不改变真实分辨率；fs 必须与数据一致。看到谱形正常（0.5–35 Hz 内能量随频率递减、无 50 Hz 尖峰）再进入批量计算。

3. MATLAB 脑电预处理与频带能量提取：从原始波形到五维特征

3.1 EDF/CSV 数据读入：采样率、单位与通道对齐

第一步先看数据里有什么通道。EDF 格式用edfread读，R2020b 之后的版本返回带单位的时间表；CSV 或 Mat 导出的数据按列读向量即可。无论哪种来源，先确认采样率 fs 和单位。常见设备输出 100 Hz、128 Hz 或 256 Hz，脑电单位是 µV。采样率不统一会导致后续所有滤波器和窗参数都要重算，所以先统一到一个值：

fs = 256; % 原始采样率，从数据头读出 eeg = double(data(:, chan_idx)); % 取出 C3 或对应通道，转 double fprintf('信号长度 %.1f min, fs=%d Hz\n', length(eeg)/fs/60, fs); if fs ~= 100 % 统一降到 100 Hz eeg = resample(eeg, 100, fs); fs = 100; end

resample 两个参数的顺序是“目标采样率在前，原始采样率在后”，写反了不会报错但输出节奏全乱。降采样不只是为省计算量，更重要的是让帧长、窗长、滤波器边界都对齐同一套约定，换数据集时不用重调。单位如果不是 µV，特征里会带一个固定偏置，建议在读入时乘回增益，否则日志里所有功率值都要做一次单位换算。

3.2 滤波链：陷波、带通、零相位

脑电有效信息基本集中在 0.5–35 Hz。低于 0.5 Hz 是基线漂移和出汗伪迹，高于 35 Hz 是肌电和无线电干扰。工频干扰是 50 Hz 的窄带噪声，要先陷波再带通：陷波会引入相位畸变，带通放最后能把带外噪声再压一次。滤波器参数按下面这套配置先跑通：

% 1) 50 Hz 工频陷波 wo = 50 / (fs/2); % 归一化中心频率 bw = wo / 35; % 相对带宽 [b_n, a_n] = iirnotch(wo, bw); eeg = filtfilt(b_n, a_n, eeg); % 2) 0.5-35 Hz 带通 [b_bp, a_bp] = butter(4, [0.5 35]/(fs/2), 'bandpass'); eeg = filtfilt(b_bp, a_bp, eeg);

filtfilt 比 filter 好在零相位：滤波器系数对信号正反各过一遍，相位延迟被抵消，纺锤波和慢波不会在时间轴上平移。iirnotch 属于 DSP System Toolbox，如果没有许可证，用butter(2, [48 52]/(fs/2), 'stop')替代，带宽略宽一点，效果接近。注意带通下界不要高于 0.5 Hz，否则 N3 的核心能量会被拦在滤波器外面，delta 占比系统性偏低。

3.3 bandpower 还是 pwelch：频带能量提取的两种写法

预处理做完把信号切帧。AASM 帧长 30 秒，不重叠，末尾不足 30 秒的残段直接丢弃——整夜记录结尾往往夹着起床动作，丢掉不心疼。切帧后逐帧算五个频带的功率。MATLAB 里最省事的是 bandpower，内部做周期图估计后在指定频带积分；想看谱形做诊断就用 pwelch。两者都要，bandpower 出特征，pwelch 出谱图：

epoch_len = 30; % 帧长，AASM 标准 n = epoch_len * fs; % 每帧 3000 点 n_epoch = floor(length(eeg) / n); epochs = reshape(eeg(1:n_epoch*n), n, n_epoch); % 列向量存帧 bands = { 'delta', 0.5, 4; ... 'theta', 4, 8; ... 'alpha', 8, 13; ... 'sigma', 12, 16; ... 'beta', 13, 30 }; n_bands = size(bands, 1); feat = zeros(n_epoch, n_bands); for k = 1:n_epoch for b = 1:n_bands feat(k, b) = bandpower(epochs(:, k), fs, [bands{b,2} bands{b,3}]); end end

bandpower 的第三个参数传给它的 [下限 上限] 是闭区间，sigma 带故意从 12 Hz 到 16 Hz，覆盖纺锤主频。输出单位受输入单位影响，输入 µV 则输出 µV²。第二种写法用 pwelch 加 trapz，效果等价但多了谱图可看：

win = hann(4*fs); % 4 秒窗 [pxx, f] = pwelch(epochs(:, k), win, round(2*fs), 1024, fs); band_p = zeros(1, n_bands); for b = 1:n_bands idxb = f >= bands{b,2} & f <= bands{b,3}; band_p(b) = trapz(f(idxb), pxx(idxb)); end

提示：pwelch 输出的 pxx 单位是 µV²/Hz，trapz 对频带积分后与 bandpower 同量纲。50% 重叠是谱估计的推荐起点，nfft 取 1024 只是补零细化频点，不会提升真实分辨率，但能让频带边界附近的积分更平滑。

3.4 坏帧剔除：峰峰值阈值与缺失值处理

睡眠脑电常见三类污染：大幅运动伪迹、电极脱落、眼动低频偏移。它们的共同表现是某段波形隆起，峰峰值能到几百甚至上千 µV，而正常睡眠脑电峰峰值一般不超过 200 µV。用帧级峰峰值做质检：

pk2pk = max(epochs) - min(epochs); % 每帧峰峰值 artifact = pk2pk' > 500; % 超过阈值标记为坏帧 feat(artifact, :) = NaN; % 特征置 NaN，按缺失处理

坏帧占比低于 5% 时，分类前用相邻帧的中位数填补；超过 20% 说明记录质量本身有问题，建议整段丢弃。不要试图把坏帧功率“修正”回来——伪迹能量会铺满整个频带，任何补全都会污染 N3 的 delta 判断，得不偿失。

4. 根据功率特征判断睡眠阶段：阈值规则与 LDA 分类实现

4.1 特征向量组装：相对功率加对数功率加比值

上一章算出的 feat 有五列绝对功率，直接丢给分类器不推荐，个体基线差异会主导距离度量。组装成三类特征再拼接：相对功率、对数功率、一个综合比值。

total = sum(feat, 2); % 总功率 rel = feat ./ total; % 相对功率，n_epoch x 5 logp = log(feat); % 对数功率，压缩动态范围 ratio = log((feat(:,1)+feat(:,5)) ./ (feat(:,2)+feat(:,3))); % (delta+beta)/(theta+alpha)，越大说明越接近深睡或清醒 feat_all = [rel, logp, ratio]; % 拼接成 11 维特征

标准化用记录内统计做 zscore，而不是全局固定值：不同受试者绝对功率可能差两倍，记录内标准化能把这部分个体差异抵扣掉。ratio 的分子分母故意不含 sigma——sigma 是 N2 的核心判别信息，把它卷进综合比值会稀释 N2 的区分力。

4.2 阈值规则式分期：从 N3 开始往下判

规则法的价值是逻辑透明、可解释、好排查。判断顺序很重要：先用最有把握的 N3，再 N2，再 N1，剩下标为待定。反过来从“清醒”开始判，整晚都会被判醒。规则阈值按经验先给初值，再根据数据分布标定：

% 规则1: 深睡 N3，delta 相对功率占绝对优势 stage = zeros(n_epoch, 1); % 0=待定 stage(rel(:,1) > 0.50 & rel(:,2) < 0.30) = 3; % 规则2: N2，sigma 相对突出，且 alpha 不高 idx = stage == 0 & rel(:,4) > 0.15 & rel(:,3) <= rel(:,4); stage(idx) = 2; % 规则3: N1，theta 抬头、alpha 消失 idx = stage == 0 & rel(:,2) > 0.25 & rel(:,3) < 0.15; stage(idx) = 1; % 剩余帧保留为 0，后续用上下文或 EOG 通道再判 W/REM stage(stage == 0) = 9; % 9=未决

规则法的三个阈值各自独立：delta 判 N3、sigma 判 N2、theta 判 N1。注意规则 2 里加了一个rel(:,3) <= rel(:,4)的条件，目的是挡住清醒闭眼时慢 alpha 峰渗进 sigma 带的误判。规则顺序不能乱，一旦 N3 和 N2 的判定条件都满足，先命中的 N3 生效。

提示：如果目标是民用产品，可以按规则输出合并成“清醒 / 浅睡 / 深睡 / REM”四态。N1 与 N2 在临床也常合并为浅睡，这样还能避开单通道脑电对 N1 识别率偏低的硬伤。

4.3 LDA 兜底：fitcdiscr 训练与五折交叉验证

规则法调到一定程度就到顶了，再往前走用简单监督学习。线性判别分析是睡眠分期里性价比最高的分类器：特征维度低、样本量大（整晚约 960 帧）、类别可分性较好。公开的 Sleep-EDF 数据集自带 hypnogram 标注，可以直接提取特征和标签来验证流程：

% feat_all 与 stage_label 来自同一提取流程，stage_label 取值 0/1/2/3/4 mdl = fitcdiscr(feat_all, stage_label, 'DiscrimType', 'pseudoLinear'); cv = crossval(mdl, 'KFold', 5); acc = 1 - kfoldLoss(cv); disp(['5-fold ACC: ', num2str(acc)]);

DiscrimType 选 pseudoLinear 是关键：睡眠特征各列高度相关，相对功率列的和本身就是 1，协方差矩阵接近奇异，普通 linear 会报错，pseudoLinear 自动改用伪逆。quadratic 二次边界在小样本下容易过拟合，等数据量积累到多条记录以后再试。kfoldLoss 返回误分率，准确率写作 1 - kfoldLoss(cv)。再往上走就是深度学习路线，用 LSTM 或 CNN 直接吃谱图，但需要大量带标注记录，阈值加 LDA 作为基线更可控，也更容易定位问题出在特征还是分类器。

4.4 阈值标定的网格搜索与跨数据集迁移

阈值初值应该来自小范围网格搜索，而不是拍脑袋。用记录前 1/3 有标签的帧做标定，后面 2/3 做验证，目标函数用整体帧准确率。注意不要用 F1 做标定目标——睡眠数据里 N3 帧占比小，F1 会把少数类权重带偏，而实际部署要的是整体不出错。标定参数的经验范围：

跨数据集迁移时先看一帧分布再决定是否重标。老年受试者的 delta 占比天然偏高，直接套青壮年标定的阈值会让 N3 比例虚高；同一设备不同导联位置也会平移整体谱能量。网格搜索步长取 0.05 足够，网格范围不超过上表给出的区间，结果稳定后记录日志并存成配置文件，不要每次重新标。

5. 睡眠分期结果的验证方法与三个高频坑

5.1 用行归一化混淆矩阵验收

拿到预测序列第一步不是报总准确率，而是看混淆矩阵。行归一化后每一行表示“真实标签的帧被分到哪去了”，比单一数值信息量大得多：

confusionchart(stage_true, stage_pred, ... 'RowSummary', 'row-normalized', 'ColumnSummary', 'column-normalized');

纯脑电功率方法整体一致率做到 75–82% 就是合格线，其中 N3 召回率应高于 80%。N2 是最容易混淆的类别，很多系统把 N1 和 N2 混在一起输出浅睡，临床并不完全反对。kappa 系数剔除了随机一致的影响，0.6 以上可以接受；MATLAB 没有内置 kappa，直接用 confusionmat 输出的混淆矩阵手算即可。

5.2 坑一：肌电污染让 beta 带全线飘红

睡眠中颈部肌电不会完全消失，翻身、磨牙、呼吸都会把 20 Hz 以上的高频分量灌进脑电通道，正好压在 beta 带。症状很典型：清醒占比异常高、N2 纺锤被抬高的 beta 盖住。对策是把带通上界从 35 Hz 收紧到 30 Hz，同时把峰峰值质检阈值从 500 µV 收紧到 300 µV。如果手头有肌电通道，用肌电包络做回归剔除是更干净的办法，单通道方案里做不到，取舍要提前说清楚。

5.3 坑二：sigma 与 alpha 重叠是设计出来的，不是 bug

12–16 Hz 与 8–13 Hz 在 12–13 Hz 处有重叠。部分受试者 alpha 峰偏慢，清醒时 alpha 会渗进 sigma 带，导致清醒帧被误判成 N2。区分要点在谱形：清醒 alpha 峰尖锐幅度高，N2 纺锤宽而平。工程处理是给规则 2 加一个 alpha 相对功率不高于 sigma 的门槛，就是 4.2 节代码里的rel(:,3) <= rel(:,4)；如果数据里 alpha 峰普遍偏慢，把 sigma 边界改成 13–16 Hz 也能绕开。

5.4 坑三：N1 与 REM 在单通道上分不开，用多数投票抹平毛刺

REM 的谱特征接近“清醒化的 theta”——theta 增强、alpha 减弱，和 N1 高度重叠。临床判 REM 靠 EOG 快速眼动和 EMG 张力消失，纯脑电功率做这件事信息不够。可行的补救是利用睡眠结构：REM 通常出现在深睡循环之后且持续数分钟，孤立的单帧 REM 判断不可信。用五帧多数投票把毛刺抹掉：

m = 5; h = (m-1)/2; s = zeros(size(stage)); for k = 1:length(stage) lo = max(1, k-h); hi = min(length(stage), k+h); s(k) = mode(stage(lo:hi)); % 窗口内众数投票 end

m=5 不会出现平票，代价是丢失 5 帧量级的真实转换边界，但睡眠分期本身以 30 秒为最小单位，这个代价可接受。手动复核时把结果叠在频谱图上：pspectrum(eeg_epoch, fs, 'spectrogram')，一眼能看出阈值是不是在乱分；训练好的模型记得把决策边界和特征分布一起存成 MAT 文件，下次换数据先跑同一条验证脚本再谈迁移。

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