Matlab实现sEMG信号处理：从原始波形到肌肉激活度%MVC曲线

简介：肌肉激活度计算是表面肌电信号处理中的常见需求，这份代码与数据包面向生物力学、康复工程或人机交互领域的研究者和学生，解决从原始肌电信号到激活度曲线的完整求解问题。压缩包共5个文件，包含主程序与激活度计算函数两个脚本、动作期肌电数据与最大自主收缩数据两个数据文件，以及一个文本说明文档，整包仅135KB，轻量易用。运行主程序后可直接输出肌肉激活度曲线，省去自行采集与预处理的步骤；数据文件覆盖常态动作与最大自主收缩两种场景，便于对照理解归一化、滤波和激活度建模等关键环节。目前已有3505人学习下载，适合需要复现算法、开展表面肌电实验或用于课程设计与论文写作的初学者快速上手。

1. 别在原始波形里硬找规律：这份代码把毛刺信号变成规则的激活度曲线

只要碰过表面肌电信号（sEMG），你就知道原始波形有多“脏”：又密又抖，里面混着工频干扰、运动伪迹和放大器噪声，直接盯屏幕根本判断不了肌肉到底用了多大力。这份资源就是一套把脏数据变成干净指标的完整工具链——从原始 sEMG 一路处理到肌肉激活度曲线，全部用 Matlab 实现，并且自带一组示例数据，下载后改一下路径就能跑通。做运动生物力学分析、康复评估、人机工效研究，或者在做肌电假肢、外骨骼控制策略验证的人，都能直接拿它当计算底座，不用再从零搭信号处理流程。

2. 项目定位与技术选型：一套从原始 sEMG 到 %MVC 的完整信号链

2.1 线性包络法：肌肉激活度不是凭空来的，是一步步处理出来的

肌肉激活度在sEMG分析里的地位，相当于“标准化后的肌肉发力程度”。最常用的算法叫线性包络法，链路大致是：原始信号 → 带通滤波 → 去工频干扰 → 全波整流 → 低通滤波 → 归一化。每一步都有明确目的，带通滤波是为了把 sEMG 有效频带（通常 20~500 Hz）从基线漂移和高频噪声里剥离出来；全波整流是把负半周翻到正半周，让信号全部落到能量侧；低通滤波则是把这些密集的尖峰“抹平”，得到一条随时间变化的光滑包络线；最后除以一次最大自主收缩（MVC）的包络值，就得到了以 %MVC 为单位的激活度。

这个链路看着简单，但每一步的参数都直接影响最终曲线形态。比如低通截止频率取 3 Hz 和取 10 Hz，包络线的平滑程度完全不同，后者会更“毛糙”但能保留更多爆发力细节。代码包里把这条链路固定成了可复现的脚本，你不用关心中间每个拐点怎么选，也能先跑通再按自己场景调。我从项目实测来看，这种“先给完整链路、再开放参数修改”的做法，比一堆孤立的滤波函数更实用。

2.2 技术栈与目录结构：Matlab 基础环境就能跑，不依赖第三方运行动态库

整个项目是纯 Matlab 脚本，依赖集中在 Signal Processing Toolbox 上，用到的核心函数是butter、filtfilt、iirnotch，这些都是信号处理工具箱里的标配。数据文件我用的是 .mat 格式存储，因为 Matlab 对 .mat 的读取效率远高于 txt 和 xlsx，也避免了 Excel 里小数精度截断的问题。如果你手头是设备导出的 csv，代码里也预留了readtable的读取分支，只要按列整理成“每列一个通道”就行。

代码包内部结构如下，拿到手不要急着到处点，先按这个脉络找文件：

muscle_activation/ ├── data/ │ ├── raw_emg.mat % 示例数据：一组静息+收缩状态的原始 sEMG │ ├── mvc_ref.mat % 示例数据：MVC 实验段信号，用于归一化基准 │ └── data_readme.txt % 数据采集条件、通道顺序、采样率说明 ├── code/ │ ├── main_activation.m % 主脚本，整条流水线在这里串起来 │ ├── compute_envelope.m % 封装的包络计算函数，滤波+整流+低通 │ ├── remove_artifacts.m % 基线漂移与运动伪迹剔除辅助函数 │ └── plot_activation.m % 画图脚本，输出对比图 └── results/ % 运行后生成的激活度曲线图

compute_envelope.m是全项目复用率最高的函数，所有通道、所有受试者的数据处理都会调它。目录里把数据、代码、结果分开，这是处理科研数据集的基本习惯，二次开发时你替换数据文件不会污染代码逻辑，改完参数重新跑也方便对比。建议你先打开data_readme.txt看一遍示例数据的采样率和通道顺序，再进代码，避免后面数据对不上号。

