简介这套MATLAB代码包面向脑电信号处理与生物医学工程研究者提供基于多尺度熵与功率谱的分析工具对应mse-analysis开源项目。代码由杨博士将Costa的C语言程序重写为MATLAB版本支持从二进制文件导入20通道、128Hz采样的脑电数据并完成粗粒度化、熵值计算、批量平均与频谱分析等流程。压缩包共40个文件体积仅23KB以26个m脚本为核心辅以7个txt数据说明、2个c源文件及README文档覆盖数据加载、熵值计算、批量处理和结果绘图等模块便于按需修改和二次开发。借助批处理与检查脚本可一键处理多个文件输出各通道、各尺度因子的熵值与频谱图附带的说明文档有助于理清依赖关系和实验设计快速搭建或复现脑电分析流程。已有342人浏览学习适合需要从零实现多尺度熵计算的初学者及从事脑电分析的中高级研究人员。 写这篇东西的起因挺有意思——我最近在整理自己的MATLAB工具包时发现一个命名叫“mse-analysis”的文件夹里面既有脑电信号处理脚本又有几段质谱数据读取的测试代码。这个文件夹的混乱命名一度让我自己都困惑它到底是做脑电功率谱的还是做质谱分析的后来仔细捋了一遍才发现这里的“mse”在脑电语境下应该被理解成Multiscale Entropy多尺度熵而“质谱分析”则是另一个完全无关的项目分支被错误归档了。这种命名歧义在实际工程项目里太常见了但恰恰是这种混乱促使我把整个脑电功率谱分析管线彻底理清了一遍。如果你也正在用MATLAB处理脑电数据、想算功率谱密度或者多尺度熵又或者你手头有一个同样命名不清的代码包需要整理那这篇东西应该能帮你少走不少弯路。1. 先把标题里的“质谱分析”四个字掰扯清楚在展开任何代码之前我建议你先停下来想清楚一件事你手里的这个“mse-analysis”到底是干什么的。因为“MSE”这个缩写本身就有多重含义在不同领域代表完全不同的算法。1.1 脑电领域的MSE多尺度熵在脑电信号处理中MSE全称是Multiscale Entropy多尺度熵。它不是直接算功率谱的而是用来衡量信号在不同时间尺度上的复杂度或自相似性。简单理解脑电信号的熵值高意味着信号更“无序”、信息量更大熵值低则意味着信号更规则、更单调。多尺度熵的核心思路是先把原始信号在不同尺度下做粗粒化coarse graining然后在每个尺度上计算样本熵Sample Entropy最后看熵值随尺度的变化曲线。这个指标在很多认知任务研究、睡眠分期、麻醉深度监测里都有应用。它和功率谱是互补的关系——功率谱告诉你“各个频率上有多少能量”多尺度熵告诉你“信号在不同时间尺度上有多复杂”。1.2 质谱领域的MS-E质量误差分析而在质谱分析Mass Spectrometry领域MS-E通常指代的是Mass Error质量误差分析。这是蛋白质组学、代谢组学里非常常规的一步操作——鉴定出的肽段或代谢物的理论质量数和实测质量数之间往往存在微小偏差这个偏差的分布特征直接关系到仪器校准状态和鉴定结果的置信度。质谱分析里的mse-analysis代码包做的应该是峰检测、质量校准、误差分布统计这些事情。1.3 命名混淆的根源在哪现在的尴尬就来了一个叫“mse-analysis”的仓库既可以装脑电多尺度熵分析代码也可以装质谱质量误差分析代码甚至可以同时装脑电功率谱和质谱分析两套不相干的东西。我在整理代码时发现真正导致混淆的往往是当初下载或创建仓库时顺手打了“mse-analysis”这个名字没有任何上下文注释。如果你是从某个学术分享链接拿到的这个包打开一看既有“power spectrum”又有“mass error”那大概率是仓库作者把多项目脚本混在一个目录下了。我的建议是先打开代码包里的README或者主脚本看它导入的数据格式。如果读入的是.edf、.set、.mat格式的多通道时序数据那这就是脑电分析包如果读入的是.mzML、.raw、.csv格式的质谱扫描数据那这就是质谱分析包。如果两类都有那你就得自己在目录层面拆分归档别指望跑通全部代码。2. 脑电功率谱分析最容易被忽略的前置工作数据与预处理明确了“mse-analysis”在脑电语境下的真实定位是功率谱加多尺度熵之后接下来就是最枯燥但最关键的一步——数据准备。我见过太多人在这一步栽跟头拿到代码包迫不及待地导入一段脑电数据就开始跑功率谱结果出来的频谱图全是毛刺和趋势项完全没法看。问题几乎都出在预处理环节。2.1 先确认你的数据长什么样脑电数据通常是一个通道数 × 采样点数的矩阵或者是采样点数 × 通道数。如果你用的公开数据集比如DEAP、SEED、BCI Competition读取后第一件事不是算功率谱而是确认采样率和通道布局。采样率直接决定你能分析的频率上限——根据奈奎斯特定理能分析的最高频率是采样率的一半。比如采样率是256Hz那你能看到的频谱范围就是0到128Hz而脑电研究的核心频段delta、theta、alpha、beta、gamma基本都在0.5到50Hz之间所以128Hz以上的信息对常规分析意义不大。2.2 预处理管线去趋势、滤波、分段缺一不可很多新手拿到的“mse-analysis”代码包预处理部分往往很简单就是直接对原始信号做FFT。但实际可用的流程应该包含以下几步去除基线漂移和线性趋势脑电采集过程中电极与皮肤之间的极化电压、受试者出汗、设备温漂都会让信号带上一个缓慢变化的趋势项。这个趋势项在频谱上表现为极低频段的巨大能量会压低其他频段的相对功率。用detrend函数做去趋势处理是最基本的操作。带通滤波常规脑电分析的频率范围取0.5Hz到50Hz或0.5Hz到100Hz。低于0.5Hz的成分主要是漂移和伪迹高于50Hz或100Hz的成分多是肌电噪声。这个滤波的作用不是让频谱“好看”而是避免无关频段能量干扰后续的频段划分。MATLAB里直接bandpass(eegData, [0.5 50], fs)就能完成零相移滤波。剔除坏导联和坏段如果某个通道的波形明显是平线、剧烈毛刺或者大范围漂移这个通道的数据就不该参与功率谱计算。通常做法是设定一个幅度阈值或者方差阈值超过阈值的通道直接置空或插值替换。2.3 分段策略直接影响频谱分辨率功率谱的频率分辨率由参与计算的信号长度决定。计算公式很简单频率分辨率 采样率 / 参与FFT的点数。如果你用1秒的数据做FFT256Hz采样率下分辨率就是1Hz用4秒的数据分辨率可以到0.25Hz。而脑电的alpha节律8-13Hz和theta节律4-8Hz分界线很近频率分辨率太低的时候频谱图上频段边界会糊在一起很难精确算各频段能量。所以一套合理的处理流程是把连续脑电按2秒或4秒一段切分每次分析的时间窗口叫epoch对每个epoch分别计算功率谱再把所有epoch的功率谱做平均。这其实就是Welch方法的思路加窗分段、分别FFT、再平均用方差换稳定。MATLAB里pwelch函数就是这么做的[pxx, f] pwelch(eegData, window, noverlap, nfft, fs)一行调用就能得到平滑的功率谱密度估计。3. 从原始EEG到可解释的频段功率谱完整代码管线预处理做完以后真正的核心计算才开始。这一段我会给出一个相对完整、可以直接照抄改写的MATLAB实现。它不是从某个现成包里拿来的而是我根据自己的项目经验整理的一套清晰流程。3.1 核心计算代码Welch功率谱估计与频段功率汇总function powerTable compute_band_power(eegData, fs, channels) % 输入 % eegData - 通道数x采样点数的矩阵 % fs - 采样率单位Hz % channels- 通道名称元胞数组 % 输出 % powerTable - 各通道在各频段的绝对功率与相对功率表 % 定义频段边界单位Hz bands struct(... delta, [0.5 4], ... theta, [4 8], ... alpha, [8 13], ... beta, [13 30], ... gamma, [30 50]); bandNames fieldnames(bands); nBands numel(bandNames); % 去掉每个通道的线性趋势 eegDetrended detrend(eegData); % detrend默认按列操作所以先转置 % 零相移带通滤波 0.5-50Hz eegFiltered bandpass(eegDetrended, [0.5 50], fs); % 分段参数窗长2秒50%重叠 windowLen 2 * fs; noverlap round(0.5 * windowLen); nfft 2^nextpow2(windowLen); % 保证FFT点数足够频率分辨率约0.5Hz nChannels size(eegFiltered, 1); powerTable table(); for ch 1:nChannels % 对单通道数据做Welch功率谱估计 [pxx, f] pwelch(eegFiltered(ch, :), hann(windowLen), noverlap, nfft, fs); % 计算各频段绝对功率对功率谱密度积分 absPower zeros(1, nBands); for b 1:nBands idx f bands.(bandNames{b})(1) f bands.(bandNames{b})(2); % 用梯形积分近似频段能量 absPower(b) trapz(f(idx), pxx(idx)); end % 相对功率各频段占总功率的百分比 totalPower sum(absPower); relPower absPower / totalPower * 100; % 存入表格 row table(); row.Channel {channels{ch}}; for b 1:nBands row.([Abs_ upper(bandNames{b})]) absPower(b); row.([Rel_ upper(bandNames{b})]) relPower(b); end powerTable [powerTable; row]; end end这段代码表面上不长但有几个值得说的设计考虑为什么要用Welch而不是直接FFT直接对整段数据做FFT得到的是整段信号的频率成分平均一旦数据里混入眨眼伪迹或运动伪迹频谱会被局部突发能量污染。Welch方法把数据切成长度较短的段计算再做平均相当于把偶发伪迹的影响“稀释”掉。代价是频率分辨率会降低但脑电频段划分本身就有一定的宽容度0.5Hz的分辨率完全够用。为什么用梯形积分而不是直接求和pwelch返回的是功率谱密度PSD单位是μV²/Hz某个频段的能量应该是这个频段范围内PSD的积分。直接求和是等间隔采样下积分的近似但用trapz更精确尤其是频段边界落在FFT频率点之间时梯形积分能更好地逼近真实面积。为什么窗函数选hannFFT默认相当于加了矩形窗矩形窗的频谱旁瓣衰减只有约13dB会让强频率成分的能量“泄漏”到邻近频段。汉宁窗hann旁瓣衰减提高到约31dB能显著减少谱泄漏。hamming窗也有类似效果但hann窗在主瓣宽度和旁瓣衰减之间更均衡是脑电功率谱计算中最常用的选择。3.2 从功率值到可解读的结论算完各频段绝对功率和相对功率后直接面对的问题就是怎么解读。绝对功率受个体差异和电极阻抗影响比较大同样的alpha节律不同受试者头皮厚度不同功率绝对值可能差好几倍。所以组间比较时通常看相对功率即各频段占总功率的百分比这相当于做了个简单的归一化。表格式的输出会非常直观通道绝对δ绝对θ绝对α绝对β相对δ(%)相对θ(%)相对α(%)相对β(%)Fz12.38.725.15.224.117.049.110.2Cz15.69.218.36.831.218.436.613.6按我的经验做认知实验时重点关注alpha和theta的相对功率变化做睡眠分析时delta和theta是核心频段做癫痫检测时则要细化到次频段比如alpha1是8-10Hzalpha2是10-13Hz。千万别拿一套固定的频段划分应付所有场景。3.3 一个实测中的意外基线段的选取第一次跑这套流程时我被一个看似简单的问题卡了很久——基线段的选取。做静息态脑电分析时基线通常是闭上眼睛、放松状态下的记录段用来和任务态做对比。但不同被试的闭眼静息态alpha功率差异巨大哪怕用了相对功率组间方差还是很大。后来我改用“任务前1分钟的静息段作为基线”并且对每个通道单独做基线归一化效果才稳定下来。这个细节在现成的“mse-analysis”代码包中往往没有得自己在实验设计阶段就定好。4. 多尺度熵MSE功率谱之外的复杂度视角如果你手头的代码包名字是“mse-analysis”那它大概率不只包含功率谱计算还会包含多尺度熵算法。既然标题里带了这个词这一节就把多尺度熵的原理和MATLAB实现完整讲透。4.1 多尺度熵到底在算什么多尺度熵的基本流程分两步粗粒化和样本熵计算。**粗粒化coarse graining**是把长度为N的信号以尺度因子τ为窗长依次取平均。尺度τ1时就是原始信号τ2时把相邻两点平均成一个点信号长度减半τ3时三点平均成一点以此类推。这一步的物理含义是把原始信号在不同时间分辨率上重新“采样”尺度越大看到的时间窗口越宽高频细节被平滑掉留下的是低频趋势性变化。**样本熵Sample Entropy**衡量的是信号中“新模式出现”的概率。具体来说设定嵌入维度m和相似容差r统计信号中长度为m的子序列和长度为m1的子序列中彼此距离小于r的比例再取负对数。样本熵越大说明信号中产生新模式的概率越高信号越复杂。把每个尺度上的样本熵连成一条曲线就是多尺度熵曲线。健康成年人在脑电上通常表现为尺度1到尺度5区间熵值较高然后随尺度增大缓慢下降而病理状态或衰老状态下这个“复杂性资源”会减少熵值整体下沉或快速衰减。4.2 MATLAB实现多尺度熵的注意点多尺度熵的MATLAB实现网上能找到很多版本但真正用起来有几个细节决定结果可靠性参数选择直接影响结论嵌入维度m一般取2容差r取0.1到0.25倍原始信号标准差。r取太大样本熵会趋近于0区分度丢失r取太小样本熵对噪声过于敏感统计波动很大。我通常取r 0.15 * std(signal)在公开数据集上效果比较稳定。信号长度要够样本熵对数据长度很敏感数据太短统计值不可靠。经验法则原始数据的长度至少是最大尺度τ的10的m1次方量级以上。换句话说如果你要算到尺度10、m2那原始数据至少需要约10*10²1000个点。在脑电分析里如果采样率是250Hz4秒的数据才1000个点勉强够用最好能用到20秒以上的连续数据。function mseCurve compute_mse(signal, maxScale, m, r) % 多尺度熵计算 % 输入 % signal - 1xN的输入信号 % maxScale - 最大尺度因子 % m - 嵌入维度通常取2 % r - 相似容差通常取0.15*signal标准差 % 输出 % mseCurve - 1xmaxScale的多尺度熵值 N length(signal); mseCurve zeros(1, maxScale); for tau 1:maxScale % 粗粒化 coarseLen floor(N / tau); coarseSignal zeros(1, coarseLen); for i 1:coarseLen segment signal((i-1)*tau 1 : i*tau); coarseSignal(i) mean(segment); end % 计算样本熵 mseCurve(tau) sample_entropy(coarseSignal, m, r); end end function se sample_entropy(signal, m, r) N length(signal); if N 10^(m1) error(信号长度不足样本熵结果不可靠); end % 构造m维和m1维模板向量 templates_m zeros(N-m, m); templates_m1 zeros(N-m, m1); for i 1:N-m templates_m(i, :) signal(i : im-1); templates_m1(i, :) signal(i : im); end % 计算匹配对数用距离小于r作为匹配条件 % 这里采用优化的归一化方法避免双重循环过慢 A 0; % 匹配长度为m1的对数 B 0; % 匹配长度为m的对数 % 用pdist2算距离矩阵注意排除自身匹配 distM pdist2(templates_m, templates_m); distM(1:N-m1:end) inf; B sum(distM(:) r) / 2; % 距离矩阵对称除以2去除重复计数 distM1 pdist2(templates_m1, templates_m1); distM1(1:N-m:end) inf; A sum(distM1(:) r) / 2; if A 0 || B 0 se NaN; % 无匹配时返回NaN避免log(0) else se -log(A / B); end end这个实现里用的pdist2是我后来优化过的方式。网上很多版本的样本熵实现用三层for循环数据长度到几千点的时候跑一个尺度序列要等半天。用距离矩阵代替代数循环计算复杂度从O(N²×m)降到O(N²m×N²)速度提升非常明显。我项目里有一段5万点的脑电数据旧代码跑完20个尺度用了将近两分钟优化后5秒多就出结果了。4.3 功率谱和多尺度熵怎么配合使用功率谱和多尺度熵并不是非此即彼的关系它们各自捕捉信号的不同特征。在实际项目中如果只做静息态差异分析功率谱的相对功率指标通常就够用了但如果你关注的是认知负荷变化、意识水平变化这类非线性动态特征多尺度熵就有不可替代的优势。举个例子我在处理一组持续注意力任务的数据时发现随着任务时间延长被试的alpha相对功率明显上升——信号变得“更同步、更放松”与此同时多尺度熵在尺度2到尺度4上显著下降——信号变得“更规则、复杂度降低”。两者反映的其实是同一个神经状态变化的不同侧面功率谱看到的是频域能量重新分配多尺度熵看到的是时间序列可预测性增加。把这两种指标联合起来建模比单独用任何一个都能更好地预测行为表现。5. 如果真的遇到质谱分析代码怎么识别和离散现在回到标题里的“质谱分析”。前面说了在脑电语境下“mse”理解成多尺度熵更合理但如果你拿到的代码包里真的有质谱相关的脚本又该怎么处理先确认质谱分析的代码特征。质谱数据的基本格式是一个二维数组扫瞄时间或保留时间一维质荷比m/z一维强度值填充在矩阵里。识别质谱代码的典型特征是脚本里会出现mz、intensity、peak picking、mass calibration这类变量名或函数名。MATLAB里做质谱分析一般会用到mspeaks、msalign这些生物信息学工具箱的函数。如果你的“mse-analysis”包里出现这些函数调用那它确实包含质谱分析模块。遇到这种情况我的经验是立即把仓库按功能拆分。脑电功率谱和质谱分析虽然都叫“信号处理”但数据来源、预处理方式、评价指标完全不同硬放在一个目录里只会让后续维护越来越痛苦。拆分的依据不是代码风格而是I/O边界——看每个脚本读入什么格式的文件输出什么类型的文件凡是读入输出互不依赖的脚本就放进独立的子目录。至于质谱分析本身如果确实需要做MATLAB里最常跑通的任务是峰检测和误差分布统计。mspeaks可以检测质谱峰msalign做质量轴校准然后用拟合误差分布的方式评估质量准确度。这部分代码本身并不复杂但它和脑电分析没有半点交集硬要在一个项目里混合只会把README写得越来越长、越来越让人看不懂。写在最后的一个小经验这段时间整理代码最大的收获是不要相信文件夹名字也不要相信下载时随手起的项目名你的代码包里很可能同时躺着“脑电功率谱”、“多尺度熵”、“质谱分析”三种截然不同的东西而它们共享了一个叫“mse-analysis”的帽子。面对这种情况先花十分钟把数据格式和函数依赖理清楚比急着跑出第一张频谱图重要得多。至于脑电功率谱和多尺度熵这两部分按文中的代码管线走一遍再根据你自己的采样率和频段需求调整参数通常就能稳定出结果了。如果你只是为了一时方便把命名搞得更乱那我可以保证三个月后的你一定会回来感谢那个现在愿意花时间拆分的自己。本文还有配套的精品资源点击获取 用MATLAB理清脑电功率谱与多尺度熵分析：从MSE命名歧义说起
简介这套MATLAB代码包面向脑电信号处理与生物医学工程研究者提供基于多尺度熵与功率谱的分析工具对应mse-analysis开源项目。代码由杨博士将Costa的C语言程序重写为MATLAB版本支持从二进制文件导入20通道、128Hz采样的脑电数据并完成粗粒度化、熵值计算、批量平均与频谱分析等流程。压缩包共40个文件体积仅23KB以26个m脚本为核心辅以7个txt数据说明、2个c源文件及README文档覆盖数据加载、熵值计算、批量处理和结果绘图等模块便于按需修改和二次开发。借助批处理与检查脚本可一键处理多个文件输出各通道、各尺度因子的熵值与频谱图附带的说明文档有助于理清依赖关系和实验设计快速搭建或复现脑电分析流程。已有342人浏览学习适合需要从零实现多尺度熵计算的初学者及从事脑电分析的中高级研究人员。 写这篇东西的起因挺有意思——我最近在整理自己的MATLAB工具包时发现一个命名叫“mse-analysis”的文件夹里面既有脑电信号处理脚本又有几段质谱数据读取的测试代码。这个文件夹的混乱命名一度让我自己都困惑它到底是做脑电功率谱的还是做质谱分析的后来仔细捋了一遍才发现这里的“mse”在脑电语境下应该被理解成Multiscale Entropy多尺度熵而“质谱分析”则是另一个完全无关的项目分支被错误归档了。这种命名歧义在实际工程项目里太常见了但恰恰是这种混乱促使我把整个脑电功率谱分析管线彻底理清了一遍。如果你也正在用MATLAB处理脑电数据、想算功率谱密度或者多尺度熵又或者你手头有一个同样命名不清的代码包需要整理那这篇东西应该能帮你少走不少弯路。1. 先把标题里的“质谱分析”四个字掰扯清楚在展开任何代码之前我建议你先停下来想清楚一件事你手里的这个“mse-analysis”到底是干什么的。因为“MSE”这个缩写本身就有多重含义在不同领域代表完全不同的算法。1.1 脑电领域的MSE多尺度熵在脑电信号处理中MSE全称是Multiscale Entropy多尺度熵。它不是直接算功率谱的而是用来衡量信号在不同时间尺度上的复杂度或自相似性。简单理解脑电信号的熵值高意味着信号更“无序”、信息量更大熵值低则意味着信号更规则、更单调。多尺度熵的核心思路是先把原始信号在不同尺度下做粗粒化coarse graining然后在每个尺度上计算样本熵Sample Entropy最后看熵值随尺度的变化曲线。这个指标在很多认知任务研究、睡眠分期、麻醉深度监测里都有应用。它和功率谱是互补的关系——功率谱告诉你“各个频率上有多少能量”多尺度熵告诉你“信号在不同时间尺度上有多复杂”。1.2 质谱领域的MS-E质量误差分析而在质谱分析Mass Spectrometry领域MS-E通常指代的是Mass Error质量误差分析。这是蛋白质组学、代谢组学里非常常规的一步操作——鉴定出的肽段或代谢物的理论质量数和实测质量数之间往往存在微小偏差这个偏差的分布特征直接关系到仪器校准状态和鉴定结果的置信度。质谱分析里的mse-analysis代码包做的应该是峰检测、质量校准、误差分布统计这些事情。1.3 命名混淆的根源在哪现在的尴尬就来了一个叫“mse-analysis”的仓库既可以装脑电多尺度熵分析代码也可以装质谱质量误差分析代码甚至可以同时装脑电功率谱和质谱分析两套不相干的东西。我在整理代码时发现真正导致混淆的往往是当初下载或创建仓库时顺手打了“mse-analysis”这个名字没有任何上下文注释。如果你是从某个学术分享链接拿到的这个包打开一看既有“power spectrum”又有“mass error”那大概率是仓库作者把多项目脚本混在一个目录下了。我的建议是先打开代码包里的README或者主脚本看它导入的数据格式。如果读入的是.edf、.set、.mat格式的多通道时序数据那这就是脑电分析包如果读入的是.mzML、.raw、.csv格式的质谱扫描数据那这就是质谱分析包。如果两类都有那你就得自己在目录层面拆分归档别指望跑通全部代码。2. 脑电功率谱分析最容易被忽略的前置工作数据与预处理明确了“mse-analysis”在脑电语境下的真实定位是功率谱加多尺度熵之后接下来就是最枯燥但最关键的一步——数据准备。我见过太多人在这一步栽跟头拿到代码包迫不及待地导入一段脑电数据就开始跑功率谱结果出来的频谱图全是毛刺和趋势项完全没法看。问题几乎都出在预处理环节。2.1 先确认你的数据长什么样脑电数据通常是一个通道数 × 采样点数的矩阵或者是采样点数 × 通道数。如果你用的公开数据集比如DEAP、SEED、BCI Competition读取后第一件事不是算功率谱而是确认采样率和通道布局。采样率直接决定你能分析的频率上限——根据奈奎斯特定理能分析的最高频率是采样率的一半。比如采样率是256Hz那你能看到的频谱范围就是0到128Hz而脑电研究的核心频段delta、theta、alpha、beta、gamma基本都在0.5到50Hz之间所以128Hz以上的信息对常规分析意义不大。2.2 预处理管线去趋势、滤波、分段缺一不可很多新手拿到的“mse-analysis”代码包预处理部分往往很简单就是直接对原始信号做FFT。但实际可用的流程应该包含以下几步去除基线漂移和线性趋势脑电采集过程中电极与皮肤之间的极化电压、受试者出汗、设备温漂都会让信号带上一个缓慢变化的趋势项。这个趋势项在频谱上表现为极低频段的巨大能量会压低其他频段的相对功率。用detrend函数做去趋势处理是最基本的操作。带通滤波常规脑电分析的频率范围取0.5Hz到50Hz或0.5Hz到100Hz。低于0.5Hz的成分主要是漂移和伪迹高于50Hz或100Hz的成分多是肌电噪声。这个滤波的作用不是让频谱“好看”而是避免无关频段能量干扰后续的频段划分。MATLAB里直接bandpass(eegData, [0.5 50], fs)就能完成零相移滤波。剔除坏导联和坏段如果某个通道的波形明显是平线、剧烈毛刺或者大范围漂移这个通道的数据就不该参与功率谱计算。通常做法是设定一个幅度阈值或者方差阈值超过阈值的通道直接置空或插值替换。2.3 分段策略直接影响频谱分辨率功率谱的频率分辨率由参与计算的信号长度决定。计算公式很简单频率分辨率 采样率 / 参与FFT的点数。如果你用1秒的数据做FFT256Hz采样率下分辨率就是1Hz用4秒的数据分辨率可以到0.25Hz。而脑电的alpha节律8-13Hz和theta节律4-8Hz分界线很近频率分辨率太低的时候频谱图上频段边界会糊在一起很难精确算各频段能量。所以一套合理的处理流程是把连续脑电按2秒或4秒一段切分每次分析的时间窗口叫epoch对每个epoch分别计算功率谱再把所有epoch的功率谱做平均。这其实就是Welch方法的思路加窗分段、分别FFT、再平均用方差换稳定。MATLAB里pwelch函数就是这么做的[pxx, f] pwelch(eegData, window, noverlap, nfft, fs)一行调用就能得到平滑的功率谱密度估计。3. 从原始EEG到可解释的频段功率谱完整代码管线预处理做完以后真正的核心计算才开始。这一段我会给出一个相对完整、可以直接照抄改写的MATLAB实现。它不是从某个现成包里拿来的而是我根据自己的项目经验整理的一套清晰流程。3.1 核心计算代码Welch功率谱估计与频段功率汇总function powerTable compute_band_power(eegData, fs, channels) % 输入 % eegData - 通道数x采样点数的矩阵 % fs - 采样率单位Hz % channels- 通道名称元胞数组 % 输出 % powerTable - 各通道在各频段的绝对功率与相对功率表 % 定义频段边界单位Hz bands struct(... delta, [0.5 4], ... theta, [4 8], ... alpha, [8 13], ... beta, [13 30], ... gamma, [30 50]); bandNames fieldnames(bands); nBands numel(bandNames); % 去掉每个通道的线性趋势 eegDetrended detrend(eegData); % detrend默认按列操作所以先转置 % 零相移带通滤波 0.5-50Hz eegFiltered bandpass(eegDetrended, [0.5 50], fs); % 分段参数窗长2秒50%重叠 windowLen 2 * fs; noverlap round(0.5 * windowLen); nfft 2^nextpow2(windowLen); % 保证FFT点数足够频率分辨率约0.5Hz nChannels size(eegFiltered, 1); powerTable table(); for ch 1:nChannels % 对单通道数据做Welch功率谱估计 [pxx, f] pwelch(eegFiltered(ch, :), hann(windowLen), noverlap, nfft, fs); % 计算各频段绝对功率对功率谱密度积分 absPower zeros(1, nBands); for b 1:nBands idx f bands.(bandNames{b})(1) f bands.(bandNames{b})(2); % 用梯形积分近似频段能量 absPower(b) trapz(f(idx), pxx(idx)); end % 相对功率各频段占总功率的百分比 totalPower sum(absPower); relPower absPower / totalPower * 100; % 存入表格 row table(); row.Channel {channels{ch}}; for b 1:nBands row.([Abs_ upper(bandNames{b})]) absPower(b); row.([Rel_ upper(bandNames{b})]) relPower(b); end powerTable [powerTable; row]; end end这段代码表面上不长但有几个值得说的设计考虑为什么要用Welch而不是直接FFT直接对整段数据做FFT得到的是整段信号的频率成分平均一旦数据里混入眨眼伪迹或运动伪迹频谱会被局部突发能量污染。Welch方法把数据切成长度较短的段计算再做平均相当于把偶发伪迹的影响“稀释”掉。代价是频率分辨率会降低但脑电频段划分本身就有一定的宽容度0.5Hz的分辨率完全够用。为什么用梯形积分而不是直接求和pwelch返回的是功率谱密度PSD单位是μV²/Hz某个频段的能量应该是这个频段范围内PSD的积分。直接求和是等间隔采样下积分的近似但用trapz更精确尤其是频段边界落在FFT频率点之间时梯形积分能更好地逼近真实面积。为什么窗函数选hannFFT默认相当于加了矩形窗矩形窗的频谱旁瓣衰减只有约13dB会让强频率成分的能量“泄漏”到邻近频段。汉宁窗hann旁瓣衰减提高到约31dB能显著减少谱泄漏。hamming窗也有类似效果但hann窗在主瓣宽度和旁瓣衰减之间更均衡是脑电功率谱计算中最常用的选择。3.2 从功率值到可解读的结论算完各频段绝对功率和相对功率后直接面对的问题就是怎么解读。绝对功率受个体差异和电极阻抗影响比较大同样的alpha节律不同受试者头皮厚度不同功率绝对值可能差好几倍。所以组间比较时通常看相对功率即各频段占总功率的百分比这相当于做了个简单的归一化。表格式的输出会非常直观通道绝对δ绝对θ绝对α绝对β相对δ(%)相对θ(%)相对α(%)相对β(%)Fz12.38.725.15.224.117.049.110.2Cz15.69.218.36.831.218.436.613.6按我的经验做认知实验时重点关注alpha和theta的相对功率变化做睡眠分析时delta和theta是核心频段做癫痫检测时则要细化到次频段比如alpha1是8-10Hzalpha2是10-13Hz。千万别拿一套固定的频段划分应付所有场景。3.3 一个实测中的意外基线段的选取第一次跑这套流程时我被一个看似简单的问题卡了很久——基线段的选取。做静息态脑电分析时基线通常是闭上眼睛、放松状态下的记录段用来和任务态做对比。但不同被试的闭眼静息态alpha功率差异巨大哪怕用了相对功率组间方差还是很大。后来我改用“任务前1分钟的静息段作为基线”并且对每个通道单独做基线归一化效果才稳定下来。这个细节在现成的“mse-analysis”代码包中往往没有得自己在实验设计阶段就定好。4. 多尺度熵MSE功率谱之外的复杂度视角如果你手头的代码包名字是“mse-analysis”那它大概率不只包含功率谱计算还会包含多尺度熵算法。既然标题里带了这个词这一节就把多尺度熵的原理和MATLAB实现完整讲透。4.1 多尺度熵到底在算什么多尺度熵的基本流程分两步粗粒化和样本熵计算。**粗粒化coarse graining**是把长度为N的信号以尺度因子τ为窗长依次取平均。尺度τ1时就是原始信号τ2时把相邻两点平均成一个点信号长度减半τ3时三点平均成一点以此类推。这一步的物理含义是把原始信号在不同时间分辨率上重新“采样”尺度越大看到的时间窗口越宽高频细节被平滑掉留下的是低频趋势性变化。**样本熵Sample Entropy**衡量的是信号中“新模式出现”的概率。具体来说设定嵌入维度m和相似容差r统计信号中长度为m的子序列和长度为m1的子序列中彼此距离小于r的比例再取负对数。样本熵越大说明信号中产生新模式的概率越高信号越复杂。把每个尺度上的样本熵连成一条曲线就是多尺度熵曲线。健康成年人在脑电上通常表现为尺度1到尺度5区间熵值较高然后随尺度增大缓慢下降而病理状态或衰老状态下这个“复杂性资源”会减少熵值整体下沉或快速衰减。4.2 MATLAB实现多尺度熵的注意点多尺度熵的MATLAB实现网上能找到很多版本但真正用起来有几个细节决定结果可靠性参数选择直接影响结论嵌入维度m一般取2容差r取0.1到0.25倍原始信号标准差。r取太大样本熵会趋近于0区分度丢失r取太小样本熵对噪声过于敏感统计波动很大。我通常取r 0.15 * std(signal)在公开数据集上效果比较稳定。信号长度要够样本熵对数据长度很敏感数据太短统计值不可靠。经验法则原始数据的长度至少是最大尺度τ的10的m1次方量级以上。换句话说如果你要算到尺度10、m2那原始数据至少需要约10*10²1000个点。在脑电分析里如果采样率是250Hz4秒的数据才1000个点勉强够用最好能用到20秒以上的连续数据。function mseCurve compute_mse(signal, maxScale, m, r) % 多尺度熵计算 % 输入 % signal - 1xN的输入信号 % maxScale - 最大尺度因子 % m - 嵌入维度通常取2 % r - 相似容差通常取0.15*signal标准差 % 输出 % mseCurve - 1xmaxScale的多尺度熵值 N length(signal); mseCurve zeros(1, maxScale); for tau 1:maxScale % 粗粒化 coarseLen floor(N / tau); coarseSignal zeros(1, coarseLen); for i 1:coarseLen segment signal((i-1)*tau 1 : i*tau); coarseSignal(i) mean(segment); end % 计算样本熵 mseCurve(tau) sample_entropy(coarseSignal, m, r); end end function se sample_entropy(signal, m, r) N length(signal); if N 10^(m1) error(信号长度不足样本熵结果不可靠); end % 构造m维和m1维模板向量 templates_m zeros(N-m, m); templates_m1 zeros(N-m, m1); for i 1:N-m templates_m(i, :) signal(i : im-1); templates_m1(i, :) signal(i : im); end % 计算匹配对数用距离小于r作为匹配条件 % 这里采用优化的归一化方法避免双重循环过慢 A 0; % 匹配长度为m1的对数 B 0; % 匹配长度为m的对数 % 用pdist2算距离矩阵注意排除自身匹配 distM pdist2(templates_m, templates_m); distM(1:N-m1:end) inf; B sum(distM(:) r) / 2; % 距离矩阵对称除以2去除重复计数 distM1 pdist2(templates_m1, templates_m1); distM1(1:N-m:end) inf; A sum(distM1(:) r) / 2; if A 0 || B 0 se NaN; % 无匹配时返回NaN避免log(0) else se -log(A / B); end end这个实现里用的pdist2是我后来优化过的方式。网上很多版本的样本熵实现用三层for循环数据长度到几千点的时候跑一个尺度序列要等半天。用距离矩阵代替代数循环计算复杂度从O(N²×m)降到O(N²m×N²)速度提升非常明显。我项目里有一段5万点的脑电数据旧代码跑完20个尺度用了将近两分钟优化后5秒多就出结果了。4.3 功率谱和多尺度熵怎么配合使用功率谱和多尺度熵并不是非此即彼的关系它们各自捕捉信号的不同特征。在实际项目中如果只做静息态差异分析功率谱的相对功率指标通常就够用了但如果你关注的是认知负荷变化、意识水平变化这类非线性动态特征多尺度熵就有不可替代的优势。举个例子我在处理一组持续注意力任务的数据时发现随着任务时间延长被试的alpha相对功率明显上升——信号变得“更同步、更放松”与此同时多尺度熵在尺度2到尺度4上显著下降——信号变得“更规则、复杂度降低”。两者反映的其实是同一个神经状态变化的不同侧面功率谱看到的是频域能量重新分配多尺度熵看到的是时间序列可预测性增加。把这两种指标联合起来建模比单独用任何一个都能更好地预测行为表现。5. 如果真的遇到质谱分析代码怎么识别和离散现在回到标题里的“质谱分析”。前面说了在脑电语境下“mse”理解成多尺度熵更合理但如果你拿到的代码包里真的有质谱相关的脚本又该怎么处理先确认质谱分析的代码特征。质谱数据的基本格式是一个二维数组扫瞄时间或保留时间一维质荷比m/z一维强度值填充在矩阵里。识别质谱代码的典型特征是脚本里会出现mz、intensity、peak picking、mass calibration这类变量名或函数名。MATLAB里做质谱分析一般会用到mspeaks、msalign这些生物信息学工具箱的函数。如果你的“mse-analysis”包里出现这些函数调用那它确实包含质谱分析模块。遇到这种情况我的经验是立即把仓库按功能拆分。脑电功率谱和质谱分析虽然都叫“信号处理”但数据来源、预处理方式、评价指标完全不同硬放在一个目录里只会让后续维护越来越痛苦。拆分的依据不是代码风格而是I/O边界——看每个脚本读入什么格式的文件输出什么类型的文件凡是读入输出互不依赖的脚本就放进独立的子目录。至于质谱分析本身如果确实需要做MATLAB里最常跑通的任务是峰检测和误差分布统计。mspeaks可以检测质谱峰msalign做质量轴校准然后用拟合误差分布的方式评估质量准确度。这部分代码本身并不复杂但它和脑电分析没有半点交集硬要在一个项目里混合只会把README写得越来越长、越来越让人看不懂。写在最后的一个小经验这段时间整理代码最大的收获是不要相信文件夹名字也不要相信下载时随手起的项目名你的代码包里很可能同时躺着“脑电功率谱”、“多尺度熵”、“质谱分析”三种截然不同的东西而它们共享了一个叫“mse-analysis”的帽子。面对这种情况先花十分钟把数据格式和函数依赖理清楚比急着跑出第一张频谱图重要得多。至于脑电功率谱和多尺度熵这两部分按文中的代码管线走一遍再根据你自己的采样率和频段需求调整参数通常就能稳定出结果了。如果你只是为了一时方便把命名搞得更乱那我可以保证三个月后的你一定会回来感谢那个现在愿意花时间拆分的自己。本文还有配套的精品资源点击获取