Matlab脑电功率谱与多尺度熵分析实战：从Welch到MSE

简介开源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也就是每个频谱点之间的间隔。比如采样率Fs1000 Hz窗口长度设为1000点频率分辨率约为1 Hz窗口长度500点分辨率就是2 Hz。分辨率越高越能区分相近的频率成分但窗口越短参与平均的段数越多谱线越平滑。脑电频率成分里alpha峰通常在8到13Hz之间相邻峰之间差距通常是几赫兹所以频率分辨率做到1 Hz左右基本够用。重叠比例的典型值在50%到75%之间。50%重叠是教科书常用值但我在处理EEG时倾向于75%因为脑电数据段往往不会无限长75%重叠比50%能产生更多分段方差进一步压低代价是计算量略增。NFFT是指每段做FFT时补零后的点数它不影响实际频率分辨率只影响谱图看起来细不细。NFFT设成大于窗口长度画出图形会变“细腻”但没有增加真实信息。下面这个表是我常用的参数参考参数我的常用值取值范围影响窗口长度Fs1秒0.5s到4s窗口越长频率分辨率越高但分段数变少重叠50%-75%0到95%重叠越高谱越平滑但相邻分段相关性越高NFFTmax(1024, 窗口长度)大于等于窗口长度只影响插值点数不影响真实分辨率窗函数HammingHamming, Hanning, Blackman抑制频谱泄漏Hamming是脑电常用折中窗口类型容易被忽略。如果信号存在很强的工频干扰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点匹配继续看m1点是否还匹配最后用匹配概率的负对数定义熵值。匹配越多信号越规则熵越低。功率谱是频域全局指标样本熵是时域非线性指标这两者不是替代关系而是互补关系所以很多分析会同时提取这两组特征。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:im-1); for j i1:N-m1 tmpr x(j:jm-1); if max(abs(tmpl - tmpr)) r B B 1; end end end % 第二轮循环计算m1维匹配数 for i 1:N-m tmpl x(i:im); for j i1:N-m1 tmpr x(j:jm); 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统计m1维模板对匹配数量。理性理解就是信号在更高维度仍然保持相似说明信号越规则熵越低。参数说明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)*s1 : 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用不重叠窗口tau1指的是窗口起点不滑动当tau2时是重叠窗口结果会多于标准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让高尺度序列因方差减小而显得规则。解决在自己文章里写清楚实现细节粗粒化用的是不重叠平均每个尺度重新计算r0.15×stdm2。如果非要和某篇文献的绝对值比找那篇论文有没有公开代码没有的话就自己用正弦信号噪声做合成测试跑通后再和文献的方向趋势对比不想对比绝对值。5.5 批处理时一堆通道算出NaN坏点和常数段捣乱现象循环处理多组数据时某个通道的sampen返回NaN或者pwelch结果里出现Inf。原因通道数据里有NaN值或长时间不变段常数段的标准差为0r变成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文件供后续统计和模型训练使用。导出时记得把被试编号、条件标签、通道名称一起存下来否则后面想按分组做置换检验又得重新打开数据才能对应上。我自己吃过这个亏特征导出时没带上条件标签后期只能按照文件名列表猜顺序浪费了一天时间。养成任何中间结果落盘都附带元数据的习惯后面会省很多事。希望这个流程能帮你在自己的脑电数据上少走这几趟弯路。本文还有配套的精品资源点击获取