简介:这套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秒的数据做FFT,256Hz采样率下分辨率就是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频率点之间时,梯形积分能更好地逼近真实面积。
为什么窗函数选hann?FFT默认相当于加了矩形窗,矩形窗的频谱旁瓣衰减只有约13dB,会让强频率成分的能量“泄漏”到邻近频段。汉宁窗(hann)旁瓣衰减提高到约31dB,能显著减少谱泄漏。hamming窗也有类似效果,但hann窗在主瓣宽度和旁瓣衰减之间更均衡,是脑电功率谱计算中最常用的选择。
3.2 从功率值到可解读的结论
算完各频段绝对功率和相对功率后,直接面对的问题就是怎么解读。绝对功率受个体差异和电极阻抗影响比较大,同样的alpha节律,不同受试者头皮厚度不同,功率绝对值可能差好几倍。所以组间比较时,通常看相对功率,即各频段占总功率的百分比,这相当于做了个简单的归一化。
表格式的输出会非常直观:
| 通道 | 绝对δ | 绝对θ | 绝对α | 绝对β | 相对δ(%) | 相对θ(%) | 相对α(%) | 相对β(%) |
|---|---|---|---|---|---|---|---|---|
| Fz | 12.3 | 8.7 | 25.1 | 5.2 | 24.1 | 17.0 | 49.1 | 10.2 |
| Cz | 15.6 | 9.2 | 18.3 | 6.8 | 31.2 | 18.4 | 36.6 | 13.6 |
按我的经验,做认知实验时,重点关注alpha和theta的相对功率变化;做睡眠分析时,delta和theta是核心频段;做癫痫检测时,则要细化到次频段(比如alpha1是8-10Hz,alpha2是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的子序列和长度为m+1的子序列中,彼此距离小于r的比例,再取负对数。样本熵越大,说明信号中产生新模式的概率越高,信号越复杂。
把每个尺度上的样本熵连成一条曲线,就是多尺度熵曲线。健康成年人在脑电上通常表现为尺度1到尺度5区间熵值较高,然后随尺度增大缓慢下降;而病理状态或衰老状态下,这个“复杂性资源”会减少,熵值整体下沉或快速衰减。
4.2 MATLAB实现多尺度熵的注意点
多尺度熵的MATLAB实现网上能找到很多版本,但真正用起来有几个细节决定结果可靠性:
- 参数选择直接影响结论:嵌入维度m一般取2,容差r取0.1到0.25倍原始信号标准差。r取太大,样本熵会趋近于0,区分度丢失;r取太小,样本熵对噪声过于敏感,统计波动很大。我通常取r = 0.15 * std(signal),在公开数据集上效果比较稳定。
- 信号长度要够:样本熵对数据长度很敏感,数据太短,统计值不可靠。经验法则:原始数据的长度至少是最大尺度τ的10的m+1次方量级以上。换句话说,如果你要算到尺度10、m=2,那原始数据至少需要约10*10²=1000个点。在脑电分析里,如果采样率是250Hz,4秒的数据才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^(m+1) error('信号长度不足,样本熵结果不可靠'); end % 构造m维和m+1维模板向量 templates_m = zeros(N-m, m); templates_m1 = zeros(N-m, m+1); for i = 1:N-m templates_m(i, :) = signal(i : i+m-1); templates_m1(i, :) = signal(i : i+m); end % 计算匹配对数(用距离小于r作为匹配条件) % 这里采用优化的归一化方法,避免双重循环过慢 A = 0; % 匹配长度为m+1的对数 B = 0; % 匹配长度为m的对数 % 用pdist2算距离矩阵,注意排除自身匹配 distM = pdist2(templates_m, templates_m); distM(1:N-m+1: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”的帽子。面对这种情况,先花十分钟把数据格式和函数依赖理清楚,比急着跑出第一张频谱图重要得多。至于脑电功率谱和多尺度熵这两部分,按文中的代码管线走一遍,再根据你自己的采样率和频段需求调整参数,通常就能稳定出结果了。如果你只是为了一时方便把命名搞得更乱,那我可以保证,三个月后的你一定会回来感谢那个现在愿意花时间拆分的自己。
本文还有配套的精品资源,点击获取