简介:面向生物医学工程、运动科学及康复工程研究者的Matlab资源包,用于从表面肌电信号(sEMG)中求解肌肉激活度。资源内含完整可运行的源代码与实测数据:5个文件涵盖2个m脚本(主程序与激活度计算函数)、2个.mat数据文件(普通运动状态下的sEMG数据与最大自主收缩MVC条件下的基准sEMG数据)及1个txt说明文档,整个包仅135KB,轻量便携。运行主程序即可自动读取数据并绘制肌肉激活度随时间变化的曲线,无需额外配置,适合快速验证算法或作为sEMG信号处理入门的参考样例。尤其对于正在学习肌电信号分析、需要理解MVC归一化流程的读者,这份资源提供了可直接对照学习的标准数据与代码实现。目前已有3505人学习,兼具教学演示与科研参考价值。
1. 肌肉激活度 Matlab 代码:把 sEMG 原始波形变成可复用的处理流水线
做表面肌电信号(sEMG)研究的人,几乎都在肌肉激活度计算上卡过:波形在屏幕上跳动,肉眼看得出发力与放松,但要回答“这块肌肉在某一时间窗内激活到百分之多少”就无从下手。这份 Matlab 代码与数据包干的就是这件事——把肌肉激活度计算拆成带通滤波、整流、线性包络和 MVC 归一化四步,配上一组可直接复现的数据,让原始 sEMG 一条龙变成 0 到 1 的激活度曲线。适合刚接触肌电处理的同学,也适合要做多通道扩展和算法改造的工程师。
2. 工程结构与数据约定:一用就懂的 .mat 数据和模块化脚本
解压后你会看到代码和数据的清晰划分:.m 文件处理逻辑,.mat 文件存放原始信号。这套代码在组织上有个很友好的约定:配置、预处理、激活计算、结果绘图各占一个脚本,入口只有一个 main_demo.m。如果你之前见过那种动辄几百行、滤波和绘图挤在一起的单脚本,应该能体会到这种拆分在换数据和调参数时有多省力。
2.1 目录结构与脚本职责
我拿到这类代码的习惯,是先扫一遍文件名,确认每个文件的边界。这套包反复出现的文件布局基本是这样:
| 文件/目录 | 职责 |
|---|---|
| parameter_config.m | 集中存放采样率、通道数、滤波频段、MVC窗口等全部参数 |
| preproc_sEMG.m | 对原始 sEMG 做带通滤波、全波整流和线性包络提取 |
| compute_mvc.m | 读入 MVC 标定试验,输出用于归一化的激活基准 |
| muscle_activation.m | 把线性包络除以 MVC 基准,再映射到 0~1 激活度 |
| plot_pipeline.m | 绘制原始信号、滤波后波形、包络和激活度的对比图 |
| main_demo.m | 主流程,调度上面所有函数并输出汇总结果 |
| data/ | 存放静息、MVC、任务三类试验的 .mat 文件 |
其中 preproc_sEMG 和 muscle_activation 是纯函数,输入输出都是向量,单独测试很方便。compute_mvc 依赖静态 MVC 试验文件,main_demo 再把三类文件按顺序组装成完整分析流程。这样拆分的原因很朴素:换一组受试者数据时,主脚本和算法函数完全不需要动,只要更新 data 目录和 parameter_config 里的对应项。
这里还有个设计取舍值得留意:为什么把静息、MVC、任务段拆成三个 .mat 文件,而不是合并成一个大矩阵?因为肌肉激活度是相对量,分母来自 MVC 标定,分子来自任务段,两类数据采集协议不同。拆开存放既避免脚本里反复切片,也降低把不同试验段搞混的风险。重新采集数据时,只需要保证文件命名和字段结构一致,整个流程就跑得起来。
2.2 数据文件内部格式:采样率、通道矩阵和单位陷阱
每个 .mat 文件内部约定存一个名为 data 的结构体,字段是 fs、channels、label。channels 是“通道数 × 采样点数”的矩阵,label 是元胞数组,和通道行一一对应。先看装进来看一眼再动手写算法:
loaded = load(fullfile('data', 'subject01_task.mat')); disp(loaded.data.fs); disp(loaded.data.label); [size(loaded.data.channels, 1), size(loaded.data.channels, 2)]逻辑说明:channels 的第一维是通道数,第二维是时间序列长度。fs 字段是所有滤波器设计的基础,每个截止频率都要先除以 fs/2 转成归一化频率,再传给 butter。label 的长度必须等于 channels 的行数,否则后面画图和统计时通道对不上。
参数说明:如果发现 channels 是 int16 类型,先别直接送进滤波器。常见做法是先把 A/D 原始码值转成毫伏,再进入后续处理:
if isa(loaded.data.channels, 'int16') loaded.data.channels = double(loaded.data.channels) / 32767 * 2.5; end这块是我踩过的一个隐藏坑:采集设备输出的是带符号整数,量程假设为 ±2.5V,不转换的话滤波器阈值、MVC 比较全都会失去物理意义。虽然最终激活度是相对值,单位换算不会改变归一化结果,但中间环节做峰谷检测时,量纲不对会让你误判噪声水平。
2.3 用 whos 和 fieldnames 预览陌生数据
拿到陌生 .mat 文件时,不建议直接双击 load 进工作区。上百 MB 的肌电数据一次载入会占用大量内存,改用 whos 先看文件头:
info = whos('-file', fullfile('data', 'subject01_mvc.mat')); disp(info); loaded = load(fullfile('data', 'subject01_mvc.mat')); disp(fieldnames(loaded));逻辑说明:whos 加 -file 选项时只读取文件元数据,不加载矩阵本体,特别适合快速确认变量名和维度。fieldnames 列出结构体的全部字段,配合 size 可以验证文件是否符合 2.2 节的数据约定。如果发现变量名不是 data,说明保存时用了别的变量名,写代码时不要 hardcode,而是用 fieldnames 取第一个字段做映射。
到这一步,你已经知道文件怎么摆放、数据长什么样了。接下来进入最关键的部分:算法核心长什么样,以及为什么是这三步而不是别的流程。
3. 核心算法拆解:带通滤波、线性包络与 MVC 归一化的代码实现
sEMG 原始信号里混着直流偏置、运动伪影、心电串扰和高频噪声,直接做归一化算出来的激活度非常不稳定。标准做法是先把原始信号变成平滑的线性包络,再用 MVC 试验得到的包络峰值做分母。整个链路拆成预处理、MVC 标定、激活度映射三个独立模块,任何一个都可以单独拆出去用在别的项目里。
3.1 预处理:零相位带通滤波与全波整流
表面肌电的有效能量大多集中在 20~450 Hz。低于 20 Hz 的成分主要是电极移动和直流漂移,高于 450 Hz 的多数是环境电磁噪声。用带通滤波把这两头切掉,再整流取绝对值,得到的就是只含幅度信息的半波信号。
function env = preproc_sEMG(raw, fs, bpFreq, lpFreq) % 输入:raw 为单通道原始sEMG序列 % fs 为采样率(Hz),bpFreq 为带通截止频率 % lpFreq 为线性包络低通截止频率 % 输出:env 为平滑后的线性包络 % 带通滤波:滤掉直流偏置和运动伪影 [b, a] = butter(4, bpFreq / (fs / 2), 'bandpass'); filtered = filtfilt(b, a, raw); % 全波整流:把负半周翻到正半周 rectified = abs(filtered); % 线性包络:对整流后的信号再做低通平滑 [bs, as] = butter(4, lpFreq / (fs / 2), 'low'); env = filtfilt(bs, as, rectified); end逻辑说明:butter 返回的分子分母系数是按归一化频率计算的,bpFreq 传 [20 450]、lpFreq 传 6 是这套包里的典型配置。整流用的 abs 是最直接的全波整流方式,缺点是局部能量波动较大。如果想要更平滑的包络,可以把整流替换成滑动 RMS:先对 filtered 做平方,再用 movmean 开 100ms 窗口滑动平均,最后开方,效果和 abs 整流接近,但在高噪声环境下对尖峰更不敏感。
参数说明:filtfilt 做的是零相位滤波,波形在时间轴上不会被拖动,代价是首尾会有边界效应。滤波器阶数 4 是常用折中,阶数再高容易出现数值不稳定,阶数太低则过渡带太宽,截止频率附近的噪声切不干净。
3.2 MVC 标定:滑动窗口提取最大激活基准
肌肉激活度是无量纲相对值,分母必须来自最大自主收缩试验。这里最容易犯的错是整段 MVC 文件取平均值,平均值会被发力前的静息段和发力后的回落到拖低,导致任务段激活度普遍超过 100%。正确做法是做短时滑动平均,再取滑动平均曲线上的峰值。
function [mvcVal, mvcIdx] = compute_mvc(envMVC, fs, winSec) % envMVC:MVC试验的线性包络 % winSec:滑动窗口长度(秒),默认0.5 winLen = round(winSec * fs); movingMean = movmean(envMVC, winLen); [mvcVal, mvcIdx] = max(movingMean); mvcVal = movingMean(mvcIdx); end逻辑说明:movmean 生成一条滑动平均值曲线,max 取出其中的峰值。0.5 秒窗口的设定是为了滤掉单点尖峰,同时保留 2 秒左右的稳定最大发力段。窗口太短会把一次电极抖动当成最大激活,窗口太长又会把强发力段的峰值摊平。
参数说明:函数把 mvcIdx 一并返回,目的是让你能在绘图时定位“MVC 峰值到底出现在哪个时间点”。如果 mvcIdx 出现在发力开始后的极短时间内,说明被试可能是爆发式发力而不是稳态最大发力,这会导致 MVC 基准偏高,任务段激活度被系统性压低。我一般会在采集时给 MVC 试验标记音,保证每段都持续 3 秒以上,只取中间稳定段。
3.3 肌肉激活度映射:线性与非线性模型的取舍
归一化后的包络 R 表示“当前发力相对最大发力的比值”,但 sEMG 幅度和肌肉力之间往往是非线性关系。低强度发力阶段,sEMG 增长慢于力增长;高强度阶段则相反。为了还原这种关系,常见的做法是在 R 上再做一次非线性映射。
function act = muscle_activation(env, mvcVal, A) % env:任务段线性包络 % mvcVal:由compute_mvc得到的MVC激活基准 % A:非线性系数,A=0时退化为线性模型 R = env ./ mvcVal; R = max(R, 0); % 防止整流问题导致负值 if nargin < 3 || A == 0 act = R; % 线性模型 else act = (exp(A * R) - 1) / (exp(A) - 1); act = max(act, 0); % 指数模型结果裁剪到非负 end end逻辑说明:指数模型在 R=0 时输出 0,R=1 时输出 1,曲线形状由 A 控制。A 越接近 0,模型越接近线性;A 为负时,低发力段的激活度被压低,高发力段的激活度被放大,这更符合肌肉募集的生理规律。代码里故意没有对 act 做 1 的上限裁剪,目的是保留异常值让你检查,而不是静默掩盖问题。
参数说明:A 的取值常见在 -3~0 之间,可以先跑一版 -1,再对比同一发力等级下线性与非线性模型的差异。对于等长收缩且发力平稳的任务,两者差距不大;对于快速动态动作,指数模型通常更贴近实际力输出。这个系数建议放进 parameter_config.m,不要写死在算法函数里。
3.4 主脚本如何串起整条链路
三个核心函数是零件,main_demo 负责把它们装成整机。主脚本的关键约束是:任务段和 MVC 段必须用同一套预处理参数,否则分母和分子的滤波特性不一致,激活度曲线会出现系统性偏差。
cfg = parameter_config(); mvcData = load(fullfile('data', 'subject01_mvc.mat')); taskData = load(fullfile('data', 'subject01_task.mat')); envMVC = preproc_sEMG(mvcData.channels(1, :), mvcData.fs, ... cfg.bandpassFreq, cfg.lowpassEnv); [mvcVal, mvcIdx] = compute_mvc(envMVC, mvcData.fs, cfg.mvcWindow); envTask = preproc_sEMG(taskData.channels(1, :), taskData.fs, ... cfg.bandpassFreq, cfg.lowpassEnv); act = muscle_activation(envTask, mvcVal, cfg.activationA); plot_pipeline(taskData.channels(1, :), envTask, act, mvcIdx, cfg);逻辑说明:这里只取了 channels 第一行作为单通道示例,实际多通道处理时外面再套一层 for 循环,每个通道独立计算 envMVC 和 mvcVal。要注意的是,每个通道的 MVC 必须来自该通道自己的 MVC 包络,不能拿 A 通道的分母去归一化 B 通道的分子。
参数说明:plot_pipeline 的最后一个参数 mvcIdx 只是在图上标注 MVC 峰值位置,方便肉眼判断标定是否合理。cfg 里的所有参数来自 parameter_config,想试不同滤波频段时,只改配置文件,不动算法函数。
4. 从解压到跑通:运行配置、参数修改和输出解释
算法看完就动手跑。这一章按实际调试顺序来写,从解压文件到解读输出曲线都走一遍。很多人跑不通并不缺在代码逻辑,而是缺在对数据的前提条件上。
4.1 环境准备与文件摆放
这套代码需要安装 Signal Processing Toolbox,其余是纯 Matlab 基础功能。拿到压缩包后先解压到英文路径,避免 load 阶段遇到编码问题。
# Windows 下用解压工具导出到指定目录,比如用 7-Zip 7z x 肌肉激活度Matlab代码及数据.rar -o"./muscle_activation" # macOS/Linux 下可以用 unar 替代 unar 肌肉激活度Matlab代码及数据.rar -o ./muscle_activation逻辑说明:解压完成后先确认 data 子目录和六个 .m 文件都在。然后在 MATLAB 里 cd 到解压目录,运行 main_demo 之前先执行 ver('signal') 确认工具箱可用。路径建议全英文的原因是 MATLAB 处理中文路径偶尔抽风,这不是代码缺陷,但能省半小时排查时间。
参数说明:如果你用的是便携版 MATLAB 或者 license 里没勾选 Signal Processing Toolbox,butter 和 filtfilt 都会报未定义函数,后文避坑章节会详细说。
4.2 三个必改参数:采样率、通道名、滤波频段
第一次跑通后换成自己的数据时,最常动的是 parameter_config.m 里的三个位置:
| 参数 | 默认值示例 | 改的时候看什么 |
|---|---|---|
| fs | 1000 | 改成采集时真实采样率,滤波器归一化频率随之变化 |
| chName | {'Biceps','Triceps'} | 改成电极安放的肌肉名,绘图和输出统计用 |
| bandpassFreq | [20 450] | 按肌肉类型和电极特征调,别盲目扩大频带 |
fs 不改是很多人翻车的点。如果实际采样率是 2000 Hz,配置里还写 1000,滤波器算出来的实际截止频率就翻了一倍,而且代码不会报错。我一般会在 parameter_config.m 开头加一个断言:如果 fs 小于 500 或大于 10000 就给出警告。
bandpassFreq 不建议低于 10 Hz,那会把运动伪影放回来;也不建议高于 500 Hz,高频噪声增益会被抬起来。在改参数前,可以用 pwelch 画一段频谱确认能量集中区:
[pxx, f] = pwelch(rawEMG, hann(512), [], [], fs); plot(f, 10*log10(pxx));逻辑说明:pwelch 输出的功率谱能直观看到信号能量落在哪个频段。如果发现 50 Hz 附近有个很强的尖峰,那是工频干扰,需要在带通之外额外加一个陷波滤波器;如果低频段能量很高,说明运动伪影没处理干净。
参数说明:chName 只是标签,不影响计算,但会影响结果可读性。多通道场景下建议把通道标签和电极位置对应关系维护在单独的文件里,避免换了电极布局后画出来的图对不上肌肉。
4.3 运行 main_demo 并解读结果曲线
配置没问题后,命令行直接运行:
main_demo逻辑说明:脚本会生成 result_single_channel.fig 和 activation_summary.txt 两个文件。fig 是交互式曲线图,包含原始 sEMG、带通滤波后信号、线性包络和最终激活度四个子图。txt 是在无图形环境下快速看数字的文本表,记录峰值激活度、平均激活度和达到峰值的时刻。
参数说明:如果 activation_summary.txt 里峰值超过 1,先别急着改算法,回去检查 MVC 标定是否偏低。激活度大于 100% 在数据层面几乎都是分母问题,不是分子太强。
4.4 失败时的第一现场:查看中间变量
结果不对时,第一步不是改算法,而是看中间变量。命令行里做两件事可以快速定位:
loaded = load(fullfile('data', 'subject01_task.mat')); plot(loaded.data.channels(1,:));逻辑说明:画出原始信号后,先确认有肉眼可见的发力段和静息段。如果原始波形平得像直线,说明数据本身有问题,可能是文件选错或者通道索引不对,算法怎么调都没用。再画 envTask 看包络是否在发力段有明显抬高,包络没起伏就去查 compute_mvc 的分母。
参数说明:排查 sEMG 代码讲究方向,永远从原始数据往输出方向查,不要反着猜滤波器参数。我习惯每次调试都先 whos 确认变量维度,再 plot 一个通道的原始波形,反复确认后再碰滤波参数。
4.5 参数调优的通用顺序
如果跑通了但结果不好看,调参顺序也很重要。先调 bandpassFreq 让滤波后噪声明显下降,再调 lpFreq 控制包络光滑度,最后才动 A 系数。很多人一上来就调 A,结果包络本身是脏的,再好的非线性映射也拉不回来。
5. 避坑记录:滤波、归一化与数据装载的 5 个常见问题
这类肌电代码跑不通,大多不是原理问题,而是细节约定没对上。把调试过程中反复踩的坑整理成五条,每一条按现象、原因、解决记录,照方抓药就能解决大部分问题。
5.1 激活度曲线全是 0 或接近 1 的直线
现象:预处理输出有波形,但最终激活度要么恒为 0,要么恒为 1。 原因:MVC 基准算错。整段 MVC 文件做平均会把基准拉高,任务段激活度整体趋近 0;反过来,MVC 文件里被试没真正发力,基准太小,任务段包络一除就超 1,被默认裁剪后全变 1。 解决:回到 compute_mvc 用滑动窗口求峰值,窗口 0.5 秒。同时采集 MVC 时给出明确发力提示音,确保在窗口内达到真实最大发力。我一般把 mvcVal 打印出来,再对比静息段包络均值,MVC 基准至少应该是静息包络的三倍,否则标定失效。
5.2 报错“未定义函数或变量 butter”
现象:跑到 preproc_sEMG 第一步就中断,提示找不到 butter。 原因:MATLAB 没有安装 Signal Processing Toolbox,或者安装后许可证没有激活。 解决:命令行执行 ver('signal') 检查工具箱版本。没有工具箱就改用 designfilt 设计 IIR 滤波器;有工具箱但提示未知函数,说明当前许可证环境没包含该模块,重启 MATLAB 并重新激活即可。顺便说一句,高阶 butter 在极端截止频率下会出现数值病态,必要时换成两个二阶滤波器串联。
5.3 滤波后首尾抖动厉害
现象:线性包络开头和结尾出现大幅振荡,中间段正常。 原因:filtfilt 做零相位滤波时会在首尾做填充,边界效应被高阶滤波器放大。 解决:在 preproc_sEMG 里对 raw 做首尾扩展,用 padarray 配 symmetric 方式扩展,扩展长度取滤波器阶数的两倍。滤波完成后再裁剪回原始长度。这样处理后的包络首尾平滑,不会影响激活度起始段的分析。
5.4 激活度输出超过 100%
现象:activation_summary 里峰值激活度达到 1.3 甚至更高,曲线明显越过 1。 原因:任务段的噪声或伪影高于 MVC 段,常见于电极在动态动作中松动、线缆摆动产生尖峰,也可能 MVC 采集时没有达到真正最大发力。 解决:先对原始波形做中值滤波剔除孤立尖峰,再进入预处理链路。同时回到 data 目录确认 MVC 文件和任务文件来自同一电极位置、同一放大倍数。换一组干净数据就能判断是算法问题还是数据采集问题。
5.5 load 进来的变量名和代码对不上
现象:main_demo 读取 data 字段时报错,但用 fieldnames 看明明有数据。 原因:.mat 文件保存时用了别的结构体名,不同批次数据混入不同命名。 解决:在 main_demo 开头加一层适配,用 fieldnames 获取实际字段名再统一映射到 data。以后接入新数据不需要改核心算法,只改映射层。从那以后我每次换数据集,都强制先 whos 确认字段名,再放行主流程。
5.6 兜底排查顺序
遇到难缠问题时,有一个通用检查顺序可以兜底:先确认数据维度是不是通道 × 采样点;再画原始波形看有没有明显噪声;然后检查滤波后基线是否归零;最后打印 MVC 基准和任务段峰值。按照这个顺序走一遍,大多数问题都能定位到具体环节。
6. 改造方向:从单通道离线走向多通道与在线估算
这套代码的定位是离线单通道处理,但把它改造成自己的工具,有三条很实用的路径。
第一条是批量多通道处理。在 main_demo 外面套一层 for 循环,逐个通道调用 preproc_sEMG 和 muscle_activation,结果存成结构数组。注意每个通道的 MVC 必须用自己的 MVC 包络,不能交叉使用。同时把输出文件的命名加上通道名,避免覆盖。
第二条是接入实时估算。把整段处理改成滑动窗口,每次取过去 0.25 秒的数据调用一次预处理和激活计算,就变成一个实时激活度指示器,可以用在生物反馈场景。关键是把滤波器状态持久化,不要每帧重新创建滤波器对象,否则会引入大量计算开销。Matlab 的 filter 函数支持初始条件参数,第一次初始化后,后续帧都沿用上次的滤波器状态。
第三条是验证改造是否引入偏差。用合成信号做测试:生成一个已知包络形状的 sEMG 仿真信号,把包络当作目标激活度,跑完看算法还原的激活度与目标的 RMS 误差。我一般把这个验证写成独立脚本,每次改完算法都跑一遍,曲线对不上就说明某一步改坏了。
从做第一个版本到现在两年多,我每次接手别人写的肌电代码,都强制走一遍“原始数据直方图 → 滤波后基线检查 → MVC 标定位置可视化 → 激活度范围检查”这条链路,宁可多花十分钟也绝不跳过任何一个中间环节。希望帮到你。
本文还有配套的精品资源,点击获取