Matlab实战：表面肌电信号sEMG肌肉激活度计算与归一化

简介：这份资源面向运动生物力学、康复工程与人体运动分析方向的学习者和研究者，提供一套用Matlab求解肌肉激活度的完整实现方案。针对表面肌电信号（sEMG）处理中激活度难以量化的问题，资源给出了可直接运行的代码与配套实验数据，帮助读者从原始信号出发得到肌肉激活度曲线。压缩包共5个文件，约135KB，包含2个mat数据文件、2个m脚本文件和1个txt说明文件，分别对应最大自主收缩（MVC）与运动过程的sEMG数据、主程序及激活度计算函数，结构紧凑、便于快速上手。目前已有3507人学习下载，说明其在sEMG处理入门与教学场景中具有较高参考价值。读者可借助MVC归一化思路与激活度求解流程，理解从信号预处理到激活度输出的关键环节，并在此基础上替换自有数据开展实验分析。

1. 从一段肌电信号到一块肌肉的发力曲线：这套 Matlab 代码到底在算什么

如果你手头有一份表面肌电信号（sEMG）采集数据，却卡在“怎么把它变成能写进论文或报告的肌肉激活度曲线”这一步，那这套 Matlab 代码和数据就是冲这个场景来的。它做的事很具体：读入原始 sEMG 信号，做滤波、整流、包络提取，再通过归一化映射成 0 到 1 之间的肌肉激活度，最后把结果画出来、存下来。适合谁？做运动生物力学、康复评估、人机交互或者外骨骼控制的同学和工程师，尤其是那些不想从零推导公式、只想拿到一份能跑通、能改参数的参考实现的人。Matlab 在这个领域依然是主力工具，信号处理工具箱和绘图能力摆在那里，这套代码就是围绕这个工作流搭的。

2. 项目结构与技术选型：为什么是 Matlab 而不是 Python

2.1 目录里有什么：从数据到脚本的对应关系

拿到压缩包解压后，目录结构通常不会太复杂，但每一层都有明确分工。常见做法是分成三个区域：数据区、函数区、脚本区。数据区放.mat或.csv格式的原始 sEMG 记录，函数区放独立的处理函数，脚本区放主运行脚本和绘图脚本。我一般会先看主脚本，因为它是整个流程的入口，能快速判断数据怎么读、参数在哪改、结果往哪存。

一个典型的目录长这样：

muscle_activation_project/ ├── data/ │ ├── raw_emg_subject01.mat │ ├── raw_emg_subject02.mat │ └── mvc_reference.mat ├── functions/ │ ├── bandpass_filter.m │ ├── rectify_signal.m │ ├── envelope_extraction.m │ └── normalize_activation.m ├── main_activation.m ├── plot_results.m └── config_params.m

这个结构的好处是职责分离。config_params.m集中管理采样率、滤波截止频率、包络窗口长度这些容易变的参数，改的时候不用翻遍所有文件。functions/里的每个函数只做一件事，方便单独测试和替换。main_activation.m负责串联流程，plot_results.m只管画图。如果你要批量处理多个受试者的数据，只需要在main_activation.m里加一层循环，遍历data/目录下的文件即可。

2.2 技术选型：Matlab 在 sEMG 处理里的真实优势

有人会问，Python 有 SciPy、NumPy、Matplotlib，为什么还要用 Matlab？这个问题在工程现场很实际。Matlab 的优势不在语言本身，而在它的信号处理工具箱和交互式调试体验。butter、filtfilt、envelope这些函数开箱即用，文档齐全，参数含义清晰。更重要的是，Matlab 的绘图窗口可以实时缩放、查看数据点，调参的时候不用反复重新运行脚本。对于 sEMG 这种对滤波参数敏感的信号，这种交互式反馈能省很多时间。

另一个现实原因是协作。很多实验室和项目组的历史代码就是 Matlab 写的，导师和同事都在用，你换 Python 反而增加沟通成本。这套代码没有依赖任何付费工具箱之外的第三方库，核心函数都在信号处理工具箱里，移植和分享都方便。如果你确实需要 Python 版本，逻辑是一样的，只是把butter换成scipy.signal.butter，把filtfilt换成scipy.signal.filtfilt，包络提取用scipy.signal.hilbert或者滑动 RMS，后面会讲具体怎么对应。