3. 核心代码拆解：从原始信号到激活度曲线的完整链路

3.1 主脚本结构：一次跑完滤波、整流、包络、归一化

主脚本main_activation.m控制整条流程，我拆成几个模块来看。第一段负责清空环境、读数据、建时间轴，代码开头长这样：

% main_activation.m % 将一段原始 sEMG 转换为肌肉激活度曲线（线性包络 + %MVC 归一化） clc; clear; close all; %% 1. 读取数据 % data/raw_emg.mat 中约定存放变量 raw_emg（double 数组）和 fs（采样率） S = load('data/raw_emg.mat', 'raw_emg', 'fs'); raw_emg = S.raw_emg; fs = S.fs; % 如果拿到的是 csv 文件，用下面这段替代 load % T = readtable('data/raw_emg.csv'); % raw_emg = T.signal_ch1; % fs = 1000; t = (0:length(raw_emg)-1) / fs; % 时间轴，单位秒

这段做了两件事：一是把数据从磁盘读进内存，二是把时间轴算出来，后面画图要用。fs是采样率，sEMG 设备常见的取值是 1000 Hz 或 2000 Hz，示例数据用的是 1000 Hz。raw_emg是一维 double 数组，如果导师给你的数据是多通道矩阵，那每列是一个通道，主脚本里要加一句raw_emg = raw_emg(:, 1);先取目标通道。这里不做任何信号处理，先把数据形态理顺，避免后面矩阵运算报维度错误。

3.2 滤波器设计：带通、工频陷波、低通三件套的配合

滤波是整个项目的核心，也是最容易出“玄学”问题的地方。代码里先设计 4 阶 Butterworth 带通滤波器，再用filtfilt做零相位滤波。零相位很重要，普通filter会引入相位延迟，导致滤波后的信号在时间轴上整体平移，激活度峰值位置全错。

%% 2. 带通滤波 % sEMG 有效能量集中在 20~500 Hz [b_bp, a_bp] = butter(4, [20 500] / (fs/2), 'bandpass'); emg_filt = filtfilt(b_bp, a_bp, raw_emg); %% 3. 工频陷波 % 国内电网 50 Hz，北美或部分设备是 60 Hz，按实际情况改 wo = 50 / (fs/2); % 陷波中心频率归一化 bw = 3 / (fs/2); % 陷波带宽，取 3 Hz 足够窄 [b_notch, a_notch] = iirnotch(wo, bw); emg_clean = filtfilt(b_notch, a_notch, emg_filt);

带通的截止频率写的是 20 和 500 Hz，这是处理 sEMG 最常用的区间：20 Hz 以下基本是运动伪迹和基线漂移，500 Hz 以上是放大器高频噪声。如果你在跑动态动作（比如跑步、跳跃），运动伪迹能量会更强，我一般会把低端从 20 抬高到 30，牺牲一点点有效信号来换干净的基线。wo和bw是陷波器的两个关键参数，分别控制中心频率和带宽。带宽取太大会伤到邻近频段的有效信号，取太小又压不住工频波动，3 Hz 是折中值。如果你的信号里 50 Hz 干扰不重，也可以注释掉陷波这一段直接往下走。

整流和低通包络在封装函数里实现。看compute_envelope.m内部：

%% 4. 全波整流 + 低通包络 function envelope = compute_envelope(emg, fs, lp_cutoff) % 默认低通截止频率 6 Hz if nargin < 3 lp_cutoff = 6; end rect = abs(emg); % 全波整流 [b_lp, a_lp] = butter(4, lp_cutoff / (fs/2), 'low'); envelope = filtfilt(b_lp, a_lp, rect); % 零相位低通 end

整流用abs一行搞定，数学上等价于把负半周折叠到正半周。低通截止频率lp_cutoff是包络线平滑度的控制旋钮：取 3 Hz，包络线很平缓，适合观察缓慢的等长收缩；取 8~10 Hz，包络线能跟上快速发力变化，适合爆发力动作。这里的nargin < 3是 Matlab 的函数参数缺省写法，调用方不传第三个参数时自动用 6 Hz。示例数据在静息和慢速收缩场景下采集，6 Hz 刚好不多不少。

3.3 归一化：用 MVC 段做基准，得到可横纵向对比的 %MVC

归一化这步直接决定激活度数值有没有可比性。代码里单独读了 MVC 参考段，走同一条包络链路再取均值：

