sEMG项目总结（5）sEMG分析

最新推荐文章于 2024-05-19 22:45:35 发布

九方先生

最新推荐文章于 2024-05-19 22:45:35 发布

阅读量2.2w

点赞数 67

分类专栏： sEMG项目

本文链接：https://blog.csdn.net/malele4th/article/details/79361526

版权

sEMG项目专栏收录该内容

7 篇文章 127 订阅

订阅专栏

sEMG分析

（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电极）
踝关节上拉、下勾；膝关节前伸、后拉；髋关节内收、外展

九方先生

关注

67
点赞
踩
344

收藏

觉得还不错? 一键收藏
20
评论
sEMG项目总结（5）sEMG分析

sEMG分析目录sEMG分析目录1肌电信号2肌电信号采集3预处理滤波去噪4特征提取与选择5分类1肌电信号肌电信号的产生表面肌电信号是由多个运动单元发放的动作电位序列，在皮肤表面呈现的时间和空间上综合叠加的结果肌电信号特点（1）sEMG是一维时间动作电位序列（2）sEMG是交流信号，幅值一般和肌肉运动力度成正比（...
复制链接

扫一扫