sEMG分析

转载最新推荐文章于 2026-09-02 10:17:04 发布 5.2k 阅读 AI写同款

收录于

sEMG信号

lihe2021

码龄10年墨衍会员 · AI 创作全网分发

264粉丝81篇原创

关注

表面肌电信号（sEMG）介绍表面肌电信号（Surface Electromyography）是众多生物电信号中的一种，也是相对来说最容易获取的一种电生理信号。将电极片放置在人体的皮肤表面，可以记录到皮肤表面因肌肉收缩而产生的微弱的电位差，而这个微弱的电位信号通过肌电采集电路的放大和转换就形成了可用作处理的表面肌电信号。sEMG是及其微弱的，下面是肌电与常见电源的电压值比较。采集sEMG的装置很多，国内外均有比较成熟的产品在售，例如美国的Delsys、Noraxon、瑞典的MYON、德国FREEEMG（下图从左到右）国内阅读详情

AI摘要

该文详细介绍了肌电信号（sEMG）的产生、特点，以及基于sEMG的动作意图识别。在肌电信号采集过程中，需要考虑信号的微弱性、干扰因素，并通过滤波和陷波技术改善。随后，文章展示了sEMG的静态和动态显示方法，并讨论了预处理滤波去噪的重要性，包括使用IIR陷波器和数字滤波器。接着，提取了多种时域特征，如IAV、RMS等，用于特征选择。最后，提到了使用LDA、KNN等分类器进行动作识别的应用，如康复系统和下肢动作识别。

1肌电信号

肌电信号的产生

表面肌电信号是由多个运动单元发放的动作电位序列，在皮肤表面呈现的时间和空间上综合叠加的结果

在这里插入图片描述

肌电信号特点

（1）sEMG是一维时间动作电位序列

（2）sEMG是交流信号，幅值一般和肌肉运动力度成正比

（3）sEMG是一种非平稳的微电信号，其幅值在 0~1.5mv，有用信号频率位于 0~500HZ，主要能量集中在 20~150HZ

（4）SEMG一般比肢体运动超前 30~150ms产生，可以进行运动提前判断
在这里插入图片描述
基于sEMG的动作意图识别

2肌电信号采集

（1）sEMG是微弱电信号，易受到心电信号和工频噪声的干扰，需要很高的共模抑制比，用模拟滤波器和双T陷波器将信号滤波放大到0~5V
（2）将电极、放大滤波电路、陷波电路封装在一起，加屏蔽层保护，降低传导耦合的干扰
（3）采集方式：探针式、电极片式
（4）已有sEMG采集系统：美国Delsys生产的Trigno、加拿大Thought生产的Myoscan
在这里插入图片描述
肌电信号采集电路图

sEMG静态显示

clc
clear

addpath(genpath(pwd));   %% 将当前文件夹下的所有文件夹都包括进调用函数的目录
data=load('201801111534zhangaidong.mat');
ch1=data.data(:,1);
ch2=data.data(:,2);

% Hd=IIR_notchFilter(2000,50);
% ch1=filter(Hd,ch1);                   % 陷波后
% ch2=filter(Hd,ch2);
N=size(ch1,1);

figure(1);
plotHandlesEMG = zeros(2,1);
axesHandlesEMG = zeros(2,1);

for i = 1:2
    axesHandlesEMG(i) = subplot(2,1,i);   
    plotHandlesEMG(i) = plot(axesHandlesEMG(i),0,'-y','LineWidth',1);
    set(axesHandlesEMG(i),'XGrid','on','YGrid','on','XLimMode', 'manual','YLimMode', 'manual');
    set(axesHandlesEMG(i),'Color',[.15 .15 .15],'XLim', [0 N],'YLim', [-.0005 .0005]);
    title(sprintf('EMG %i', i)); 
end

set(plotHandlesEMG(1), 'Ydata', ch1); 
set(plotHandlesEMG(2), 'Ydata', ch2);

sEMG动态显示

function dynamicShow_6ch
addpath(genpath(pwd));   %% 将当前文件夹下的所有文件夹都包括进调用函数的目录
load('20170927_xinhongyu.mat');
global plotHandlesEMG flag ch1 ch2 ch3 ch4 ch5 ch6;
ch1 = p3_open_ch1;
ch2 = p3_open_ch2;
ch3 = p3_grasp_ch1;
ch4 = p3_grasp_ch2;
ch5 = p3_snooze_ch1;
ch6 = p3_snooze_ch2;

