简介：开源MATLAB多尺度熵（MSE）分析工具包，面向脑电信号处理研究者，聚焦EEG数据的熵特征提取与功率谱分析，适合生物医学工程、认知神经科学领域具备一定MATLAB基础的研究人员。包内共40个文件，其中26个m脚本/函数承担核心计算，7个txt文档记录数据清单与说明，另有c辅助程序及电极位置文件，整体仅23KB，结构紧凑。工具实现了从Brain-atlas III二进制数据加载、20通道微伏级EEG读取、粗粒度化、多尺度熵计算到批量MSE与频谱分析、结果绘图核验的完整流程，并附有示例列表与缺失数据说明。目前已有342人学习下载，可直接借鉴其批处理架构与函数组织方式，用于多通道脑电信号的复杂度分析与频域考察。

1. 这个标题最迷惑人的地方：质谱分析为什么和脑电功率谱挂在了一起

刚看到这个标题，我愣了一下：Matlab脑电功率谱代码、mse-analysis、质谱分析，这三个词怎么会出现在同一个项目名里？先说结论：脑电的谱分析和质谱分析完全是两回事。质谱是分子电离后按质荷比做分离检测，而脑电功率谱是看EEG在不同频率上（delta、theta、alpha、beta、gamma）的能量分布。标题里真正要落地的是功率谱分析，mse-analysis指的是多尺度熵分析（Multiscale Entropy），它是脑电非线性特征里的常用指标，跟质谱没有任何关系。这篇文章我就把这条链路上的东西讲清楚：功率谱怎么用Matlab算出来，mse-analysis怎么和功率谱搭配，以及那些一跑就崩、一画就糊的坑在哪。适合正在整EEG数据、需要把脑电转成可量化特征做分组对比的开发者。

2. 先立住功率谱理论：Welch估计里的窗口、重叠、NFFT是怎么共同决定结果的

2.1 功率谱密度和FFT的关系：不能直接拿FFT幅值谱当功率谱

很多新手拿到脑电信号后的第一反应是做FFT，然后直接把幅度谱画出来。这个做法在脑电数据分析里不够严谨。FFT得到的是离散谱线，它反映的是每个频率分量的幅度和相位，但脑电本质上是一个随机过程，幅度谱对噪声和截断非常敏感，拿一段几秒钟的EEG做FFT，画出来的谱线会非常毛糙，而且不同时间段结果差异很大。功率谱密度（PSD）的物理含义是信号在该频率附近单位频率带宽上的功率，它需要对多段信号做平均，才能得到稳定的估计。

脑电功率谱分析的目标不是追求“精确画出一条谱线”，而是获得一个可比较、可重复的频率能量分布。所以正确的做法不是做一次FFT，而是做基于分段平均的谱估计。Matlab里最常用的就是pwelch函数，它实现的正是Welch方法：把长信号切成多段，每段加窗后做FFT，再对所有段的功率谱取平均，从而降低方差。pwelch的调用参数里，窗口长度、窗口重叠比例、FFT点数这三个参数直接影响谱的形状和分辨率。

2.2 Welch法为什么是脑电分析的主流选择

EEG作典型特征是非平稳，躺着静息态和睁眼任务态的波形能差出一大截。Welch法用分段平均的方式把“局部波动”抹平，得到的是整个记录片段内的平均功率分布，这对静息态EEG分析、睡眠分期、认知任务对比都更稳定。另一个选择是AR模型（自回归模型）谱估计，它的好处是在短数据上也能得到平滑谱，但AR阶数选择比较玄学，不同通道和不同被试可能需要不同阶数，批处理时容易翻车。还有多窗口谱估计（Multitaper），它对信号的高斯性有假设，EEG离这个假设还有距离，而且Multitaper在Matlab基础工具箱里没有现成函数，需要额外工具包，落地成本高。

