简介这份资源围绕表面肌电信号sEMG的肌肉激活度求解展开面向运动生物力学、康复工程与人体工学研究方向的学习者和工程人员帮助其用 Matlab 完成从原始肌电采集数据到激活度曲线的完整计算。压缩包共 5 个文件包含 2 个 m 脚本、2 个 mat 数据文件与 1 个 txt 说明整体约 135KB脚本与数据配套齐全数据中既有常规动作的 sEMG 记录也提供最大自主收缩MVC状态下的参考数据说明文档用于交代调用方式。资源已有 3509 人学习下载被较多初学者用作入门练习。运行主程序即可直接得到肌肉激活度曲线读者可借此理解 MVC 归一化、滤波与包络计算的处理流程并可将脚本迁移到自己的实验数据上验证与修改适合作为信号处理课程练习或课题预处理的起步模板。1. 从一段跑不出包络的 sEMG 说起这份肌肉激活度代码到底解决什么拿到一份肌肉激活度 Matlab 代码最常见的翻车点不是跑不起来而是图跑出来了、曲线也平滑却没人说得清纵轴那个 0.37 到底代表什么。这套面向表面肌电信号 sEMG 的代码包做的事情其实很聚焦把毫伏级的原始电压序列一步一步压成 0 到 1 之间、能跨受试者、跨通道、跨天比较的激活度曲线。它适合刚接手 sEMG 数据、两周内要出一版能写进报告的结果的人也适合手上已有自己一套处理流程、想找份实现对照着查错的人。原始信号里混着工频、直流漂移、心电串扰和运动伪迹这些不处理干净后面所有统计都是空中楼阁。2. 四段信号链的参数怎么定带通 20–450 Hz 之后的每一步从电压到激活度是一条固定链路去直流 → 陷波 → 带通 → 全波整流 → 低通包络 → 归一化。这条链路上每个环节都有参数而参数错了未必表现为报错更常见的是曲线形状看着还挺像幅值却已经不可比了。这一章把每个环节的取值依据讲清楚后面拆代码时才知道每一行在防什么。2.1 采样率决定了带通上限别把 450 Hz 硬套在低采样率数据上sEMG 的有效能量主要集中在 20–450 Hz 这个区间低于 20 Hz 的部分大量是运动伪迹和电极-皮肤界面的慢漂移高于 450 Hz 的部分能量已经很低再往上留只会把高频噪声一起放进来。但这条 20–450 Hz 能不能用取决于采集时的采样率butter的归一化截止频率必须严格小于 1也就是 Nyquist 频率工程上通常要求上限不超过 0.45 倍采样率。常见采样率下能安全使用的带通上限采样率 fsNyquist建议带通上限说明1000 Hz500 Hz400 Hz450 已到 0.9 归一化数值不稳定2000 Hz1000 Hz450 Hz最常用的配置余量充足4000 Hz2000 Hz450 Hz上限不必跟着提高噪声带宽反而变大滤波器阶数也要克制。butter的 4 阶带通在数学上是 8 阶filtfilt正反各滤一次等效 16 阶过渡带够陡了再往上加振铃和边界效应会明显盖过收益。真正容易忽略的是长度约束filtfilt要求数据长度大于3*(max(length(b),length(a))-1)一段 0.5 s 的短 MVC 数据在 2000 Hz 采样下只有 1000 点4 阶带通系数长度是 9约束是 24 点还过得去但如果为了追求陡峭把阶数提到 6 阶、8 阶短段数据立刻报错。% 设计带通并检查归一化截止频率是否安全 fs 2000; % 采集时定死后续所有截止频率都依赖它 bp [20 450]; % sEMG 有效带宽 ord 4; % Butterworth 阶数 Wn bp / (fs/2); % 归一化到 Nyquist assert(all(Wn 0 Wn 1), 带通截止频率超出 Nyquist检查 fs 或 bp); [b_bp, a_bp] butter(ord, Wn, bandpass); % 打印频率响应肉眼确认 -3dB 点落在哪 figure; freqz(b_bp, a_bp, 4096, fs); title(带通频率响应); xlim([0 500]);这段代码的重点不在butter本身而在那行assert。数据来源五花八门有人给的是降采样到 500 Hz 的版本有人给的是原始 2000 Hz脚本里写死 450 Hz 上限跑到 500 Hz 数据上就直接抛Wn must be less than 1。把采样率作为参数传进来、先断言再设计比事后一行行 debug 省事得多。2.2 陷波和带通的先后50 Hz 工频在带内必须单独处理工频干扰落在 50 Hz部分采集环境是 60 Hz正好在 20–450 Hz 的带通范围内所以带通滤波器拦不住它必须单独做陷波。顺序上常规做法是先去掉直流分量detrend或直接减均值再陷波再带通。虽然线性时不变系统理论上可以交换顺序但数值实现上先处理大幅值的直流偏移再进窄带陷波边界瞬态会小一些。陷波的品质因数 Q 决定阻带宽度。Q 越大阻带越窄对 50 Hz 附近有用信号的误伤越小但对工频频率漂移的容忍度也越低Q 取 30–35 是常见折中对应阻带宽度大约 ±1.5 Hz。如果采集设备本身有硬件陷波软件这一步可以保留但把 Q 调大双重保险如果设备没有那这一步就不能省。% 50 Hz 陷波Q 值控制阻带宽度 f0 50; % 工频50 或 60 wo f0 / (fs/2); % 归一化中心频率 bw wo / 35; % 归一化带宽Q 35 [b_n, a_n] iirnotch(wo, bw); % 验证用 pwelch 对比陷波前后的功率谱 [pxx0, f] pwelch(x, hann(1024), 512, 2048, fs); x_n filtfilt(b_n, a_n, x); pxx1 pwelch(x_n, hann(1024), 512, 2048, fs); figure; plot(f, 10*log10(pxx0), f, 10*log10(pxx1)); legend(陷波前,陷波后); xlim([40 60]);这里bw用的是归一化带宽传入iirnotch后实际阻带约为f0/35Hz 的一半宽度写参数时按这个量级估就行。验证环节不要省有些采集系统导出的数据里工频已经被硬件压到很低再叠一层软件陷波反而会在 50 Hz 附近挖出一个凹陷如果后续要做频谱类指标比如中值频率这个凹陷会让结果失真。2.3 全波整流加低通截止频率决定了包络的平滑度和峰值保留整流就是取绝对值把双极性的电压信号变成单极性这样低通滤波后的结果才是包络而不是归零。接下来的低通截止频率是整个链路里最需要斟酌的参数定得太高包络跟着原始信号抖动激活度曲线全是毛刺定得太低上升沿被拉平峰值被削动作起始时间判不准。4–6 Hz 是文献和工程实践里都常见的取值。6 Hz 对应的时间常数大约是 26 ms对于大多数肢体动作肌肉激活上升时间在几十毫秒量级足够跟得上如果研究对象是慢肌主导的姿势维持任务可以降到 4 Hz 换取更平滑的曲线如果是快速发力类任务可以提到 8–10 Hz代价是曲线噪声变大。阶数用 2 阶就够阶数高了在阶跃处会产生明显振铃包络上会出现本不存在的过冲。% 全波整流 低通包络 lp 6; % 包络低通截止4~6 Hz 常用 [b_lp, a_lp] butter(2, lp/(fs/2), low); x_rect abs(x_clean); % 全波整流 env filtfilt(b_lp, a_lp, x_rect); % 零相移低通无延迟 env(env 0) 0; % 数值上可能出现极小负值兜底filtfilt在这里的作用是消除相位延迟——如果用filter单向滤波包络会比真实激活晚几十毫秒做多通道同步分析时这种延迟是致命的。至于最后那行把负值截断看起来多余实际上 2 阶低通在剧烈变化处会轻微过冲出现过 −0.001 这样的值画对数坐标或者算比值时会报错。2.4 归一化基准用 MVC 峰值还是整段最大值结论完全不同原始包络的单位是毫伏受电极阻抗、皮下脂肪厚度、贴片位置影响极大同一个人换一天贴电极幅值就可能差一倍。所以激活度必须归一化标准做法是除以最大随意收缩MVC时的包络峰值得到 %MVC。关键细节在峰值怎么取。直接对整段数据取max是最省事也最容易出问题的做法如果这段数据里恰好有一个动作伪迹或者电极被碰到那个尖峰会把整个归一化基准抬高之后所有值都被压扁。稳妥的做法是先对包络做 1 s 滑动平均再取最大值或者取 MVC 平稳段的均值而不是瞬态峰值。这两种做法在代码包里通常都留了开关跑之前确认用的是哪一种不然和别人的结果对不上。% 归一化先在 MVC 平稳段上取基准再作用于全部数据 mvc_seg env(idx_mvc_start : idx_mvc_end); % 人工标定的 MVC 区间 mvc_smooth movmean(mvc_seg, round(1.0*fs)); % 1 s 滑动平均压掉尖峰 mvc_ref max(mvc_smooth); act env / mvc_ref; % 得到 %MVC0~1 小数形式 act min(max(act, 0), 1); % 截断到 [0,1] act_pct act * 100; % 需要百分数时再乘 100movmean的窗口取 1 s 是个经验值太短压不住尖峰太长会把整个上升段抹平导致基准偏低。代码里idx_mvc_start、idx_mvc_end这种人工标定索引在批量处理时最好单独存一个标定表不要写死在脚本里——受试者一多改索引改到怀疑人生。3. 数据怎么喂进去三类文件格式与批量处理的目录组织代码包能不能直接用一半取决于算法一半取决于数据接口。sEMG 数据的来源太杂有的采集系统直接导出.mat有的是带表头的 CSV有的是第三方运动学格式里嵌着模拟通道。这一章讲怎么把不同来源的数据对齐成同一种结构以及多受试者、多通道时目录该怎么组织才能让一套脚本跑完全部数据。3.1 三类常见格式的读取方式与对应的坑数据来源典型文件读取方式常见坑采集系统导出 matraw.mat内含data、fs变量load变量名不统一采样率藏在嵌套结构体里文本或 CSV 表emg.csv首行是通道名readmatrix/readtable首行为非数值readmatrix会整行丢弃导致少一行第三方格式.c3d/.edf开源读取脚本单位是 V 不是 mV通道顺序与电极贴法不一致第二种格式的坑最隐蔽。readmatrix遇到首行是字符串时默认行为是跳过看起来一切正常但如果你按行号去对齐同步信号比如一个单独的光电标记文件整段数据就错位了一行而这种错位在图上根本看不出来只会在做事件相关平均时表现为莫名其妙的相位偏移。稳妥写法是先用detectImportOptions看清列类型再显式读。% 读 CSV 前先看结构避免首行被吞 opts detectImportOptions(emg.csv); opts.VariableNamesLine 1; % 表头在第 1 行 opts.DataLines [2 Inf]; % 数据从第 2 行开始 T readtable(emg.csv, opts); emg table2array(T(:, 2:end)); % 第 1 列通常是时间戳单独拿出来 t T{:, 1}; fs round(1 / median(diff(t))); % 从时间戳反推采样率比手填可靠最后那行从时间戳反推采样率是我很推荐的一个习惯。采集系统的配置界面上写的采样率和实际写进文件的时间戳偶尔会不一致尤其是中途改过配置的会话。用时间戳反推如果算出来是 1998 或者 2003至少你知道有问题直接写死 2000这个误差会被带进所有截止频率里。3.2 用一个 params 结构体收口所有可调参数脚本满天飞的项目里最怕的就是参数散落在十几个文件里改一个忘一个。代码包如果只有一个入口脚本通常会把参数集中在一个结构体里往下传。这种写法不优雅但极其实用跑新数据集时只改一处出问题回溯时只查一处。function P emg_params(fs) % 所有可调参数集中在这里函数入参 fs 保证一致性 P.fs fs; P.bp [20 min(450, 0.45*fs)]; % 上限自动适配低采样率 P.notch 50; P.q_notch 35; P.ord_bp 4; P.lp_env 6; P.ord_lp 2; P.mvc_win 1.0; % MVC 基准的平滑窗口秒 P.clip [0 1]; % 激活度截断范围 endP.bp那行的min(450, 0.45*fs)是整套参数里最值得抄的一句。它把 2.1 节讲的 Nyquist 约束固化进了参数生成逻辑换一套低采样率数据进来带通上限自动降到 0.45 倍采样率不用手动改也不会触发断言失败。3.3 批量处理的目录约定一个受试者一个文件夹多受试者处理最容易乱的地方不是算法是文件对应关系。比较省心的组织方式是按受试者 / 任务 / 文件三层放脚本用dir递归收集从路径里解析出受试者编号和任务名写进结果表里。这样任何一条结果都能追溯到原始文件出问题不用猜。root fullfile(pwd, data); subjDirs dir(fullfile(root, S*)); % 受试者目录以 S 开头 allAct table(); for i 1:numel(subjDirs) sid subjDirs(i).name; files dir(fullfile(root, sid, *.mat)); for j 1:numel(files) S load(fullfile(files(j).folder, files(j).name)); fs S.fs; P emg_params(fs); act process_emg(S.data, P); % 单文件处理函数见 4.1 % 结果按行追加保留来源信息 tmp table(repmat({sid}, size(act,1), 1), ... repmat({files(j).name}, size(act,1), 1), ... act, VariableNames, {subject,file,activation}); allAct [allAct; tmp]; end end writetable(allAct, activation_all.csv);两个细节dir返回的顺序在不同系统上不保证一致如果结果要和别的表按行合并先排序[allAct; tmp]这种逐次拼接在数据量大时会很慢几百个文件以内无所谓上千个文件建议先存 cell 最后vertcat。4. 逐段拆代码陷波、包络、归一化与出图的可复现写法前面讲的是参数依据这一章把完整处理函数拆开从原始矩阵进去到激活度矩阵出来中间每一步都能单独拿出来验证。拆代码的目的不是背下来而是知道哪一行出错时该去查什么。4.1 单通道处理函数六个步骤每一步都能单独打断点function act process_emg(x, P) % 输入 x: 单通道原始电压列向量 % 输出 act: 与 x 等长的激活度曲线0~1 x double(x(:)); x x - mean(x); % 1. 去直流 [b_n, a_n] iirnotch(P.notch/(P.fs/2), (P.notch/(P.fs/2))/P.q_notch); x filtfilt(b_n, a_n, x); % 2. 工频陷波 [b_bp, a_bp] butter(P.ord_bp, P.bp/(P.fs/2), bandpass); x filtfilt(b_bp, a_bp, x); % 3. 带通 20-450 Hz x abs(x); % 4. 全波整流 [b_lp, a_lp] butter(P.ord_lp, P.lp_env/(P.fs/2), low); env filtfilt(b_lp, a_lp, x); % 5. 低通取包络 env(env 0) 0; env_s movmean(env, round(P.mvc_win * P.fs)); % 6. 平滑后归一化 act env / max(env_s); act min(max(act, P.clip(1)), P.clip(2)); end这段代码里有两个容易被忽略的点。第一x double(x(:))把输入强行转成双精度列向量——有些采集系统导出的.mat里数据是int16原始码值如果不转换后续filtfilt会报类型不支持或者更糟在某些版本上静默按整数运算导致精度损失。第二归一化用的是max(env_s)1 s 平滑后的最大值而不是max(env)这两者在有尖峰的数据上结果能差 20% 以上而func内部做归一化意味着每一段数据都用自己当基准跨段比较的意义就被削弱了——真正严谨的做法是把 MVC 段的基准单独传进来改成act env / mvc_ref这里的写法适合单段探索性分析。阶数、截止频率这些参数传给butter时必须是归一化频率P.bp/(P.fs/2)这个除法一旦漏掉butter会把 450 当成归一化值直接报错如果误写成P.bp/P.fs不报错但滤波器完全跑偏20 Hz 高变成了 10 Hz 低曲线看着平滑了其实把有用信号全滤掉了。4.2 包络提取的三种写法什么时候该换方法实现优点代价整流 低通filtfilt(butter, abs(x))平滑、零相移、参数直观截止频率选择影响峰值保留滑动 RMSmovrms(x, win)有明确物理意义能量需要写窗口函数边界处理麻烦Hilbert 包络abs(hilbert(x))瞬时幅值分辨率高对噪声敏感宽带信号上包络抖动大默认用第一种就行它和大多数文献流程对得上复现别人的结果时争议最小。滑动 RMS 在需要和力信号做幅值-力关系拟合时更好用因为 RMS 本身就是能量度量窗口宽度直接对应时间分辨率。Hilbert 适合信号相对干净、需要看清瞬时变化的场景比如研究肌肉激活的快速调制但 sEMG 是宽带随机信号Hilbert 包络会保留大量高频起伏后面还得再低通一次等于绕了一圈。换方法时最容易犯的错是窗口长度和截止频率混着改。RMS 窗口 100 ms 在 2000 Hz 采样下是 200 点对应的等效低通截止大约 5 Hz和 6 Hz 低通量级接近如果同时又改窗口又改截止频率最后发现曲线不对根本分不清是哪个参数导致的。4.3 激活起始点检测阈值法够用但阈值怎么定有讲究很多分析要做事件相关平均需要把每次激活的起始时刻找出来。代码包里通常是阈值法包络超过基线的某个倍数比如静息均值加 3 倍标准差就记为起始低于同一阈值持续一段时间就记为结束。% 基于静息段的阈值法检测激活区间 rest env(1 : round(2*P.fs)); % 前 2 s 视为静息 th mean(rest) 3*std(rest); % 阈值 静息均值 3σ above env th; d diff([0; above; 0]); % 找上升沿和下降沿 onset find(d 1); offset find(d -1) - 1; % 剔除过短的激活段抖动引起的假激活 minLen round(0.2 * P.fs); % 最短 200 ms keep (offset - onset) minLen; onset onset(keep); offset offset(keep);阈值里的 3 倍标准差是个经验值。静息段如果本身就 noisy比如受试者紧张、电极接触不稳std(rest)会很大阈值被抬得过高小的激活就被漏掉反过来静息段太干净阈值过低基线波动会被判成激活。用diff找沿时[0; above; 0]首尾补零这一步不能少否则第一个采样点就是激活起点时会被漏掉。4.4 批量出图与导出图件和数值要能对上分析结果最终要落到两份东西上一份数值表给统计用一份图给报告用。图件的横轴必须是秒而不是采样点否则换个采样率图就变形纵轴如果是 %MVC要标清楚是相对于哪次 MVC。t (0:numel(act)-1) / P.fs; fig figure(Visible,off); plot(t, act*100, LineWidth, 1.2); hold on; yline(mean(act(rest_idx))*100, --); % 静息基线参考线 xlabel(Time (s)); ylabel(Activation (%MVC)); ylim([0 110]); grid on; exportgraphics(fig, fullfile(outDir, [sid _act.png]), Resolution, 300); close(fig); % 数值同步导出列名带单位说明 writetable(table(t, act*100, VariableNames, {time_s,activation_pctMVC}), ... fullfile(outDir, [sid _act.csv]));exportgraphics的 300 dpi 是投稿常用值Visible,off让批量出图时不会弹出一堆窗口把内存占满。数值表和图的横轴用同一个t生成保证两者严格对齐后期如果审稿要求补数据直接拿 CSV 就能对上。5. 避坑与排查五个把曲线画歪的具体问题5.1 幅值类异常曲线整体贴地或者整体顶到 1现象画出来的激活度曲线长期在 0.05 以下偶尔有小鼓包或者反过来整段都在 0.9 以上看不出起伏。原因归一化基准取错了。前者通常是mvc_ref取到了某个异常尖峰把分母放大了一个量级后者是基准取自静息段或者一段没发力的数据分母太小。解决把max(env_s)打印出来和静息段包络均值对比。正常情况下 MVC 峰值应该比静息高一个数量级以上具体倍数因肌肉和电极位置而异。如果差距不到 3 倍基本可以确定这段数据里没有真正的 MVC得回去找标定数据不要硬凑。另外确认截断之前的值分布min(max(act,0),1)会把超过 1 的部分压平掩盖掉基准偏低的问题。5.2 时域类异常包络比原始信号晚了一截现象把原始整流信号和包络叠在一张图上包络的上升沿明显滞后多通道之间还不同步。原因用了filter而不是filtfilt单向滤波引入了群延迟延迟量随截止频率变化通道间滤波器参数不同的话延迟还不一样。解决全链路统一用filtfilt。如果因为数据太短filtfilt报错优先考虑缩短滤波器阶数或者对数据两端做镜像延拓而不是退回filter。检查方法很简单喂一个理想阶跃信号进去看输出有没有延迟一眼就能看出来。5.3 边界类异常首尾各有一段明显的假波动现象激活度曲线在开头 0.2 s 和结尾 0.2 s 出现幅度异常的振荡形状像衰减的正弦。原因filtfilt的边界处理假设信号在端点之外是零实际数据端点不为零阶跃式的突变激发了滤波器的瞬态响应。阶数越高、截止频率越低这段瞬态越长。解决分析时把首尾各截掉 0.5 s 不要用于统计或者在滤波前对数据做镜像延拓处理完再裁回原长度。批量处理时我一般直接在流程里加一步裁边比事后逐个检查省事。5.4 通道类异常左右侧结果反了或者通道对不上现象同一受试者的左右腿激活度曲线互换或者某块肌肉的激活和动作明显不符。原因原始数据的通道排列顺序和电极贴法表不一致。不同采集系统的通道编号习惯不同有的是从 1 开始按贴片顺序有的是按放大器接口编号。解决每次处理前用一段已知动作的数据做单通道验证比如让受试者只抬左腿看哪个通道响应最大据此建立通道映射表写进配置而不是靠记忆。批处理脚本里加一条断言如果某通道在全段数据上的方差为 0说明该通道没接好直接报错而不是静默输出一条直线。5.5 可比性类异常同一个人不同天的结果对不上现象同一个受试者、同一块肌肉、同样的动作两天的激活度曲线形状相似但幅值差了一大截。原因电极重新贴过阻抗变了或者两次的 MVC 标定力度不一致也可能第二次数据里工频干扰更强陷波后残余噪声抬高了基线。解决跨天比较必须用当天的 MVC 做归一化不能跨天共用一个基准。同时记录每次的静息段包络均值和信噪比作为数据质量指标一并输出事后筛数据时有据可依。这部分血泪经验换来的一句话就是采集记录表比算法更影响结论。6. 验证激活度是否可信三个自检技巧与一条硬习惯拿到任何一份肌肉激活度代码在信任它之前先做三件事。第一件是合成信号反推。构造一段已知包络的调幅信号频率成分放在 100 Hz 附近包络用半波整流的低频信号采样率设为 2000 Hz。跑完整链路把输出和已知包络对比相关系数应该接近 1幅值比例应该稳定。fs 2000; t (0:1/fs:10-1/fs); trueEnv 0.5*(1 square(2*pi*0.5*t, 30))/2; % 已知的方波状包络 carrier sin(2*pi*100*t); x trueEnv .* carrier 0.01*randn(size(t)); % 叠一点噪声 P emg_params(fs); act process_emg(x, P); % 对比忽略两端边界比较形状和幅值 idx round(0.5*fs) : numel(t)-round(0.5*fs); fprintf(相关系数 %.3f\n, corr(act(idx), trueEnv(idx)));如果相关系数低于 0.9说明链路里有环节把包络改形了优先怀疑包络低通截止频率太低把方波的上升沿抹平或者陷波把载波削掉了。第二件是量纲校验。静息段的激活度应该在 0.02–0.1 之间也就是 2%–10% MVCMVC 段峰值接近 1。如果静息段就在 0.3 以上说明基准偏低或者基线没去干净如果 MVC 段只有 0.4说明基准取到了异常值。这个区间不是绝对标准但偏离太远一定有问题。第三件是抽查原始波形。随机挑三个文件把原始信号、整流信号、包络、激活度四条曲线画在同一张图上肉眼过一遍。自动化测试能查出报错查不出包络跟着心电在跳这种问题——心电串扰在胸腹部肌肉的 sEMG 上很常见表现为规律的尖峰一旦进到包络里激活度曲线就会出现和心跳同频的锯齿这种情况得回头做心电伪迹去除而不是调参数。校验项通过标准不通过时的第一排查点合成信号相关性 0.9包络低通截止频率静息段激活度2%–10% MVC去直流与归一化基准MVC 段峰值0.8–1.0MVC 区间标定原始波形肉眼检查无规律性尖峰心电串扰、电极接触从那以后我每处理一批新数据都会先跑一遍合成信号校验再挑三个文件手动画图看一遍最后才放到批处理里跑全量。多花二十分钟能省掉事后发现整批归一化基准错了、重跑两天的麻烦。希望帮到你。本文还有配套的精品资源点击获取
