简介:面向需要处理表面肌电信号(sEMG)的科研与工程人员,这份MATLAB源码针对肌电信号幅值差异大、难以直接比较的问题,实现了归一化处理,并将归一化前后的波形以图形方式清晰显示,方便用户观察信号整体变化。适合新手及有一定经验的开发人员快速掌握肌电信号的预处理与绘图方法,也可作为生物医学信号处理课程设计或毕业设计的实用参考。资源包大小仅770B,共1个文件,核心为一个MATLAB脚本,代码经过测试校正,百分百可运行,涵盖了从表面肌电信号读取、归一化参数计算、数值映射到最终对比图绘制的完整流程,步骤清晰,即使初学者也能轻松理解关键逻辑。用户还可在此基础上灵活调整归一化方式(如极值归一化、均值方差标准化)或添加滤波、特征提取等处理,便于二次开发与算法比较,资源由达摩老生出品,质量有保障。目前已有1721人学习下载,适用于康复评估、人机交互、运动分析等场景,是一份小而精的肌电信号处理入门资源。
1. 表面肌电归一化不是可选项:为什么幅值必须折算到同一尺度
表面肌电信号(sEMG)的原始幅值,只要电极位置偏移 1 cm、贴合松紧不同,或者皮下脂肪厚度有差异,同一个受试者在相同力量水平下测到的电压就可能相差一倍以上。直接把两条 sEMG 曲线叠在一起谈“谁的激活更强”是没有意义的,除非先折算到同一个参考尺度上。归一化处理解决的就是这个问题:把实验或握力任务中采集到的原始表面肌电信号,通过基准收缩或信号自身幅值区间,转换为可比较的相对值。这个资源里的emg_jizhangli.m就是围绕“先消除噪声与漂移、再做 min-max 归一化、再用统一坐标显示”这条线完成处理的。适合做实验数据预处理、康复评估以及需要批处理肌电数据的 MATLAB 使用者,不管现在是在校做课题还是在写仪器采集程序,这套处理思路都能直接移植到自己的信号链路上。
2. MVC 归一化与 min-max 归一化:方法选择和信号预处理
2.1 不同归一化方法的实际含义
表面肌电归一化在生理信号处理里常用两种参考:最大自主收缩(maximum voluntary contraction, MVC)和信号自身范围。MVC 归一化需要让受试者先做几次持续 3 秒左右的最大力量收缩,取稳定段的 RMS 或平均整流值作为 100% 基准,之后所有信号都表示成“MVC%”。它的优点是结果有明确的生理含义,适合比较不同肌肉、不同受试者的激活程度。缺点是需要额外采集 MVC 数据,并且受试者如果不会正确发力,MVC 值容易偏低,后续所有归一化值都会超过 100%,图形直接失效。
min-max 归一化则直接以当前这一段信号的幅值最小值作为 0、最大值作为 1,计算简单,适合同一段信号内部相对变化的分析。emg_jizhangli.m中比较稳妥的做法是:先用带通滤波和陷波把噪声清掉,再做 min-max 归一化,这样可以不需要额外测量 MVC,又能把握力过程中从静息到最大强度的变化范围压缩到 0-1 之间。要注意的是,min-max 结果会被偶发的尖峰伪迹拉偏,所以滤波和平滑步骤的质量,直接决定归一化曲线是否可信。
2.2 滤波参数怎么取:带通 20-500 Hz 与 50 Hz 陷波
表面肌电的有效能量主要分布在 20-500 Hz。低于 20 Hz 的部分主要是运动伪迹和基线漂移,高于 500 Hz 的部分大多是测量电路的高频噪声。常见做法是先用零相位带通滤波器处理一次,再用窄带陷波器把 50 Hz 工频干扰去掉。MATLAB 中filtfilt比filter更合适,因为它是零相位滤波,不会造成波形整体延迟,多通道叠加显示时也不会出现通道之间的时间错位。
fs = 2000; % 采样率,根据采集设备修改 f_low = 20; % 高通截止频率,单位 Hz f_high = 500; % 低通截止频率,单位 Hz [b, a] = butter(4, [f_low/(fs/2), f_high/(fs/2)], 'bandpass'); emg_filtered = filtfilt(b, a, raw_emg); % 零相位带通滤波执行完这段代码后,emg_filtered中仍可能残留 50 Hz 工频成分,尤其是台式采集设备和未屏蔽线缆。再用iirnotch把 50 Hz 及附近 3 Hz 带宽内的成分衰减掉:
wo = 50 / (fs/2); bw = 3 / (fs/2); [b_notch, a_notch] = iirnotch(wo, bw); emg_clean = filtfilt(b_notch, a_notch, emg_filtered);参数上有几点经验:四阶巴特沃斯在大多数 sEMG 场景下够用,阶数太高会在通带边缘引入振铃;如果采样率只有 1000 Hz,f_high应降到 400 Hz,因为 500 Hz 已经接近奈奎斯特频率的一半,滤波器过渡带会变得不可控。陷波带宽 3 Hz 是折中值,太窄收敛慢,太宽会把 48-52 Hz 附近的真实肌电一起削掉。做完滤波后再执行emg_clean - mean(emg_clean)去掉直流偏置,后续归一化才不会被静息基线抬高。
2.3 归一化之前是否要整流和平滑
严格来说,min-max 归一化既可以直接作用于滤波后的双极性信号,也可以作用于整流后的信号。直接作用于双极性信号的问题在于,负半周的幅值会让最小值接近负的峰值,归一化结果会呈现以 0.5 为中心的剧烈振荡,看不出肌肉收缩的包络。因此,更合理的流程是先把滤波后的信号取绝对值,再做平滑,得到反映肌肉激活程度的包络线。
平滑窗口通常选 20-100 ms。窗口太短包络仍然抖动,窗口太长会丢掉动作电位发放的细节。对于握力这类力量变化较慢的动作,我一般用 50 ms 窗口,对应代码是movmean(rect_signal, round(0.05*fs))。如果最终目的是做频谱分析,就不要整流和包络;如果目的是做图形显示和力量等级对比,就一定要做整流和平滑。这里不要把平滑和低通滤波混为一谈,movmean是滑动平均,能用很直观的方式把包络拉出来。
| 步骤 | 常用参数 | 作用 | 注意事项 |
|---|---|---|---|
| 带通滤波 | 20-500 Hz,4 阶 | 去除运动伪迹与高频噪声 | 采样率低时降低 f_high |
| 陷波滤波 | 50 Hz,3 Hz 带宽 | 去除工频干扰 | 确认当地电网频率为 50 Hz |
| 去直流 | mean 或 detrend | 消除基线偏移 | 静息段较短时优先 detrend |
| 整流 | abs | 把负相位翻正 | 用于包络显示 |
| 平滑 | 50 ms 移动平均或 RMS | 得到激活包络 | 不要用于频谱分析 |
3. 实现 emg_jizhangli.m:从原始 EMG 到归一化曲线
3.1 数据读取与变量命名
资源中的emg_jizhangli.m是一个可直接运行的 MATLAB 脚本,它的核心数据流可以抽象为加载、滤波、归一化、绘图四个环节。假设数据文件和脚本放在同一目录,常见做法是通过load命令把信号读入工作区:
raw_data = load('grip_emg.mat'); raw_emg = raw_data.emg; % 假设数据文件中变量名为 emg fs = 2000;load返回一个结构体,访问字段名要和grip_emg.mat中保存的变量名一致。如果数据是 CSV 格式,改用readmatrix('grip_emg.csv'),第一列是时间,第二列是肌电。读者如果拿到的是没带数据文件的纯脚本,也可以先用一段模拟信号验证归一化流程:
t = 0:1/fs:5; % 生成 5 秒时间轴 raw_emg = 0.45*sin(2*pi*3*t) + 0.2*randn(size(t));这里的目的是把关注点放在归一化本身,而不是卡在数据格式上。
3.2 滤波和归一化的完整代码
下面是一段与脚本思路一致的参考实现,把上一章的预处理步骤串联起来。实际采集时,放大器的直流偏置可能比较高,运动伪迹也会带来缓慢漂移,所以先取静息段的前 0.5 秒平均值并减去,保证 min 值反映的是真实基底噪声:
raw_emg = raw_emg - mean(raw_emg(1:round(0.5*fs))); % 去直流 [b, a] = butter(4, [20/(fs/2), 500/(fs/2)], 'bandpass'); emg_f = filtfilt(b, a, raw_emg); [bn, an] = iirnotch(50/(fs/2), 3/(fs/2)); emg_f = filtfilt(bn, an, emg_f); rect_emg = abs(emg_f); smooth_emg = movmean(rect_emg, round(0.05*fs)); min_val = min(smooth_emg); max_val = max(smooth_emg); normalized_emg = (smooth_emg - min_val) / (max_val - min_val + eps);这段代码的关键是分母中的eps。当一段信号完全静息、包络几乎为零时,max_val - min_val可能接近 0,加入eps可以避免出现Inf或NaN。处理完成后,normalized_emg被限制在 0 到 1 之间,0 表示当前段的最低激活水平,1 表示本次任务中的最大激活水平。
3.3 为什么不能用原始信号直接除最大值
有些教程会把归一化写成norm = raw_emg / max(abs(raw_emg))。这种写法在双极性信号上会出现两个问题:一是正峰是 1,负峰是 -1,图形看上去是幅度调制而不是激活程度;二是基线漂移会让最大值偏高,整条曲线被压低,力量差异看不出来。
即使先取abs再除以最大值,也只是做了等比例缩放,并没有把基线放到 0,静息阶段的噪声仍然会在图上占据较大视觉比例。因此在emg_jizhangli.m中,先整流、平滑,再通过 min 和 max 做线性映射,才能得到一条从 0 附近开始、峰值接近 1 的激活曲线。不同受试者之间的差异,此时会更直观地反映在曲线的上升斜率、平台宽度和回落陡度上,而不是绝对幅值。
| 平滑窗口 | 适用场景 | 说明 |
|---|---|---|
| 20 ms | 快速挥动、力量突变 | 包络抖动大,峰值偏高 |
| 50 ms | 握力、等长收缩 | 包络稳定,适合演示 |
| 100 ms | 慢速力量变化 | 趋势清楚,但细节丢失 |
4. 绘图显示归一化结果:坐标轴、对比与伪影识别
4.1 基础绘图与坐标轴设置
归一化完成后的图形显示,不是简单一句plot就结束。表面肌电归一化曲线的横轴是时间,纵轴刻度在 0 到 1 之间;纵轴留白过大会让曲线显得很平,留白过窄又会让曲线顶部被截断。我会用下面这种方式一次画出两行子图,上一行是带通和陷波后的原始肌电,下一行是同一时段的归一化包络,方便核对信号对齐情况:
t = (0:length(normalized_emg)-1) / fs; figure('Color', 'w'); subplot(2,1,1); plot(t, emg_f, 'Color', [0.30, 0.30, 0.30]); ylabel('原始肌电 (mV)'); xlim([0, t(end)]); grid on; subplot(2,1,2); plot(t, normalized_emg, 'LineWidth', 1.2, 'Color', [0.85, 0.33, 0.10]); ylim([-0.05, 1.05]); ylabel('归一化幅值'); xlabel('时间 (s)'); grid on;ylim([-0.05, 1.05])留出 5% 边距,既能完整看到 0 和 1,又不会让曲线贴在边框上。上一行的原始肌电不要用强饱和色,波形密集时会变成一片黑色,[0.30, 0.30, 0.30]这样的深灰已经足够。
4.2 多通道和多段收缩的对比
实际握力实验经常同时采集多个通道,比如指浅屈肌、桡侧腕屈肌、肱桡肌。把每个通道的normalized_emg放到同一张图里,配合图例,可以快速看出发力时各肌肉的启动顺序:
channels = {ch1_normalized, ch2_normalized, ch3_normalized}; labels = {'指浅屈肌', '桡侧腕屈肌', '肱桡肌'}; colors = lines(3); figure('Color', 'w'); hold on; for k = 1:length(channels) plot(t, channels{k}, 'LineWidth', 1, 'Color', colors(k, :)); end hold off; legend(labels, 'Location', 'northeast'); ylim([-0.05, 1.05]); ylabel('归一化幅值'); xlabel('时间 (s)');这里有一个容易忽视的细节:不同通道的 min 和 max 是各自独立计算的,因此两个通道即使原始幅度差 3 倍,归一化后也都能到 1,不能据此判断哪块肌肉更“强”。要比较肌肉之间的强度差异,需要所有通道共用一个参考值,比如后面提到的 MVC。如果只是观察启动顺序,各通道各自归一化是合理的;但描述“谁发力更大”时,就会产生误导。所以emg_jizhangli.m的纵轴只标“归一化幅值”,避免把相对值直接当成绝对值解释。
4.3 从图形上判断归一化是否成功
图形是质量检查的第一道关卡。如果归一化后的曲线在静息段不断出现窄尖峰,说明陷波带宽太窄或运动伪迹没有去干净,需要回头检查带通滤波器的f_low。如果曲线在最大发力处出现平顶,也就是有一段持续等于 1,说明该段信号已经发生过载削波,min-max 归一化的最大值不是真实肌电峰值,而是放大器饱和电压,必须回退到原始采集数据重新检查。
反过来,如果发现静息段整段都是 0,不是“无激活”的理想表现,通常是平滑窗口太短,静息噪声没有被纳入包络,min_val取到了接近零的极小值。这时把平滑窗口加大到 80 ms,曲线通常会变得更连续。这些判断做完之后,再保存图形才有意义。
| 绘图参数 | 推荐值 | 作用 |
|---|---|---|
| ylim | [-0.05, 1.05] | 给 0 和 1 留出边距 |
| LineWidth | 1 至 1.5 | 防止多通道重叠时看不清 |
| 网格 | on | 辅助比较时间点和峰值位置 |
| 图例位置 | northeast | 标注通道或实验状态 |
5. 验证 MVC 区间:自动定位最大自主收缩的小技巧
min-max 归一化虽然能把曲线压缩到 0-1,但当你想把结果解释为“接近最大自主收缩的百分比”时,只取整段信号的最大值并不可靠。如果最大收缩只是短暂尖峰,中等力量收缩会被显示成接近 100%;如果受试者没有真正发力,最大值又没有达到生理上限,归一化曲线整体会被抬高。更实用的做法是识别出收缩的稳定平台段,用平台段的平均值作为 100% 基准。
先从包络信号计算移动 RMS,找到最高平台的中心位置:
rms_win = sqrt(movmean(emg_clean.^2, fs)); % 1 秒窗口 RMS [~, idx] = max(rms_win); % 峰值位置 half = round(fs / 2); % 前后各取 0.5 秒 mvc_idx = max(1, idx-half) : min(length(rms_win), idx+half); mvc_level = mean(rms_win(mvc_idx)); % 作为 100% 基准得到mvc_level后,把所有时刻的短窗 RMS 转换成百分比:
rms_short = sqrt(movmean(emg_clean.^2, round(0.05*fs))); pct_mvc = rms_short / mvc_level * 100; plot(t, pct_mvc); ylabel('MVC %'); ylim([0, 120]);这里的参数选择和普通 min-max 有本质区别:MVC 用的是峰值附近的邻域平均值,而不是瞬时最大值。瞬时最大值叠加了噪声和电极运动干扰,峰值中心 0.5 秒的窗口才能代表“能维持住”的最大能力。两个移动窗口长度要分开:检测 MVC 用 1 秒窗口,能看到稳定力量平台;绘制百分比曲线用 50 ms 窗口,保留动作细节。
实际使用中容易遇到一个坑:当数据里有多次重复握力时,max(rms_win)会落在第一次或最后一次发力处,而第一次发力往往偏小,最后一次可能因疲劳而撑不住。这时只搜索整个时长的中间 80% 区域:
start_idx = round(0.1 * length(rms_win)); end_idx = round(0.9 * length(rms_win)); [~, idx] = max(rms_win(start_idx:end_idx)); idx = idx + start_idx - 1;这样识别出的 MVC 区间更稳定,再把它单独存成mvc_segment.mat,下一次对同一批数据进行归一化时直接调用,避免手工截取带来的批次差异。
本文还有配套的精品资源,点击获取