flag=1;
t = timer('Period', 0.1, 'ExecutionMode', 'fixedSpacing', 'TimerFcn', {@updatePlots, plotHandlesEMG});
figure('CloseRequestFcn', {@localCloseFigure,t});

plotHandlesEMG = zeros(6,1);
axesHandlesEMG = zeros(6,1);
for i = 1:6
    axesHandlesEMG(i) = subplot(3,2,i);   
    plotHandlesEMG(i) = plot(axesHandlesEMG(i),0,'b','LineWidth',1);
    set(axesHandlesEMG(i),'XGrid','on','YGrid','on','XLimMode', 'manual','YLimMode', 'manual');
    set(axesHandlesEMG(i),'XLim', [0 2000],'YLim', [-.00005 .00005]);
    title(sprintf('EMG %i', i)); 
end
drawnow
start(t);

function updatePlots(obj, Event,tmp)
global plotHandlesEMG flag ch1 ch2 ch3 ch4 ch5 ch6;
data_ch1 = ch1(flag:flag+2000);   
data_ch2 = ch2(flag:flag+2000);
data_ch3 = ch3(flag:flag+2000);   
data_ch4 = ch4(flag:flag+2000);
data_ch5 = ch5(flag:flag+2000);   
data_ch6 = ch6(flag:flag+2000);

flag=flag+200;
set(plotHandlesEMG(1), 'Ydata', data_ch1); 
set(plotHandlesEMG(2), 'Ydata', data_ch2); 
set(plotHandlesEMG(3), 'Ydata', data_ch3); 
set(plotHandlesEMG(4), 'Ydata', data_ch4); 
set(plotHandlesEMG(5), 'Ydata', data_ch5); 
set(plotHandlesEMG(6), 'Ydata', data_ch6); 
drawnow

function localCloseFigure(figureHandle,~,t)
stop(t);
delete(t);
delete(figureHandle);

3预处理滤波去噪

经由肌电采集系统得到的原始肌电信号如图所示，肌电采集的频率为2000HZ，EMG1表示一秒内肌肉的动作电位序列，EMG2表示一秒内肌肉静息时的电位序列。由肌电采集系统得到的肌电信号虽然经过了硬件上的滤波、去噪等一系列的处理，但是仍然存在工频噪声、尖峰幅值等干扰，现在对得到的原始肌电信号进行软件滤波处理。
在这里插入图片描述
对原始肌电信号进行FFT变换，得到信号的单边振幅谱如图所示，可以看到50HZ的工频干扰很大，现在采用IIR数字陷波器，得到陷波后的信号单边振幅谱如图所示，消除工频干扰后，肌电信号的主要频率集中在0-500HZ，尤其是在10-300HZ。
在这里插入图片描述
由于肌电信号的主要频率集中在10-300HZ，现在对肌电信号进行软件数字滤波，常用的数字滤波器有巴特沃兹滤波器和切比雪夫滤波器。

巴特沃兹滤波器比切比雪夫滤波器具有更平滑的滤波特性
切比雪夫滤波器在过渡带的衰减更快，虽然频率响应的幅频特性没有巴特沃兹滤波器的好，但是与理想滤波器频率响应曲线之间的误差小。
在这里插入图片描述
经过滤波后的肌电信号更加平滑，去除了高频尖峰的干扰。

4特征提取与选择

为了保证特征的连续性，通常采用设置时间窗＋增量窗的方式提取特征
在这里插入图片描述

时域特征