2.3 数据格式：原始 sEMG 长什么样

原始 sEMG 数据通常是一个二维矩阵，行是采样点，列是通道。采样率常见的是 1000 Hz 或 2000 Hz，采集设备不同会有差异。数据里往往还包含一个时间向量，或者时间信息隐含在采样率里。MVC 参考数据是另一份记录，用来做归一化，通常是受试者在最大随意收缩时采集的几秒信号。

读数据的时候要注意一点：有些设备导出的.mat文件里变量名不统一，可能是emg_data、data、signal或者别的。主脚本里一般会有一行load然后直接取变量，如果你换了自己的数据，记得改这一行。我见过有人直接把数据拖进 Matlab 工作区，变量名对不上，脚本报错“未定义变量”，折腾半天。所以第一步永远是whos看一下工作区里到底有什么。

3. 核心处理链路：从原始信号到激活度的每一步

3.1 带通滤波：去掉运动伪迹和高频噪声

sEMG 的有效频段通常在 20 Hz 到 450 Hz 之间，具体取决于肌肉和电极位置。低于 20 Hz 的部分主要是运动伪迹和基线漂移，高于 450 Hz 的部分主要是电子噪声。带通滤波的目的就是保留有效频段，去掉这两头。代码里一般用 Butterworth 滤波器，因为它在通带内平坦，相位失真可控。

function filtered = bandpass_filter(raw, fs, low_cut, high_cut, order) % raw: 原始信号，列向量或矩阵，每列一个通道 % fs: 采样率，单位 Hz % low_cut: 低截止频率，通常 20 Hz % high_cut: 高截止频率，通常 450 Hz % order: 滤波器阶数，通常 4 阶 % 归一化频率，Nyquist 频率是 fs/2 Wn = [low_cut, high_cut] / (fs / 2); % 设计 Butterworth 带通滤波器 [b, a] = butter(order, Wn, 'bandpass'); % 使用 filtfilt 做零相位滤波，避免相位偏移 filtered = filtfilt(b, a, raw); end

这里的关键是filtfilt而不是filter。filter会引入相位延迟，不同频率的成分延迟不一样，导致波形变形。filtfilt做的是前向-反向两次滤波，相位延迟相互抵消，结果零相位。代价是计算量翻倍，但对于离线处理来说完全不是问题。参数方面，low_cut和high_cut要根据你的采集设备和肌肉类型调整。如果信号里运动伪迹特别重，可以把low_cut提到 30 Hz；如果发现信号能量损失太多，可以降到 15 Hz。order一般 4 阶够用，太高会导致数值不稳定，尤其是filtfilt两次滤波后等效阶数翻倍。

3.2 整流与包络提取：把振荡信号变成平滑曲线

滤波后的 sEMG 还是围绕零线振荡的，直接看幅值没有意义，因为正负抵消。整流就是取绝对值，把所有成分变成正向。然后做包络提取，也就是低通滤波或者滑动 RMS，得到一条平滑的曲线，反映肌肉激活的整体趋势。

function envelope = envelope_extraction(rectified, fs, window_ms) % rectified: 整流后的信号 % fs: 采样率 % window_ms: 滑动窗口长度，单位毫秒，通常 50-200 ms % 将窗口长度转换为采样点数 window_len = round(window_ms / 1000 * fs); % 使用滑动 RMS 提取包络 % 方法：对平方后的信号做滑动平均，再开方 squared = rectified .^ 2; kernel = ones(window_len, 1) / window_len; mean_squared = conv(squared, kernel, 'same'); envelope = sqrt(mean_squared); end

滑动 RMS 的物理意义是信号在窗口内的平均功率的平方根，比简单的低通滤波更能反映能量变化。窗口长度是个关键参数：太短，包络抖动大，看起来毛刺多；太长，包络过于平滑，丢失快速变化的细节。50 到 200 毫秒是常见范围，具体看你的动作有多快。如果是慢速等长收缩，200 毫秒没问题；如果是快速动态动作，50 到 100 毫秒更合适。代码里用conv做滑动平均，'same'保证输出长度和输入一致。注意边界效应，两端会有半个窗口的失真，如果在意的话可以裁剪掉。

3.3 归一化：把不同受试者的数据放到同一把尺子上