我平时做脑电分析，默认就是Welch。它不要求信号严格平稳，计算效率也够，pwelch一行调用就能得到PSD，后面接频带积分、分组统计都很顺手。如果只是为了某个几十秒的静息态片段做粗略看图，窄带FFT也可以接受，但只要涉及多被试、多通道的批量处理，Welch是更安全的选择。

2.3 三个参数怎么定：窗口长度、重叠比例、NFFT三者的关系

先看频率分辨率。Welch法的频率分辨率约等于Fs / window_length，也就是每个频谱点之间的间隔。比如采样率Fs=1000 Hz，窗口长度设为1000点，频率分辨率约为1 Hz；窗口长度500点，分辨率就是2 Hz。分辨率越高，越能区分相近的频率成分；但窗口越短，参与平均的段数越多，谱线越平滑。脑电频率成分里alpha峰通常在8到13Hz之间，相邻峰之间差距通常是几赫兹，所以频率分辨率做到1 Hz左右基本够用。

重叠比例的典型值在50%到75%之间。50%重叠是教科书常用值，但我在处理EEG时倾向于75%，因为脑电数据段往往不会无限长，75%重叠比50%能产生更多分段，方差进一步压低，代价是计算量略增。NFFT是指每段做FFT时补零后的点数，它不影响实际频率分辨率，只影响谱图看起来细不细。NFFT设成大于窗口长度，画出图形会变“细腻”，但没有增加真实信息。下面这个表是我常用的参数参考：

窗口类型容易被忽略。如果信号存在很强的工频干扰，pwelch默认的Hamming窗旁瓣衰减只有约-45dB，可能压不住；这种情况我建议先用陷波滤波把干扰去掉，而不是指望换窗函数。Blackman窗旁瓣衰减更彻底，但主瓣更宽，会把靠近的alpha峰和theta峰左右抹开，实际效果不一定好。

2.4 一个鲜为人知的习惯：先肉眼检查PSD形状再信频段数字

写代码之前如果先想清楚这几件事，后面少走不少弯路。算出来的PSD曲线主要形状应该在低频有一个缓慢上升的背景，theta和alpha区域出现凸起，高频呈近似幂律下降。如果看到全频段平得一条直线，很可能信号本身是白噪声；如果低频0到1Hz处有一根冲天的柱子，说明直流漂移和缓慢基线漂移没有去除干净。这些经验判断虽然不“定量”，但能帮你第一时间发现预处理是否出了问题。

3. Matlab脑电功率谱代码实现：从原始通道到delta/theta/alpha/beta频带功率

3.1 读入EEG数据与去漂移预处理

假设现在手上有一份.mat格式的EEG数据，里面是名为data的通道×时间矩阵，采样率Fs存为变量。第一步不是直接上pwelch，而是先做基本的预处理：去除均值、去掉低频频漂、剔除坏段。下面这段代码把数据从文件加载出来，做一次基本清洗：