function feature = fARC(data)  
    % - ARC - Auto Regression Model Coefficients

    order = 6;              
    f = real(lpc(data,order)');    %real()取复数的实数部分 lpc():Linear prediction filter coefficients
    feature = -f(order+1,:);       %首项默认为1，从第二项开始
end

function feature = fIAV(data)
    % - IAV - Integrated Absolute Value
    feature = sum(abs(data))/length(data);
end

function feature = fLogD(data)
    feature = exp(mean(log(abs(data))));
end

function feature = fMAV(data)
    % - MAV - Mean Absolute Value
    feature = mean(abs(data));
end

function feature = fMAX(data)
    feature = max(abs(data));
end

function feature = fNZM(data)
    % - NZM - Non Zero Median 非0中位数
    data_nonzero = data(find(data ~= 0));
    feature = median(data_nonzero); 
end

function feature = fRMS(data)
    % - RMS - Root Mean Square
    feature = sqrt(mean(data.^2));  
end

function feature = fSSC(data)
    % - SSC - Slope Sign Change, number times
    DeadZone = 10e-7;
    data_size = length(data);
    feature = 0;

    if data_size == 0
        feature = 0;
    else
        for j=3:data_size
            difference1 = data(j-1) - data(j-2);
            difference2 = data(j-1) - data(j);
            Sign = difference1 * difference2;
            if Sign > 0
                if abs(difference1)>DeadZone || abs(difference2)>DeadZone
                    feature = feature + 1;
                end
            end
        end
        feature = feature/data_size;
    end
end

function feature = fVAR(data)
    feature = var(data);
end

function feature = fWA(data)
    % - WA - Willison Amplitude
    Threshold = 10e-7;
    data_size = length(data);
    feature = 0;
    if data_size == 0
        feature = 0;
    else
        for i=2:data_size
            difference = data(i) - data(i-1);
            if abs(difference)>Threshold
                feature = feature + 1;
            end
        end
        feature = feature/data_size;
    end
end

function feature = fWL(data)
    % - WL - Waveform Length
    feature = sum(abs(diff(data)))/length(data);
end

function feature = fZC(data)
    % - ZC - Zero Crossing
    DeadZone = 10e-7;
    data_size = length(data);
    feature = 0;

    if data_size == 0
        feature = 0;
    else
        for i=2:data_size
            difference = data(i) - data(i-1);
            multy      = data(i) * data(i-1);
            if abs(difference)>DeadZone && multy<0
                feature = feature + 1;
            end
        end
        feature = feature/data_size;
    end
end

频域特征

function feature = fFMD(data)
    % - MDF - Median Frequency
    % 未陷波滤波
    Fs = 2000;           
    feature = medfreq(data, Fs);    % medfreq()计算中值频率，2015a版本才有
end

function feature = fFME(data)
    % - MNF - Mean Frequency

    %未陷波滤波
    Fs = 2000; 
    feature = meanfreq(data, Fs);
end

function feature = fPSD(data)
    % PSD  Power Spectral Density 功率谱密度估计

   [b,a]=butter(2,[0.01,0.5]);           
   x= filter(b,a,data);

   Nfft     =  512;
   window   = hanning(256);
   noverlap = 128;
   Fs       = 2000;

   [pxx,f]=pwelch(x,window,noverlap,Nfft,Fs);
   feature = pxx;    %得到的是一组数据257维
end

在这里插入图片描述
特征提取

clear
clc

addpath(genpath(pwd));  

% 获取一维数据
load('mll.mat');
CH1   = double(mll_snooze_ch1);
CH2   = double(mll_snooze_ch2);
N = size(CH1,1);

%选择时间窗1秒（2000个点），增量窗0.1秒（200个点）
N = N/200;
for i=0:N-10
    ch1 = CH1(200*i+1:200*i+2000,:);
    ch2 = CH2(200*i+1:200*i+2000,:);

    % 提取常规时域特征,每个通道11个特征，总共22个特征
    A(1) = fIAV(ch1);   B(1) = fIAV(ch2);
    A(2) = fLogD(ch1);  B(2) = fLogD(ch2);
    A(3) = fMAV(ch1);   B(3) = fMAV(ch2);
    A(4) = fMAX(ch1);   B(4) = fMAX(ch2);
    A(5) = fNZM(ch1);   B(5) = fNZM(ch2);
    A(6) = fRMS(ch1);   B(6) = fRMS(ch2);
    A(7) = fSSC(ch1);   B(7) = fSSC(ch2);
    A(8) = fVAR(ch1);   B(8) = fVAR(ch2);
    A(9) = fWA(ch1);    B(9) = fWA(ch2);
    A(10)= fWL(ch1);    B(10)= fWL(ch2);
    A(11)= fZC(ch1);    B(11)= fZC(ch2);

    % 保存特征在TXT文件中
    fid = fopen('features/timef/zzzzzz.txt','a');
    for j=1:10
        fprintf(fid,'%g\t',A(j));
        fprintf(fid,'%g\t',B(j));
    end
    fprintf(fid,'%g\t',A(11));
    fprintf(fid,'%g\n',B(11));
    fclose(fid);
end

5分类

在这里插入图片描述

常用：LDA、KNN、SVM、ANN、MLP
康复系统、各个演示系统的分类模型均是ANN
网络结构：三层（输入层：通道数*特征数、隐含层、输出层：类别数）

完成的工作：
（1）偏瘫患者手部（张手、闭手、放松）三类识别
（2）正常人（12动作在线识别）6电极
（3）下肢6动作在线识别（6电极）
踝关节上拉、下勾；膝关节前伸、后拉；髋关节内收、外展
————————————————
版权声明：本文为优快云博主「九方先生」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.youkuaiyun.com/malele4th/article/details/79361526