不同人、不同肌肉、不同电极位置采集到的 sEMG 幅值差异很大，直接比较没有意义。归一化就是解决这个问题的。最常见的方法是用 MVC 参考值做分母，把激活度映射到 0 到 1 之间。

function activation = normalize_activation(envelope, mvc_envelope) % envelope: 当前信号的包络 % mvc_envelope: MVC 参考信号的包络，取最大值作为参考 % 取 MVC 包络的最大值作为归一化基准 mvc_max = max(mvc_envelope); % 避免除以零 if mvc_max == 0 error('MVC 参考信号最大值为零，请检查数据'); end % 归一化到 0-1 activation = envelope / mvc_max; % 可选：截断到 0-1，防止异常值 activation(activation > 1) = 1; activation(activation < 0) = 0; end

这里有个细节：MVC 包络的最大值不一定出现在同一时刻，也不一定代表真正的最大激活。常见做法是取 MVC 信号包络的均值或者 95 百分位数，而不是最大值，因为最大值可能是噪声尖峰。代码里用max是最简单的版本，你可以改成prctile(mvc_envelope, 95)更稳健。另外，归一化后的激活度理论上应该在 0 到 1 之间，但实际信号可能有超出 1 的情况，截断一下更安全。

3.4 主脚本串联：一次跑通完整流程

把上面的函数串起来，主脚本的逻辑就很清晰了。读数据、滤波、整流、包络、归一化、绘图、保存，每一步对应一个函数调用。

% main_activation.m clear; clc; close all; % 加载配置参数 config_params; % 读取原始数据 load('data/raw_emg_subject01.mat'); % 变量名假设为 emg_raw load('data/mvc_reference.mat'); % 变量名假设为 mvc_raw % 带通滤波 emg_filtered = bandpass_filter(emg_raw, fs, low_cut, high_cut, filter_order); mvc_filtered = bandpass_filter(mvc_raw, fs, low_cut, high_cut, filter_order); % 整流 emg_rectified = abs(emg_filtered); mvc_rectified = abs(mvc_filtered); % 包络提取 emg_envelope = envelope_extraction(emg_rectified, fs, window_ms); mvc_envelope = envelope_extraction(mvc_rectified, fs, window_ms); % 归一化 activation = normalize_activation(emg_envelope, mvc_envelope); % 绘图 plot_results(activation, fs); % 保存结果 save('data/activation_result.mat', 'activation');

这个脚本跑通的前提是config_params.m里定义了fs、low_cut、high_cut、filter_order、window_ms这几个变量。如果你换了自己的数据，只需要改load那一行的文件名和变量名，以及config_params.m里的采样率。绘图函数plot_results一般会画两条曲线：原始信号和激活度，方便对比。保存结果用.mat格式，后续可以直接加载进其他脚本做统计分析。

4. 避坑与排查：跑这套代码时最容易翻车的几个地方

4.1 滤波后信号幅值异常小或者全是零

现象：滤波后的信号看起来几乎是一条直线，或者幅值比原始信号小了几个数量级。原因通常是归一化频率算错了。butter函数的Wn参数是相对于 Nyquist 频率的，不是相对于采样率。如果你把low_cut和high_cut直接除以fs而不是fs/2，截止频率就偏了一倍，可能把有效信号全滤掉了。解决方法是检查Wn = [low_cut, high_cut] / (fs / 2)这一行，确认fs是采样率而不是 Nyquist 频率。另外，如果high_cut超过了fs/2，butter会报错，因为数字滤波器无法处理超过 Nyquist 频率的成分。

4.2 包络曲线抖动严重，像没滤波一样

现象：包络提取后曲线还是毛刺很多，看不出平滑趋势。原因可能是窗口长度太短，或者整流前忘了做带通滤波。如果原始信号里的高频噪声没被滤掉，整流后这些噪声变成正向成分，滑动 RMS 也压不住。解决方法是先确认带通滤波生效了，然后增大window_ms，从 50 毫秒逐步加到 200 毫秒，看曲线什么时候变得平滑。但也不要加太大，否则快速动作的细节会被抹掉。另一个可能的原因是conv的'same'参数用错了，导致输出长度不对齐，画图时看起来错位。

4.3 归一化后激活度超过 1 或者全是 1

