简介面向表面肌电信号分析与肌肉激活度计算的Matlab实践资料适合生物力学、康复工程、人机交互及运动科学方向的学生和研究人员使用可用于从原始sEMG数据出发完成激活度估计与曲线绘制。压缩包共5个文件包含2个mat数据文件、2个m脚本文件和1个txt说明文件整体约135KB已提供普通sEMG数据与最大自主收缩MVC时的参考数据。代码以MVC归一化思路组织运行主脚本即可得到肌肉激活度曲线便于读者对照理解信号预处理、归一化与激活度求解流程。已有3509人学习下载适合课程实验、论文复现和算法入门时快速搭建可运行的Matlab示例并在此基础上调整滤波参数、数据分段和绘图细节形成自己的分析流程。1. 从一段肱二头肌的原始信号说起这套 MATLAB 代码到底能算出什么实验室里最常见的场景是这样受试者做完一组动作采集卡里丢下一堆.mat导师路过问一句「这几块肌肉的激活度是多少」。这份肌肉激活度 Matlab 代码及数据处理的就是这件事它把表面肌电信号 sEMG 从原始电压一路算成 0 到 100% 的激活度曲线。它适合刚接手肌电数据、想先把整条流程跑通再谈算法改进的人也适合要给运动学或力学数据配一条肌肉发力参考曲线的开发者。数据是现成的代码是可读的脚本不是封装好的黑匣子改参数的地方都摊在明面上。2. 目录结构与数据流从解压到激活度曲线的完整链路拿到一个压缩包最忌讳的操作是双击main.m直接 Run然后对着一屏红色报错发呆。正确的顺序是先看目录、再看数据格式、最后才动脚本。这一章先把结构讲清楚后面的代码讲解才有落脚点。2.1 解压后先看什么脚本职责与文件清单这类项目的目录通常不会太深一般拆成「入口脚本 函数库 数据 输出」四块。不同人整理的习惯不一样但职责划分基本逃不出下面这张表路径典型内容第一遍要不要动main_activation.m主入口串起读取、滤波、归一化、出图只改路径func/preprocess_emg.m带通 工频陷波返回滤波后信号一般不动func/emg_envelope.m整流 包络提取换 RMS 策略时动func/normalize_mvc.mMVC 归一化输出百分比换参考值就动这里data/raw/*.mat原始 sEMG通常带fs和通道名换成自己的data/mvc/*.mat最大随意收缩试次没有就得自己补采results/导出的csv与png自动生成真正需要你先确认的是数据文件的内部结构。.mat只是一个容器里面可能是一个double矩阵也可能是一个带data、fs、channelName字段的结构体甚至可能是别人用旧版本存的table。用whos -file先把变量名、类型、尺寸打出来比猜快得多% 先摸清 .mat 里存了什么别急着写滤波 info whos(-file, data/raw/trial01.mat); for k 1:numel(info) fprintf(%-12s %-10s %s\n, info(k).name, info(k).class, mat2str(info(k).size)); end S load(data/raw/trial01.mat); disp(fieldnames(S)); % 结构体字段一览没有 fs 就得去问采集参数这段代码做两件事whos(-file, ...)不加载数据就能读出变量清单适合大文件load之后看fieldnames确认采样率字段叫什么。参数上没有玄学info(k).size返回的是维度向量比如[60000 8]就是 60 秒、8 通道、采样率 1000 Hz 的常见形态。如果size是[8 60000]说明通道在行上后面所有按列滤波的调用都得先转置这一步不做滤波结果会完全错位——这种坑我见过不止一次。2.2 一次完整变换原始电压 → 滤波 → 包络 → 百分比从电压到激活度中间要过四道工序每一道都改变了信号的物理含义不能省略阶段目的关键参数输出量纲去均值去掉电极偏置带来的直流无mV零均值带通 20–450 Hz保留肌电主频段压掉运动伪迹和高频噪声4 阶零相位mV50 Hz 陷波压工频干扰及其谐波Q ≈ 30mV整流 包络把双极性振荡变成单边强度窗 50–100 msmVMVC 归一化消除个体与电极差异MVC 峰值%这里有个容易被忽略的权衡包络的平滑程度直接决定时间分辨率。窗口取 50 ms能看清快速发力时的上升沿但曲线抖取 200 ms趋势漂亮可是一次 100 ms 的短促收缩会被抹平。做步态或者吞咽这类快节奏动作我一般从 50 ms 起步做等长收缩的疲劳分析100–150 ms 更稳。写论文时这个参数必须交代清楚否则别人复现不出你的曲线形状。3. 核心算法拆解滤波、整流、包络与 MVC 归一化这一章把每个环节的代码摊开讲包括为什么这么选参数、换参数会带来什么后果。整套流程本身不复杂难的是每个参数背后的取舍。3.1 带通与陷波先把 50 Hz 工频和基线漂移清掉表面肌电的有效带宽大致在 20–450 Hz低于 20 Hz 的成分主要是电极与皮肤的机械漂移、运动伪迹高于 450 Hz 多半是放大器噪声。工频干扰在 50 Hz 及其谐波上幅度可能比肌电本身还大不陷波的话包络会被抬起来一个台阶。%% 带通 陷波先把无关频段清掉 fs 1000; % 采样率必须和采集设备一致 fc [20 450]; % 带通上下限单位 Hz wn fc / (fs/2); % 归一化到 Nyquist 频率 [b,a] butter(4, wn, bandpass); % 4 阶 Butterworth [sos,g] tf2sos(b, a); % 转二阶节数值更稳 emg double(emg(:)); % 强制列向量防止通道方向搞反 emg emg - mean(emg); % 先去直流否则 filtfilt 边缘瞬态很大 emg_bp filtfilt(sos, g, emg); % 零相位滤波不引入群延迟 % 50 Hz 工频陷波Q 越大陷波越窄 wo 50 / (fs/2); [bw, aw] iirnotch(wo, 30); emg_bp filtfilt(bw, aw, emg_bp);逻辑上分三步先归一化频率再设计滤波器最后做零相位滤波。filtfilt是关键它正反各滤一次相位响应抵消代价是边缘会有瞬态所以后面出图时首尾各裁掉 0.5 秒比较稳妥。butter的阶数不要贪高4 阶滚降已经够用8 阶以上再配合filtfilt容易在低频段出现数值不稳定真要用高阶就统一走sos形式。iirnotch的第二个参数是品质因数取 30 对应带宽约 1.7 Hz够窄不会把 50 Hz 附近的肌电成分一起吃掉如果工频漂移明显可以把 Q 降到 20 放宽一点。采样率这一项必须亲口确认。有人拿着 2000 Hz 采的数据按 1000 Hz 算归一化频率全部翻倍滤波器实际通带跑到 40–900 Hz结果就是曲线看着有模有样数值却对不上任何文献。3.2 整流加包络滑动 RMS 与低通滤波两条路肌电是零均值的双极性信号直接平均会归零所以要先整流再取包络。工程上有两条常见路线效果接近取舍不同。第一条是滑动 RMS物理意义明确就是窗口内的均方根对应肌肉发力水平%% 路线一滑动 RMS emg_rec abs(emg_bp); % 全波整流 win round(0.05 * fs); % 50 ms 窗口 emg_rms sqrt(movmean(emg_rec.^2, win)); % 平方-均值-开方movmean是矩形窗等效于一个低通但旁瓣抑制差曲线会带一点台阶感。想要更平滑可以把movmean换成movmedian去尖峰或者用高斯核做加权平均。第二条是整流后低通滤波本质是把包络当作低频信号滤出来%% 路线二整流 低通包络 emg_rec abs(emg_bp); [b2, a2] butter(4, 6/(fs/2), low); % 截止 6 Hz emg_env filtfilt(b2, a2, emg_rec); emg_env max(emg_env, 0); % 滤完可能出现极小负值削掉两条路的参数可以互相换算50 ms 矩形窗的等效噪声带宽约 10 Hz对应低通截止大概在 4–6 Hz 区间。截止频率越高包络越贴近原始起伏但也越毛躁降到 3 Hz 以下曲线会明显滞后于真实发力时刻。我一般先用 6 Hz 看一眼原始包络形状再决定要不要降到 4 Hz 做正式结果。提示包络之后再接任何非线性处理比如阈值分割找发力起始点都要考虑包络本身带来的延迟别把延迟算进生理反应时间里。3.3 MVC 归一化激活度百分比是怎么落到 0~100 的到这一步信号还是 mV不同人、不同电极位置、不同皮肤阻抗下没有可比性。归一化到最大随意收缩MVC是肌电领域的通行做法让受试者对该肌肉做最大力收缩 3–5 秒取这段的包络峰值作为参考其他试次除以它。function ref mvc_reference(mvcFile, cfg) %MVC_REFERENCE 从一次 MVC 试次里提取参考包络峰值 % mvcFile MVC 数据文件路径 % cfg 含 fs / ch / sos / g 的配置结构体 S load(mvcFile); x double(S.data(:, cfg.ch)); x x - mean(x); x filtfilt(cfg.sos, cfg.g, x); % 与正式数据同一套滤波 env sqrt(movmean(abs(x).^2, round(0.05*cfg.fs))); % 取最高的 0.2 s 求均值避免单个尖峰把参考值抬高 n round(0.2 * cfg.fs); ref mean(maxk(env, n)); end这里的细节比代码本身重要。第一MVC 试次必须走和正式数据完全相同的滤波与包络参数否则参考值和分子不在一个尺度上。第二取峰值不要用max用一个短窗的均值抗尖峰maxk拿前 0.2 秒再平均是常见做法。第三MVC 一般重复 2–3 次取各次参考值的平均单次收缩受被试配合度影响很大。最容易翻车的地方是用试次自身的最大值做归一化。那样每条曲线都会被拉到 100%看起来漂亮实际上把发力强弱的信息全抹掉了横向比较彻底失效。4. 跑通实操从加载数据到出图算法讲完接下来把它拼成能跑的东西。这一章的脚本可以直接抄改两个路径就能出图。4.1 环境与路径三行代码先把工程挂起来MATLAB R2018b 及以上版本即可butter、filtfilt、iirnotch需要 Signal Processing Toolboxmovmean、maxk是基础函数。路径不要写死绝对路径用mfilename拿脚本所在目录做根换电脑、换盘符都不用改代码clear; clc; close all; root fileparts(mfilename(fullpath)); % 以脚本所在目录为根 addpath(genpath(fullfile(root, func))); % 递归加入函数库 dataDir fullfile(root, data, raw); resDir fullfile(root, results); if ~exist(resDir, dir); mkdir(resDir); end cfg struct( ... fs, 1000, ... % 采样率 Hz务必和设备一致 band, [20 450], ... % 带通范围 Hz notch, 50, ... % 工频 Hz envWin, 0.05, ... % 包络窗口 s ch, 1); % 当前处理通道cfg集中放参数是后面批量处理和二次开发的基础所有魔法数字都从这里取改一处全篇生效。addpath(genpath(...))会把子目录一起挂上函数多了也不用逐个addpath。4.2 单试次跑通一段脚本走完全流程下面这段把前面的环节串起来从读文件到出图一次走完%% 单试次全流程 S load(fullfile(dataDir, trial01.mat)); emg double(S.data(:, cfg.ch)); emg emg(:) - mean(emg); % 列向量 去均值 % 1) 带通 [b, a] butter(4, cfg.band/(cfg.fs/2), bandpass); [sos, g] tf2sos(b, a); emg_bp filtfilt(sos, g, emg); % 2) 工频陷波 [bw, aw] iirnotch(cfg.notch/(cfg.fs/2), 30); emg_bp filtfilt(bw, aw, emg_bp); % 3) 整流 滑动 RMS 包络 env sqrt(movmean(abs(emg_bp).^2, round(cfg.envWin*cfg.fs))); % 4) MVC 归一化 ref mvc_reference(fullfile(root,data,mvc,mvc01.mat), cfg); act 100 * env / ref; act min(act, 100); % 削顶防止超过 100% % 5) 出图 t (0:numel(act)-1) / cfg.fs; figure(Color, w); subplot(2,1,1); plot(t, emg_bp); ylabel(幅值 (mV)); title(sprintf(通道 %d 滤波后 sEMG, cfg.ch)); grid on; subplot(2,1,2); plot(t, act, LineWidth, 1.2); ylabel(激活度 (%)); xlabel(时间 (s)); ylim([0 110]); grid on; saveas(gcf, fullfile(resDir, trial01_activation.png));五段各自独立任何一段出问题都能单独拎出来看。第 4 步的min(act, 100)是后悔药MVC 采集不充分时正式试次的包络可能超过参考值不削顶会画出 130% 这种没法解释的数。第 5 步的ylim([0 110])固定纵轴多条曲线对比时不会因为自动缩放造成误判。跑通之后先别急着批量把emg_bp单独画出来看一眼如果幅值集中在 ±0.02 mV 这种量级多半是单位问题原始数据可能是 V如果波形整体贴着一条斜线说明去均值之后还有漂移带通下限要往上抬。4.3 批量处理与结果落盘一个试次跑通剩下的就是循环。批量处理的关键是把每个试次的结果整理成表方便后面导进统计软件files dir(fullfile(dataDir, *.mat)); T table(); for k 1:numel(files) S load(fullfile(dataDir, files(k).name)); emg double(S.data(:, cfg.ch)); emg emg(:) - mean(emg); [b,a] butter(4, cfg.band/(cfg.fs/2), bandpass); [sos,g] tf2sos(b,a); x filtfilt(sos, g, emg); [bw,aw] iirnotch(cfg.notch/(cfg.fs/2), 30); x filtfilt(bw, aw, x); env sqrt(movmean(abs(x).^2, round(cfg.envWin*cfg.fs))); ref mvc_reference(fullfile(root,data,mvc,mvc01.mat), cfg); act min(100 * env / ref, 100); T.name(k,1) string(files(k).name); T.meanAct(k,1) mean(act); T.peakAct(k,1) max(act); T.durS(k,1) numel(act) / cfg.fs; fprintf([%d/%d] %s 峰值 %.1f%%\n, k, numel(files), files(k).name, T.peakAct(k,1)); end writetable(T, fullfile(resDir, activation_summary.csv));fprintf的进度输出不是装饰批量处理最怕中途静默失败有了每行的峰值就能立刻发现某个文件是不是算出了 0.0%。writetable输出 CSV 而不是 MAT是为了让不做 MATLAB 的合作者也能打开。注意ref在循环里重复计算了实际用的时候提到循环外这里保留是为了让每一步都自洽。5. 常见问题排查曲线发毛、数值爆表、对不上文献时看这几处跑不通的时候报错反而是好事最怕的是不报错但结果不对。下面五条是我自己和身边人踩过的按「现象 → 原因 → 解决」记下来。5.1 信号层面的三类异常现象一激活度曲线是一条几乎为零的直线但原始信号看起来有起伏。原因是量纲或滤波参数不对。一种常见情况是原始数据单位是 V包络算出来是 0.0001 量级除以同样是 V 的 MVC 参考后数值正常但画图时纵轴被自动缩放看起来像直线另一种是fs填错导致归一化频率越界butter直接报错或者通带跑到别处。解决方法是先fprintf(std%.4f\n, std(emg))打印原始标准差确认数量级再核对fs和whos里的采样点数与时长是否自洽。现象二曲线毛刺很重像心电图一样一格一格的。包络窗口太短或者陷波没做干净。50 ms 窗口在 1000 Hz 采样下只有 50 个点movmean的矩形窗旁瓣抑制差起伏会原样透出来。把cfg.envWin调到 0.1或者改用低通 4 Hz 的包络路线毛刺基本能压下去。如果毛刺是等间隔的先看频谱上 50 Hz 有没有残留Q 值调到 20 试试。现象三曲线首尾各有一小段剧烈跳变中间正常。这是filtfilt的边缘瞬态几乎所有零相位滤波都有。4 阶带通在 1000 Hz 下瞬态长度大约是几十到上百毫秒。解决办法是分析时裁掉首尾各 0.5 秒或者在滤波前用padarray镜像延拓再滤滤完截回原长度。裁剪是最省事的代价是丢掉一点数据做短时动作分析时要在实验设计阶段就预留这段时间。5.2 数值与流程层面的两类异常现象四所有试次的激活度峰值都接近 100%完全分不出强弱。归一化用错了参考。要么用了试次自身的最大值要么 MVC 文件里混进了静息段。正确的做法是每个肌肉单独采 MVC参考值取多次收缩峰值的平均且 MVC 数据必须走完全相同的滤波链路。检查方法很简单把ref和每个试次的max(env)都打印出来如果两者比值恒为 1说明参考就是自身。现象五和文献里的激活度曲线对不上形状相似但时间轴错位几十毫秒。采样率不一致或者没做机电延迟补偿。肌电从肌肉兴奋到产生力矩有 30–100 ms 的延迟如果你的激活度要和力矩、关节角对齐需要把激活度曲线后移这段时间。另一种情况是肌电和运动学数据采样率不同比如 1000 Hz 对 100 Hz直接按索引对齐会累积误差正确做法是resample到同一时间基再用时间戳插值。注意五条里最容易误判的是第一条。看到直线先去怀疑算法往往查了半天才发现是数据单位或者采样率的问题先把输入摸清楚再动滤波器。6. 二次开发换成自己的采集设备、加通道、接运动学数据换成自己采的数据只需要对齐三件事采样率、物理单位、通道顺序。这三项任一对不上后面的滤波和归一化全是白算。加通道最简单的做法是把整个流程封成一个单通道函数输入列向量输出激活度外层用循环按列扫一遍而不是在脚本里复制粘贴八份。function act emg2act(emg, fs, mvcRef, cfg) %EMG2ACT 单通道 sEMG 转激活度百分比方便批量与二次封装 emg double(emg(:)) - mean(emg); % 去直流 emg filtfilt(cfg.sos, cfg.g, emg); % 带通参数外部传入 emg filtfilt(cfg.bn, cfg.an, emg); % 陷波 env sqrt(movmean(abs(emg).^2, round(cfg.envWin*fs))); act min(100 * env / mvcRef, 100); % 归一化并削顶 end函数化的好处是可测。改完参数后拿一段已知静息段跑一遍激活度应该贴近 0拿 MVC 段跑一遍峰值应该落在 95–100 之间。这两个断言过了说明链路没断。接运动学数据时把激活度resample到 100 Hz再按机电延迟后移 50 ms 左右比在原始采样率上硬对齐稳得多。改算法参数时有个习惯我一直保持每次只动一个动完立刻回看静息段和 MVC 段两条基准曲线。曾经为了压毛刺把包络窗口一路加到 300 ms结果曲线是好看了一次 150 ms 的短促发力直接消失白跑了一轮实验。从那以后我每次调envWin和低通截止都强制先把这两条基准曲线过一遍再去看正式数据。希望帮到你。本文还有配套的精品资源点击获取