企业官网建设流程全解析

简介：面向表面肌电信号分析与肌肉激活度计算的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 解压后先看什么：脚本职责与文件清单

这类项目的目录通常不会太深，一般拆成「入口脚本 + 函数库 + 数据 + 输出」四块。不同人整理的习惯不一样，但职责划分基本逃不出下面这张表：

真正需要你先确认的是数据文件的内部结构。.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 一次完整变换：原始电压 → 滤波 → 包络 → 百分比

从电压到激活度，中间要过四道工序，每一道都改变了信号的物理含义，不能省略：

这里有个容易被忽略的权衡：包络的平滑程度直接决定时间分辨率。窗口取 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 Toolbox，movmean、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和低通截止，都强制先把这两条基准曲线过一遍，再去看正式数据。希望帮到你。

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