简介面向脑电信号分析的多尺度熵计算MATLAB代码包主要服务神经科学、生物医学工程等领域的科研人员与工程师可用于处理多通道EEG数据完成从二进制文件读取、预处理、多尺度熵计算到频谱分析、批量统计的完整流程。代码由杨博士根据C程序改写并附带C语言辅助工具适合需要自定义参数进行脑电复杂度分析的进阶用户。压缩包共40个文件大小仅23KB以26个m脚本为主体覆盖数据加载、粗粒化、熵计算、批处理及可视化等模块另有2个c源文件、电极位置文件、文本说明及README文档便于理解程序结构和电极对应关系。该资源已有342人学习属于轻量级工具包。使用者可直接运行主脚本或借助批量平均熵计算函数实现多文件批量处理输出各时间片段、通道与尺度因子的熵值同时附带结果核对与绘图辅助模块便于双重检查。整体目录组织清晰适合作为脑电多尺度熵分析的参考实现或二次开发基础。1. 脑电功率谱与多尺度熵分析这份 MATLAB 工具箱能同时啃下频域特征和质谱数据预处理先抛一个可能颠覆认知的结论标题里的「mse-analysis」不是质谱Mass Spectrometry专用库而是一个以 MATLAB 为宿主、主攻脑电信号频域特征与非线性熵特征的代码包只是顺便把质谱数据读取和基线校正也做进去了。我拆这类科研工具箱的习惯是先看它的 IO 边界输入是脑电通道矩阵输出是功率谱密度曲线和多尺度熵曲线中间穿插滤波、去伪迹、批处理脚本。对搞脑机接口、睡眠分期、情绪识别或者认知负荷评估的人来说这套代码能省掉大量「从零搭分析流水线」的时间。目标人群很明确手里有 EEG 采集设备、但 MATLAB 信号处理工具箱用得不够溜的从业者以及需要把「功率谱 多尺度熵」做成联合特征喂给分类器的研究者。它不是点开即用的软件而是需要你自己按数据格式微调的半成品框架——这恰恰是它最有价值的地方。2. 功率谱计算选 Welch 而不是周期图参数设置全在脚本里2.1 为什么 Welch 法是脑电频段分析的默认解脑电信号是非平稳、低信噪比的随机过程直接对整个数据段做傅里叶变换周期图法得到的频谱方差极大单次估计根本不能代表真实频段能量分布。Welch 法把长信号切成重叠段、逐段加窗、计算周期图后取平均用「牺牲频率分辨率换取方差降低」的思路换稳定结果。这套代码里默认窗长是采样率的 2 倍——比如 500 Hz 采样率就取 1000 点窗频率分辨率大约 0.5 Hz足够区分 Delta0.5–4 Hz、Theta4–8 Hz、Alpha8–13 Hz、Beta13–30 Hz这些经典频段。如果你要分析慢皮层电位或者极低频成分需要把窗长调到采样率的 4–8 倍但代价是平均段数变少方差回升。这套代码把window和noverlap定义在脚本头部改起来比在 GUI 里点选直观得多。2.2 从原始通道矩阵到频段能量占比的可复现步骤假设你已经有一个channels × samples的二维矩阵eeg_data采样率fs已知下面这段代码是这个工具箱的核心起点% 核心参数定义 fs 500; % 采样率单位 Hz window fs * 2; % 窗长 2 倍采样率频率分辨率 0.5 Hz noverlap window * 0.5; % 重叠率 50%常用折中值 nfft 2^nextpow2(window); % FFT 点数取 2 的幂加速计算 % 逐通道计算功率谱密度 for ch 1:size(eeg_data, 1) [pxx, f] pwelch(eeg_data(ch, :), window, noverlap, nfft, fs); psd_all(ch, :) pxx; % 行 通道列 频率点 end % 定义经典频段边界单位 Hz bands [0.5 4; 4 8; 8 13; 13 30; 30 45]; band_names {Delta, Theta, Alpha, Beta, Gamma}; % 计算每个通道各频段的相对能量占比 for ch 1:size(psd_all, 1) total_energy sum(psd_all(ch, :)); for b 1:size(bands, 1) idx f bands(b, 1) f bands(b, 2); band_ratio(ch, b) sum(psd_all(ch, idx)) / total_energy; end end这段代码的逻辑链是三层递进的第一层pwelch做的是单通道功率谱密度估计返回频率向量f和对应功率值pxx第二层把所有通道的pxx堆叠成矩阵方便后续批量处理第三层把连续的频谱按经典频带分段求和并归一化得到「相对能量占比」这种分类器友好的特征。需要留意nfft的取法2^nextpow2(window)会自动把 FFT 点数对齐到大于窗长的最小二次幂实际计算频率分辨率由采样率和nfft共同决定而不是由窗长单独决定。如果你发现频谱曲线毛刺很多优先把noverlap提高到 75%代价是计算时间增加如果发现峰被展宽、频段边界模糊说明窗太长频率分辨率过剩可以缩短窗长。2.3 功率谱结果的可视化与导出陷阱工具箱里附带了一个简单绘图脚本本质上就是plot(f, 10*log10(psd_all(ch,:)))配合xlim([0 50])但有一个值得注意的工程细节不同受试者头部电阻不同、采集设备增益不同绝对功率值个体差异大到离谱直接用plot看原始 PSD 曲线几乎得不出横向结论。这套代码内置了对数转换和归一化选项默认输出的是10*log10(pxx / sum(pxx))这种相对值。我在实际项目里会把每个频段的绝对功率、相对功率、峰值频率、谱熵这四个量全部落盘因为投论文时审稿人经常问「你们的 Alpha 功率是怎么定义的」留着原始指标才能应对追加分析。还有导出环节的坑dlmwrite写 CSV 是旧习惯新项目建议直接用writematrix配合cell2table把频率点和通道名带上避免下游 Python 脚本读进来还得重新对位。3. 多尺度熵与质谱数据预处理同一个工具包里的两条特征主线3.1 多尺度熵背后的粗粒化逻辑和参数选择多尺度熵Multiscale Entropy, MSE是衡量信号在多个时间尺度上复杂度的指标核心思想是先对原始信号做粗粒化——把每τ个点平均成一个新点——再计算每个尺度的样本熵。脑电信号在粗粒化后如果熵值能保持稳定说明该频段存在长程相关性如果熵值快速下降说明信号偏向随机噪声。这套代码里的mse_analysis.m接受三个关键参数嵌入维数m默认 2、相似容差r默认 0.15 倍标准差、最大尺度因子max_tau默认 20。m取 1 对短数据友好但容易漏掉非线性结构取 3 以上需要的数据量按指数增长r取 0.15 是样本熵领域的经验公约数取太小噪声会被当有效信号取太大所有序列都显得相似。运行时间随max_tau线性增长20 个尺度在 500Hz 采样率、30 秒脑电数据上大概要跑半分钟到一分钟做好心理准备。3.2 一组形态正确的模拟验证代码直接拿真实脑电跑 MSE 有个问题——你不知道「正确答案」长什么样。我一般先用合成信号验证这段代码再切真实数据% 生成模拟信号白噪声 正弦波 1/f 噪声混合 rng(42); % 固定随机种子 t 0:1/fs:30; % 30 秒数据 white_noise randn(size(t)); sine_wave sin(2*pi*10*t); % 10 Hz 正弦成分 alpha_osc sin(2*pi*10*t) .* (0.5 0.5*sin(2*pi*0.5*t)); % 幅度调制振荡 sig 0.5*white_noise sine_wave alpha_osc; % 提取连续 10 秒片段分析 seg sig(1:fs*10); mse_curve zeros(1, 20); for tau 1:20 % 粗粒化每个非重叠窗口取均值 seg_reshape reshape(seg(1:floor(length(seg)/tau)*tau), [], tau); coarse mean(seg_reshape, 2); % 调用工具箱中的样本熵函数 mse_curve(tau) sample_entropy(coarse, 2, 0.15); end % 绘制熵曲线 figure; plot(1:20, mse_curve, LineWidth, 2); xlabel(Scale Factor); ylabel(Sample Entropy);这段代码至少验证了三件事固定随机种子保证结果可复现合成的三种成分白噪声、固定频率正弦、调幅振荡对应真实脑电中的噪声、稳态节律、事件相关同步化单尺度熵算完后再画曲线看趋势。运行后你会看到白噪声的样本熵初始值很高、随后几乎水平直线正弦波的样本熵在低尺度就跳到接近零——因为周期信号在粗粒化后变得完全可预测调幅振荡的熵曲线则会在特定尺度出现峰值。如果你的sample_entropy函数算出来纯白噪声的熵值低于 1.5说明r的取值有问题正常白噪声在m2, r0.15条件下样本熵应该在 1.5–2.2 之间。这一步验证通过之前不要拿真实脑电数据跑任何下游统计。3.3 附带质谱工具的真实定位与适用场景标题中「质谱分析」这部分工具箱里实际实现的是质谱数据格式读取与基线校正不是峰识别更不是化合物鉴定。它包含一个read_mzxml.m函数能把 mzXML 格式的质谱文件读成(m/z, intensity)两列矩阵附带一个基于局部最小值拟合的基线扣除脚本。这个模块适用的场景是质谱仪导出的原始数据需要先梳理出干净谱线、再做后续的峰对齐或 PCA 之类多元分析。注意它不包含数据库检索功能也不擅长处理高分辨质谱的谱图去卷积。我在某实验室见过他们把脑电和质谱数据放在同一个工具包里管理分析流程是「质谱找代谢标志物 → 脑电找神经反应特征 → 两类特征做相关性分析」这套代码正好卡在这条流程的上游位置。基线校正的强度参数写在脚本第三行默认是 0.15 倍峰高如果你处理的样品基线漂移很厉害把这个值提到 0.3代价是弱峰可能被一并扣掉。4. 避坑与排查五个高频翻车现场及修复办法我从拆包到实际跑通这套代码前后踩了不少坑主要集中在五个方向上。把坑写出来比贴十页源码更能说明问题。坑位一pwelch 报错「Window length must be less than or equal to signal length」现象代码跑到pwelch直接红字报错检查不出逻辑问题。 原因这一步几乎都是因为数据短于窗长。如果采样率 1000 Hz、数据只有 1 秒默认窗长 2000 点直接超过信号长度。脑电采集时导联脱落、实验提前终止都可能产生这类短段。 解决在调用pwelch前加一行保护性判断if size(eeg_data, 2) window, eeg_data [eeg_data, eeg_data]; end或者直接把窗长改成min(fs*2, size(eeg_data, 2)/2)。工具箱原版没有这个保护逻辑你需要自己加。坑位二多尺度熵的计算时间长得像死机一样现象mse_analysis跑 20 个尺度在旧电脑上整整五分钟没出结果CPU 占用却很低。 原因原始代码用的嵌套 for 循环计算样本熵复杂度是 O(N²) 级别数据长度 10 秒 × 500 Hz 5000 点逐点算距离矩阵是内存密集型操作。 解决把max_tau从 20 降到 10低频成分仍然能覆盖或者把样本熵函数里计算距离矩阵那部分改成向量化操作用pdist2配合sum替换内层 for 循环速度能提升 5–10 倍。如果数据超长还可以先用resample把采样率降到 250 Hz 再跑。坑位三所有通道的功率谱曲线一模一样现象画图时看到全部通道的频谱轮廓完全相同以为是算错了。 原因我没注意到代码里有个隐藏的mean(eeg_data, 1)操作——它把多通道信号压成了一个平均导联再算 PSD。这对某些参考电极方案可能合理但如果你想分析的是 C3/C4/Cz 这些特定通道的空间差异这个预处理就是毁灭性的。 解决检查脚本头部是否有avg_ref mean(eeg_data, 1)之类的赋值有就注释掉然后把pwelch里的输入参数改为eeg_data(ch, :)。坑位四质谱基线校正后强度出现负值现象baseline_correct跑完谱图里的强度列有接近一半是负数下游分析没法继续。 原因局部最小值拟合的基线被噪声拉到了真实零线以下扣除时把统计噪声也当成信号减掉了。 解决在基线校正前先对强度列做一次 7 点滑窗平滑或者在脚本里给基线加一个最小阈值baseline_thresh max(baseline) * 0.05让低于阈值的量直接截断为零。这一步在原代码里没有属于质谱分析的标准工程处理。坑位五批处理脚本只处理了第一段数据就停了现象循环遍历 50 个文件只输出了一个结果文件日志显示处理完第一个就报错跳出。 原因文件读取函数里有[~, name, ext] fileparts(filepath)这行代码如果文件名包含中文或路径带空格解析结果异常循环退出时抛错被catch吞掉了日志。 解决把所有输入文件名改成英文 下划线不要用空格同时把循环包在try...catch...end里但把错误信息写入日志文件这样哪一行依赖的配置文件错了都能查。5. 把特征批量落盘并和 Python 下游对接一个顺手多年的收尾习惯代码跑通不是终点真正痛苦的是分析完以后的特征矩阵怎么交出去。脑电领域几乎必然要和 Python 的 scikit-learn 或 PyTorch 打交道所以我在这套工具箱里养成了一个固定习惯用 MATLAB 把特征全部整理成一份长表 CSV再用 Python 读回做模型训练。具体做法是设计一个统一出口行是样本片段列是「被试编号、通道名、频段名、尺度序号、特征值」五个字段。这样后续做混合效应模型、画箱线图、做分类器输入都不用回 MATLAB 再导一次。用代码描述就是% 整理特征表的结构 FeatureTable table(); for subj 1:n_subjects for ch 1:size(psd_all, 1) for b 1:length(band_names) FeatureTable [FeatureTable; table(subj, ch, band_names{b}, band_ratio(ch,b), ... VariableNames, {Subject,Channel,Band,Value})]; end end end writetable(FeatureTable, features.csv);Python 端读回来就一行pd.read_csv(features.csv)然后groupby([Subject,Band])想怎么透视怎么透视。这个做法有效的前提是在 MATLAB 端就把通道名、频段定义、采样率这些元信息都打进列名里而不是只丢一个矩阵让下游猜。如果你要走得更深可以把 MSE 的 20 个尺度值也拉成一行宽表——行是样本、列是MSE_Scale1到MSE_Scale20——这样标准化之后直接进分类器符合大多数机器学习框架的输入习惯。我现在的项目流程固定是「跑 Welch 出频段特征 → 跑 MSE 出熵特征 → MATLAB 拼宽表 → Python 训练模型」三十分钟一个迭代效率比之前两边轮流改脚本高出不少。这套工具包把前两步封装好了你只需要在出口上对齐格式。关于质谱数据这套我后来也习惯把基线校正前后的强度值同时导出因为复现实验时经常需要判断校正强度是否合理留原始值等于给自己留了后悔药。希望这个收尾习惯帮到你至少能省掉一次「模型结果不对但不知道特征是哪个环节算错」的排查经历。本文还有配套的精品资源点击获取