%% 5. 归一化：除以 MVC 基准 % 同样先算 MVC 段的包络，取收缩最稳定的一段做均值 M = load('data/mvc_ref.mat', 'mvc_raw'); mvc_env = compute_envelope(M.mvc_raw, fs, 6); mvc_base = mean(mvc_env(fs*1 : fs*2)); % 取收缩后第 1~2 秒稳定段 activation = envelope ./ mvc_base * 100; % 结果单位是 %MVC

mvc_base的选取有个细节：不是把整段 MVC 信号全平均，而是取发力稳定后的第 1 到第 2 秒。开头那一下波形还没顶上去，平均进去会拉低基准，后面所有激活度都会被高估。除以mvc_base之后乘 100，得到的是百分数，便于后面的统计分析和画图标注。activation就是最终肌肉激活度曲线，每一帧对应时间轴t上的一个值。

画图脚本plot_activation.m会输出三张图：原始信号、滤波后信号、激活度曲线，方便一眼看到处理前后的差异。你用这个项目跑完，重点看第三张图——激活度曲线是否干净、有无异常尖峰，有尖峰说明滤波链路里某个环节的参数没匹配上数据特征，回头调参数而不是去改画图代码。

4. 避坑与常见问题排查：滤波相位、工频干扰和归一化翻车实录

4.1 激活度曲线带周期性毛刺：工频干扰没滤干净

现象：跑出来的激活度曲线每隔 0.02 秒（50 Hz 周期）出现一次小鼓包，整条曲线看起来像锯齿叠加在平滑线上。原因：陷波滤波器带宽太窄（比如取了 1 Hz），无法完全抑制频率波动的工频干扰；或者陷波那一步被整体注释掉了。解决：把陷波带宽从 1 Hz 扩大到 3 Hz，先看滤波后信号的频谱里 50 Hz 处有没有明显凸起再决定是否收窄。如果数据采集时用了差分电极和右腿驱动电路，工频干扰本身很小，跳过陷波也行，但前提是你确认过频谱，而不是默认“设备好就没干扰”。

4.2 激活度算出来超过 100%：归一化基准选错了

现象：某些收缩瞬间激活度飙到 150% 甚至 200%，曲线整体偏高。原因：最常见的是 MVC 段没有做同样参数的带通和低通处理，直接拿原始abs后的峰值当分母，导致分母偏小；或者 MVC 段本身采集时用力不充分，没有达到真正的最大自主收缩。解决：先确认mvc_ref.mat的采集条件是否真的包含最大发力，再检查compute_envelope是否对 MVC 段和普通段用了完全一致的参数。这类问题排到最后往往是数据采集环节的锅——受试者没有收到足够的发力激励，MVC 偏小，你只能重新采集。

4.3 滤波后波形整体往后平移：用了普通 filter 替代 filtfilt