% 加载数据 load('eeg_test.mat', 'data', 'Fs'); % data: 通道数×采样点数, Fs: 采样率 data = data - mean(data, 2); % 每个通道减均值，去掉直流分量 % 高通滤波：截止频率0.5Hz，去除基线漂移 [b, a] = butter(2, 0.5/(Fs/2), 'high'); data_f = filtfilt(b, a, data')'; % filtfilt零相位，避免相位偏移

这段代码的逻辑是先去除直流分量，再用0.5Hz的高通滤波器消除大体上的基线漂移。注意我把数据转置后做filtfilt，是因为filtfilt默认沿第一个非单一维度滤波，原始数据通道在行方向，转置后时间轴变到行方向，滤波结果再转置回来。

参数说明：高通截止频率不要设得太高，比如直接5Hz高通会砍掉delta频段，后面算delta带功率全是失真值。0.5Hz是一个常用折中，能滤掉慢漂移但保留delta波。如果是睁眼任务态，也可以设1Hz，但静息态和睡眠分析保留delta很重要，所以0.5Hz更保险。

3.2 主力pwelch代码：一键输出全频段功率谱密度

清洗完成后开始计算PSD。用pwelch时我习惯把窗口长度、重叠、NFFT定义成变量，方便后面批量调参：

% 参数区 win_len = floor(Fs * 1); % 1秒窗口 overlap = floor(win_len * 0.75); % 75%重叠 nfft = 2^nextpow2(win_len); % 补零到1024或更大 [pxx, freq] = pwelch(data_f', win_len, overlap, nfft, Fs); % px: 频率点数×通道数，每列是一个通道的PSD % freq: 对应的频率轴

这段代码里data_f是通道×时间，pwelch需要把每一列当作一个观测信号，所以我把数据转置成时间×通道输入。返回值pxx是频率点数×通道数，freq是频率向量。pwelch内部的逻辑是把每列信号切段、加Hamming窗、做FFT、再平均，输出PSD单位是uV^2/Hz，不是绝对功率值。

参数说明：窗口长度设成1秒，频率分辨率是1Hz；如果某个通道的数据只有几秒钟，窗口可以缩短到0.5秒，但频率分辨率会降到2Hz，细小的alpha峰可能看不明显。重叠75%是我处理30秒静息态EEG的常用配置，能得到大约117段参与平均；如果数据长度到几分钟，重叠降到50%也能得出稳定结果，还能节约计算时间。

3.3 提取频带功率并整理成特征表

拿到PSD之后，真正用于后续统计的是各频段的积分功率，而不是整条曲线。下面这段函数把delta、theta、alpha、beta、gamma的功率分别算出来，并除以频率带宽得到平均功率密度，保证不同频段之间可比较：

function bands = eeg_band_power(pxx, freq) % 输入: pxx 频率×通道; freq 频率轴 band_def = struct(... 'delta', [0.5 4], ... 'theta', [4 8], ... 'alpha', [8 13], ... 'beta', [13 30], ... 'gamma', [30 50]); fields = fieldnames(band_def); for i = 1:numel(fields) fname = fields{i}; f_lo = band_def.(fname)(1); f_hi = band_def.(fname)(2); idx = freq >= f_lo & freq < f_hi; % 积分功率: 对PSD求和再乘频率分辨率 band_power = sum(pxx(idx, :), 1) .* (freq(2)-freq(1)); band_density = band_power / (f_hi - f_lo); % 平均功率密度 bands.(fname) = band_density; end end

这里我用梯形积分近似，如果PSD曲线变化剧烈，sum平均法会低估真实积分，但脑电PSD频带内相对平滑，这种简化在组间对比时影响不大。重点是我把gamma上限定在50Hz，不做50Hz以上部分，因为工频干扰和肌电伪迹集中在50Hz附近，没有做严格陷波时，gamma上限再高会把噪声也带进特征里。

调用方式是：band_result = eeg_band_power(pxx, freq)，返回的bands是一个结构体，每个字段是通道×1的向量。直接再拼成矩阵就能当分类特征或统计输入。频带划分没有绝对标准：delta上界有从4Hz到4.5Hz的划分，beta有从12Hz起的版本，写死参数前先看文献里自己领域的主流定义。

3.4 画图：频谱曲线和能量topoplot的取舍

如果不做topoplot，只画一条通道的PSD曲线，下面这段就够用：

figure; plot(freq, 10*log10(pxx(:, 1)), 'LineWidth', 1.2); xlabel('Frequency (Hz)'); ylabel('Power (dB)'); xlim([0 60]); grid on;

把PSD转成dB是为了高动态范围看得清楚，脑电功率谱在低频和高频相差可达几个数量级，线性坐标下高频段几乎压成一条直线。坐标上限限制在60Hz，避免把不必要的频率噪声画进去。如果后续要发论文，建议用parpool或者循环把每个通道的PSD画到同一张图上，观察是否有某个通道明显偏离群体，那是坏通道的常见信号。

4. mse-analysis代码：多尺度熵和功率谱怎么搭配着用

4.1 为什么功率谱还不够：mse-analysis补的是什么

功率谱给出的是“能量多少”，但脑电还有另一个维度——信号的不规则性和复杂度。比如两种状态下alpha带功率可能没差异，但波形规律程度明显不同。mse-analysis（多尺度熵分析）就是用来量化这种复杂度的，它通过多尺度粗粒化，计算每个尺度下的样本熵，刻画脑电在不同时间尺度上的规则性。认知负荷、睡眠阶段、癫痫发作前后的EEG复杂度变化，都是功率谱不容易直接区分的。

样本熵的计算思想是看时间序列里有多少“模板向量”在容差r内互相匹配，m点匹配继续看m+1点是否还匹配，最后用匹配概率的负对数定义熵值。匹配越多，信号越规则，熵越低。功率谱是频域全局指标，样本熵是时域非线性指标，这两者不是替代关系，而是互补关系，所以很多分析会同时提取这两组特征。

4.2 样本熵计算的核心代码

Matlab没有内置样本熵函数，需要自己实现。下面是基础版本，输入是单通道信号，输出是样本熵值：

function se = sampen(x, m, r) % x: 行向量信号 % m: 嵌入维数，EEG常用1或2 % r: 容差，通常取信号标准差的0.15~0.25倍 x = x(:)'; N = length(x); B = 0; A = 0; % 第一轮循环，计算m维匹配数 for i = 1:N-m tmpl = x(i:i+m-1); for j = i+1:N-m+1 tmpr = x(j:j+m-1); if max(abs(tmpl - tmpr)) <= r B = B + 1; end end end % 第二轮循环，计算m+1维匹配数 for i = 1:N-m tmpl = x(i:i+m); for j = i+1:N-m+1 tmpr = x(j:j+m); if max(abs(tmpl - tmpr)) <= r A = A + 1; end end end if B == 0 || A == 0 se = NaN; else se = -log(A / B); end end

这个实现用了双重循环，数据长会非常慢，30秒EEG按1000Hz采样就是30000个点，跑一个通道要数分钟。实际使用时我一般先对信号做降采样到250Hz，再把数据切到10秒以内，或者写向量化版本。代码的原理：B统计所有m维模板对在容差r内的数量，A统计m+1维模板对匹配数量。理性理解就是信号在更高维度仍然保持相似，说明信号越规则，熵越低。

参数说明：m取2是EEG的常用值，m取1对长信号更稳定但会丢失部分非线性信息。r取0.15乘以信号标准差是比较标准的，取得太小匹配对数少，熵值波动大；取得太大会让所有信号都显得规则，区分度下降。注意r是随信号标准差变化的，所以原始信号只要不是被放大缩小，r会自动适应，这个性质对多被试比较很友好。

4.3 多尺度粗粒化与MSE曲线

多尺度熵的“多尺度”指的是把原始序列每隔τ个点不重叠求平均，形成新的粗粒化序列，再在每个序列上算样本熵，得到一条MSE曲线。粗粒化的代码如下：

function mse = compute_mse(x, max_scale, m, r, tau) % x: 原始信号 % max_scale: 最大尺度，一般10~20 % tau: 按文献可设1或2，不设则用默认1 mse = zeros(1, max_scale); for s = 1:max_scale % 粗粒化 len_coarse = floor(length(x) / s); y = arrayfun(@(k) mean(x((k-1)*s+1 : k*s)), 1:len_coarse); y = y(:)'; % 每尺度下重新计算r，因为粗粒化后标准差变了 r_used = r * std(y); mse(s) = sampen(y, m, r_used); end end

这里的坑在于：r如果不随尺度更新，会在粗粒化后的低幅序列上用原来的绝对值r，导致匹配数剧增、熵值整体变小，MSE曲线变成单调递减，失去区分意义。所以我在每个尺度里重新算r_used。tau参数在很多开源实现里是粗粒化时相邻窗口之间的滑动步长，普通MSE用不重叠窗口，tau=1指的是窗口起点不滑动，当tau=2时是重叠窗口，结果会多于标准MSE曲线，论文里必须写明，否则别人复现不出你的曲线。

计算MSE的耗时比样本熵还高，建议通道多时先用循环逐个跑，并打印进度，一个20尺度、数据长度1万点的单通道大约需要几十秒，长度翻倍耗时指数上升。如果有并行计算工具箱，用parfor替换for是最省事的加速方式。

4.4 把MSE和功率谱拼成特征向量

把功率谱特征和MSE特征拼接成一行特征向量，是分类和统计分析的常用做法。下面是整合代码：

% band_result: 5个频段的通道×1特征 % mse_result: 通道×max_scale，每个通道一段MSE曲线 [m, s] = size(mse_result); features = zeros(m, 5 + s); features(:, 1:5) = [band_result.delta, band_result.theta, ... band_result.alpha, band_result.beta, band_result.gamma]; features(:, 6:end) = mse_result; % 把特征矩阵统一缩放到零均值单位方差，防止频带功率和熵值量纲差异过大 features_z = zscore(features, [], 1);

拼接前先留意数值尺度：delta带功率可能是几uV^2，熵值是0到2的小数，直接丢给分类器会默认权重偏向功率。zscore按每列标准化是基本操作。如果后续做统计检验而不是分类，则不一定需要标准化，但分组对比时也要先画箱线图看分布形态，避免用错检验方法。

5. 脑电谱分析避坑指南：5个高频问题与排查办法

5.1 低频能量异常高：0到1Hz处突起一个尖峰

现象：pwelch画出来的PSD曲线在0Hz到1Hz之间有一个明显的能量尖峰，delta带功率异常大，其它频段的特征都被压低。

原因：脑电信号里存在直流偏置和电极缓慢漂移，这个分量经过FFT后会散布到超低频区域，形成假的高能量。另一个常见原因是高通滤波器用的截止频率低于0.5Hz，没能压住漂移。

解决：进入pwelch之前先做0.5Hz高通滤波，filtfilt零相位滤波优先于任何谱估计操作。如果数据里有大幅运动伪迹，高速滤波截止频率还可以提到1Hz，但前提是后续不做delta带分析。还有一个习惯：pwelch之前先plot原始波形，看到明显的“掉毛”式漂移就说明预处理没到位。

5.2 频谱毛刺太多：谱线像锯齿，频带功率不稳定

现象：同一段数据重复计算，功率谱曲线很粗，频带功率数值在两次运行之间差别很大。

原因：分段太少，Welch平均没有发挥降方差的功效。常见是窗口设置过长，30秒数据只分到2-3段；另一可能是没有设置重叠比例，默认不重叠导致段数直接减半。

解决：把窗口长度缩到1秒，重叠提到75%，这样30秒数据大约能产生117个分段，谱线平滑度明显改善。如果数据本身只有5秒，那就不要追求高频率分辨率，窗口缩短到0.5秒，分辨率降到2Hz，但谱会更稳定。

5.3 alpha峰不在8到13Hz：峰顶偏移或直接消失

现象：睁眼闭眼实验设计里，闭眼状态应该在约10Hz处看到alpha峰，但算出来峰位置飘到7Hz或14Hz，或者根本没有峰。

原因：工频干扰和肌电污染是最常见的干扰源；另外谱泄漏会把alpha峰的能量向两侧展宽，看起来峰顶变平。有些数据集里alpha峰本来就不在10Hz，年龄和认知状态都会改变峰频率。

解决：先画PSD看有没有50/60Hz工频尖刺，有就加陷波滤波。加窗时用Hamming窗而不是默认矩形窗。如果你用的是固定频带汇总法，不看峰位置，那就保证8到13Hz的积分包含实际峰能量。想精确定位peak位置，可以用findpeaks在8到13Hz窗口内找局部最大点，再加峰值前后的左右半功率宽度特征。

5.4 MSE尺度0的结果和文献对不上：数值差一倍或方向相反

现象：用别人的论文参数计算MSE，第一个尺度的熵值和文献相差很大，或者前人曲线随尺度先升后降，自己算出来一直降。

原因：粗粒化各自实现方式不一样。有些论文用滑动平均加降采样，有些用不重叠平均，还有的用“尺度因子”和“时间延迟”区分，这三种方法在同一数据上结果完全不同。另一个原因是没有随尺度更新r，让高尺度序列因方差减小而显得规则。

解决：在自己文章里写清楚实现细节：粗粒化用的是不重叠平均，每个尺度重新计算r=0.15×std，m=2。如果非要和某篇文献的绝对值比，找那篇论文有没有公开代码；没有的话，就自己用正弦信号+噪声做合成测试，跑通后再和文献的方向趋势对比，不想对比绝对值。

5.5 批处理时一堆通道算出NaN：坏点和常数段捣乱

现象：循环处理多组数据时，某个通道的sampen返回NaN，或者pwelch结果里出现Inf。

原因：通道数据里有NaN值或长时间不变段，常数段的标准差为0，r变成0，匹配条件永远不成立，样本熵无定义。

解决：做样本熵前先检测：x = x(~isnan(x)); 如果去NaN后长度不到原长度一半，直接丢弃该通道。对常数段，用std(x) == 0跳过。批量循环时把坏通道索引存下来，分析完统一看，不要直接删数据，否则后面电极位置表和特征矩阵对不上。

6. 让结果更可信：合成信号验证和批处理导出

6.1 用已知频率的合成信号验证功率谱

写好的功率谱代码第一次跑出一堆数字，你并不知道对不对。我的习惯是先造一个已知频率成分的信号做回归测试：

t = 0:1/1000:30-1/1000; y = 10*sin(2*pi*10*t) + 5*sin(2*pi*20*t) + 2*randn(size(t)); [pxx, freq] = pwelch(y, 1000, 750, 1024, 1000); [~, idx] = max(pxx(freq > 8 & freq < 13)); peak_freq = freq(freq > 8 & freq < 13); fprintf('Peak within alpha: %.2f Hz\n', peak_freq(idx));

这段代码造了一个含10Hz和20Hz正弦的混合信号，加了一点白噪声。如果功率谱在10Hz附近没有出现明显峰，说明滤波、窗口或频率轴有问题。MSE同样可以做验证：白噪声的MSE曲线应该在各尺度都较高，正弦信号的MSE应该在低尺度较低，若结果反了，多半是粗粒化实现有bug。

6.2 批处理多段EEG并导出特征矩阵

到了实际项目阶段，我一般把全流程封成函数，输入是数据矩阵和参数结构体，输出是特征表：

function feat = eeg_feature_pipeline(data, Fs, params) data = data - mean(data, 2); [b, a] = butter(2, params.highpass/(Fs/2), 'high'); data_f = filtfilt(b, a, data')'; [pxx, freq] = pwelch(data_f', params.win, params.overlap, params.nfft, Fs); bands = eeg_band_power(pxx, freq); mse_result = zeros(size(data, 1), params.max_scale); for ch = 1:size(data, 1) mse_result(ch, :) = compute_mse(data_f(ch, :), params.max_scale, 2, 0.15, 1); end feat = [bands.delta, bands.theta, bands.alpha, bands.beta, bands.gamma, mse_result]; end

跑完所有被试后，把每个被试的特征矩阵按行拼成一个大矩阵保存成一个CSV或.mat文件，供后续统计和模型训练使用。导出时记得把被试编号、条件标签、通道名称一起存下来，否则后面想按分组做置换检验又得重新打开数据才能对应上。我自己吃过这个亏：特征导出时没带上条件标签，后期只能按照文件名列表猜顺序，浪费了一天时间。养成任何中间结果落盘都附带元数据的习惯，后面会省很多事。希望这个流程能帮你在自己的脑电数据上少走这几趟弯路。

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