现象：归一化后的激活度曲线大部分值都在 1 附近，或者明显超过 1。原因通常是 MVC 参考值选得太小。如果 MVC 信号采集时受试者没有真正用全力，或者电极位置和测试时不一致，MVC 包络的最大值就会偏小，导致归一化时分母太小。解决方法是检查 MVC 数据的质量，看看包络最大值是否合理。如果 MVC 数据本身有问题，可以考虑用同一受试者同一肌肉在测试中的最大激活值作为参考，而不是单独的 MVC 采集。代码里加截断是一种补救，但根本问题还是参考值的选择。

4.4 换了自己的数据后脚本报“未定义变量”

现象：把load那一行改成自己的文件名后，脚本在下一行报错，说某个变量不存在。原因是.mat文件里的变量名和你代码里用的名字不一致。Matlab 的load会把文件里的变量原样加载到工作区，变量名是什么就是什么。解决方法是先用whos('-file', '你的文件名.mat')查看文件里有哪些变量，然后把代码里的变量名改成一致的。或者用data = load('你的文件名.mat');然后通过data.变量名访问，这样更明确。我一般会在脚本开头加一行whos('-file', data_file)作为检查，避免这种低级错误。

4.5 绘图时横坐标是采样点而不是时间

现象：画出来的图横坐标是 0 到几万，而不是 0 到几秒。原因是绘图时直接用了plot(activation)，没有指定时间向量。解决方法是构造时间向量t = (0:length(activation)-1) / fs;，然后plot(t, activation)。这样横坐标就是秒。如果数据里有时间戳，直接用时间戳更好。这个坑很小，但每次换数据都容易忘，尤其是从别人那里拿来的脚本，绘图部分往往写得很随意。

5. 二次开发与改造：把这套代码变成你自己的工具

5.1 批量处理多个受试者

原始代码一次只处理一个文件，但实际项目里往往有几十个受试者的数据。改造思路很简单：把主脚本里的文件加载和处理逻辑包进一个循环，遍历data/目录下所有.mat文件。用dir('data/*.mat')获取文件列表，然后逐个加载、处理、保存结果。保存的时候文件名加上受试者编号，避免覆盖。如果数据量大，可以在循环里加一个进度条waitbar，不然跑起来不知道进度。

files = dir('data/raw_emg_*.mat'); for i = 1:length(files) load(fullfile('data', files(i).name)); % ... 处理流程 ... save(fullfile('results', ['activation_', files(i).name]), 'activation'); end

5.2 参数扫描与敏感性分析

滤波截止频率和包络窗口长度对结果影响很大，但很多时候你不知道选多少合适。一个实用的改造是写一个参数扫描脚本，固定其他参数，只变一个，看激活度曲线的变化。比如window_ms从 50 到 200，步长 25，画在一张图上对比。这样你能直观看到窗口长度对平滑度和细节保留的权衡。滤波截止频率也可以这样扫。这个改造不需要改核心函数，只需要在外面加一层循环，把结果存到 cell 数组里，最后统一绘图。

5.3 从 Matlab 迁移到 Python 的对应关系

如果你需要把代码迁移到 Python，核心逻辑完全一样，只是函数名和调用方式不同。下面是一个对照表：

迁移时最容易出问题的地方是数组维度。Matlab 的conv对矩阵默认按列处理，Python 的numpy.convolve只支持一维数组，多通道需要循环或者用scipy.ndimage.uniform_filter1d。另外，Matlab 的索引从 1 开始，Python 从 0 开始，切片的时候要注意。

5.4 一个我踩过的坑：归一化基准的选择

早期我做肌肉激活度分析的时候，直接用了 MVC 包络的最大值做归一化，结果发现有些受试者的激活度曲线整体偏高，有些偏低，对比起来很别扭。后来才意识到，MVC 采集时受试者的发力程度、电极粘贴位置、皮肤阻抗都会影响最大值，单次采集的最大值波动很大。从那以后我每次做归一化之前，都会先检查 MVC 信号的包络曲线，确认最大值出现在合理的发力阶段，而不是噪声尖峰。如果 MVC 数据质量不好，我会改用测试数据中同一肌肉的最大激活值作为参考，或者用 MVC 包络的 95 百分位数代替最大值。这个习惯让我少了很多返工。希望这套代码和这些经验能帮到你，拿到数据后先跑通流程，再根据你的具体场景调参数、改逻辑。

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