现象：激活度曲线和原始信号在时间轴上错位，发力起始点明显滞后，爆发力峰值位置对不上。原因：filter是因果滤波，会引入与滤波器阶数相关的相位延迟，阶数越高延迟越大；而filtfilt是零相位滤波，先正向再反向处理，抵消了相位偏移。解决：全项目搜索filter(，凡是信号处理一律换成filtfilt。这条我踩过一次，当时对比两通道发力时序，差了几十毫秒，查了半天才发现是相位延迟造成的假性差异，换掉滤波方式后时序对上了。

4.4 读取 .mat 文件报错，提示 “Unable to load”

现象：load('data/raw_emg.mat', 'raw_emg', 'fs')直接报错，说不存在变量或文件。原因：代码包的 .mat 文件可能是从某实验室设备软件里直接另存的旧版本格式，或者是 v7.3 以上的 HDF5 格式，低版本 Matlab 不认；还有一种可能是文件路径里带了中文或空格，load解析失败。解决：先用whos -file data/raw_emg.mat查看文件里实际存了哪些变量名，确认变量名是raw_emg不是data_raw之类；再用matlab -v7在旧环境里重新保存一份。路径问题优先把整个项目放到纯英文目录下，避免文件路径里的中文干扰。

4.5 批处理多个受试者时结果对不上：参数被循环变量覆盖

现象：循环处理 10 个文件，前 5 个结果正常，后 5 个激活度普遍偏低，曲线形态也不对。原因：循环里某个滤波参数用了循环变量名，比如把fs误用作截止频率赋值，后续迭代把参数越改越偏；或者是filtfilt内部缓存了上一轮数据长度，导致边界效应污染。解决：每轮循环结束用assert(size(activation,1) == length(emg_filt))自检长度一致性，再把滤波器设计代码放到循环体外，只循环读写数据。核心原则是：滤波器系数只算一次，循环里只做数据切片和调用。

5. 运行方法与参数调整：数据格式约定和 Matlab 代码改造点

5.1 跑通步骤：从环境检查到出图

先把项目放到纯英文路径下，我习惯直接放D:\sEMG_activation\，避免中文路径触发 Matlab 老毛病。打开 Matlab，当前目录切到项目根目录，先确认装了 Signal Processing Toolbox，运行ver('signal')，如果是空输出就需要先装工具箱，这一步不解决，后面跑哪段都报错。

接着在命令行里执行主脚本：

run('code/main_activation.m')

顺利的话results/目录下会生成三张图。如果第一次跑就报错，优先看是不是data路径没写对——主脚本里用的是相对路径，当前目录必须在项目根目录，而不是code子目录。建议顺手把主脚本里cd到项目根目录的语句加上，或者统一改成用fileparts(mfilename('fullpath'))定位，这样不管从哪启动都不会跑丢路径。运行成功后，你先看激活度曲线的量级：静息段应该在 5% MVC 以下，收缩段能看到明显的爬升和回落曲线，如果静息段都快 20%，回去检查滤波参数或数据文件是否选错通道。

5.2 参数调整对照表：不同实验场景怎么改滤波值

不同动作类型的 sEMG 处理参数差异很大，代码包的默认值只适合慢速收缩场景，换成动态动作必须主动改。下面是我总结的参数对照：

带通低端调高的逻辑是压掉运动伪迹：动态实验里肢体摆动会产生大振幅低频干扰，不抬高低端，包络线会被伪迹顶起来，激活度虚高。低通截止频率调高则是为了让包络线跟上快速发力节奏。滤波器阶数要谨慎，阶数越高滚降越陡，但相位畸变也越明显，我用filtfilt时最多上到 6 阶，再高就会出现振铃现象，包络线出现不存在的波浪。改完参数一定要重新看一次原始信号和包络线的叠加图，确认包络没有超前或滞后于原始信号的能量变化。

5.3 用自己的数据替换示例：通道格式和采样率的适配方法

自己的数据是 csv 或 xlsx 时，先按列整理好，第一列是时间戳就删掉，Matlab 里直接用时间轴变量重建。主脚本里替换读取段：

T = readtable('data/my_emg.csv'); raw_emg = table2array(T(:, 1)); % 取第一列通道 fs = 1000; % 改成设备真实采样率

这里最坑的是通道顺序：有些设备导出的通道顺序和电极位置对应关系写在采集软件里，不写在数据里。先把data_readme.txt的通道顺序说明看一遍，确认第一列到底是哪块肌肉的信号，别把肱二头肌当肱三头肌分析。另外，如果设备采样率是 2000 Hz，但代码里fs还写 1000，滤波器的所有截止频率会直接翻倍错位，出来的曲线全是乱的。换数据源的第一件事就是核对采样率，这是我的血泪经验。

6. 二次开发方向：把激活度计算封装成函数与实时化改造

整体代码结构已经为二次开发留好了口子，compute_envelope本身就是独立函数。如果要做批处理，把主脚本里从滤波到归一化的过程封装成calc_activation(raw_emg, fs, mcv_raw, params)，参数用 struct 传，循环读文件往里丢就行。我一般会把所有参数收敛到一个struct里，像params.lp_cutoff = 6，避免函数参数一多就传错位。

想要更贴合生理模型，可以在线性包络之后接一个非线性映射层。经典做法是 Lloyd 和 Besier 的激活模型，把线性包络通过一个带指数系数的传递函数转换成激活度，这样包络值和激活度之间不是简单的线性比例，更接近肌肉的力-电关系。代码库里没有现成实现，但改造点很明确：在compute_envelope返回之后，加一行activation = (exp(A * envelope) - 1) / (exp(A) - 1)，其中 A 是肌肉非线性系数，取值范围 -3 到 0，试几次就能找到适合自己数据的值。

如果想做实时计算方向，把filtfilt换成filter配合延迟补偿，再改成逐帧滑动窗口处理。实时场景下不能用filtfilt，因为它要等整段数据到齐才能回卷，延迟不可接受。我常用的是dsp.FilteredAudio或手写一个环形缓冲区，每 100 ms 推一次新窗口进去算包络，延迟控制在 200 ms 以内，足够做实时视觉反馈。从那以后我每次处理 sEMG 数据，都先把采样率、工频频率、动作类型三个信息写进文件名再动代码，强制走一遍确认流程，才不会再掉进参数错位的坑。希望帮到你。

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