sEMG的时域特征

梁小憨憨

已于 2023-09-01 15:27:13 修改

阅读量3.5k

点赞数 4

分类专栏：信号处理文章标签：信号处理 matlab

于 2021-07-24 15:36:05 首次发布

本文链接：https://blog.csdn.net/qq_41990294/article/details/119058722

版权

信号处理专栏收录该内容

30 篇文章 125 订阅

订阅专栏

本篇为《信号处理》系列博客的第十二篇，该系列博客主要记录信号处理相关知识的学习过程和自己的理解，方便以后查阅。

sEMG的时域特征

绝对均值（MAV）
均方根值（RMS）
方差（VAR）
过零点数（ZC）
willision赋值（WAMP）
特征提取总程序

绝对均值（MAV）

请添加图片描述

%% 绝对均值(MAV);反映sEMG信号的平均强度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的绝对均值的和
        MAV_fn = MAV_fn+sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
    end
    MAV(ch_i) = MAV_fn/FN;% 计算一个通道的绝对均值，并赋值
    MAV_fn = 0;% 累计变量清零，用于下一次累计计算
end

均方根值（RMS）

请添加图片描述

%% 均方根值(RMS);反映sEMG的有效值，表示个肌肉群在动作过程中的贡献力度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的均方根值的和
        RMS_fn = RMS_fn+sqrt(sum((channal_frame(ch_i, frame_i, :)).^2)/(FrameLen-1));
    end
    RMS(ch_i) = RMS_fn/FN;% 计算一个通道的均方根值，并赋值
    RMS_fn = 0;% 累计变量清零，用于下一次累计计算
end

方差（VAR）

请添加图片描述

%% 方差(VAR);反映信号的变化的剧烈程度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算当前通道上当前窗口的绝对均值
        meandate = sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
        % 计算一个通道上的所有窗的方差的和
        VAR_fn = VAR_fn+sum((channal_frame(ch_i, frame_i, :)-meandate).^2)/FrameLen;
    end
    VAR(ch_i) = VAR_fn/FN;% 计算一个通道的均方根值，并赋值
    VAR_fn = 0;% 累计变量清零，用于下一次累计计算
end

过零点数（ZC）

请添加图片描述

%% 过零点数(ZC)
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算当前通道上当前窗口的绝对均值
        meandate = sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
        % 计算一个通道上的所有窗的方差的和
        channal_frame_f = channal_frame(ch_i, frame_i, 1:FrameLen-1);
        channal_frame_b = channal_frame(ch_i, frame_i, 2:FrameLen);
        for signal_i = 1:FrameLen-1
            if abs(channal_frame_b(signal_i)-channal_frame_f(signal_i))>th_zc && (channal_frame_f(signal_i)-meandate)*(channal_frame_b(signal_i)-meandate)<0
                ZC_fn = ZC_fn+1;
            else
                ZC_fn = ZC_fn+0;
            end
        end
    end
    ZC(ch_i) = ZC_fn/FN;% 计算一个通道的均方根值，并赋值
    ZC_fn = 0;% 累计变量清零，用于下一次累计计算
end

willision赋值（WAMP）

请添加图片描述

%% willision幅值(WAMP);反映相应肌肉的收缩水平
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的方差的和
        channal_frame_f = channal_frame(ch_i, frame_i, 1:FrameLen-1);
        channal_frame_b = channal_frame(ch_i, frame_i, 2:FrameLen);
        for signal_i = 1:FrameLen-1
            if abs(channal_frame_f(signal_i)-channal_frame_b(signal_i))>th_wamp
                WAMP_fn = WAMP_fn+1;
            else
                WAMP_fn = WAMP_fn+0;
            end
        end
    end
    WAMP(ch_i) = WAMP_fn/FN;% 计算一个通道的均方根值，并赋值
    WAMP_fn = 0;% 累计变量清零，用于下一次累计计算
end

特征提取总程序

%%% 时域特征提取 %%%
clear;
clc;

%% 读取数据 %% 
% filename = '/home/al007/Matlab/聪姐论文程序/聪数据/放松/放松/放松数据/放松_Plot_and_Store_Rep_2.1.csv';
filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲总/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_15.csv';
P=csvread(filename,1,0);%从第1行第0列开始读取数据
[row, column]=size(P);%获取数据维度大小
sEMG = zeros(row,8);
for ch_i  = 2:2:column
    sEMG(:,ch_i/2)=P(:,ch_i);%获取一个通道的肌电数据
end
sEMG_L=length(sEMG(:,1));%获取第一个通道的数据量

%% 对信号进行IWD滤波 %%
for ch_i = 1:column/2
    [XD,CXD,LXD] = wden(sEMG(:, ch_i), 'heursure', 's', 'mln', 4, 'sym6');%对信号进行小波变换阈值去噪

    [C,L] = wavedec(XD, 5, 'sym6');%多级一维小波分解 改分解层数
    ca5 = appcoef(C,L,'sym6',5);%一维近似系数，计算一维信号的近似系数
    [cd1,cd2,cd3,cd4,cd5] = detcoef(C,L,[1,2,3,4,5]);%一维细节系数，获取小波分解的细节系数
    %利用数字滤波阈值 去除基线漂移 和 高频信号降噪
    %将第一层细节系数设置为0，进行高频信号降噪，将最后一层的近似系数设置为0，进行基线漂移的去除
    C5 = [zeros(size(ca5)); cd5; cd4; cd3; cd2; zeros(size(cd1))];
    sEMG(:, ch_i) = waverec(C5,L,'sym6');%多层次一维小波重构
end

%% 截取活动段 %% 
FrameLen = 64;% 设置帧的长度
FrameInc = 32;% 设置帧的步长

[action_start, energy_mean] = endpoint_detection(sEMG, sEMG_L);% 获取动作的起点

for ch_i = 1:column/2
    action_signal(:,ch_i) = sEMG(action_start*FrameInc:action_start*FrameInc+6000-1, ch_i);%截取3s长度的信号
end

action_signal_L = length(action_signal(:,ch_i));%获取第一个通道的数据量

%% 变量定义 %%
channal_sum = 8;% 肌电采集的通道总数

FrameLen = 256;% 的设置帧的长度，重叠窗处理使用
FrameInc = 32;% 的设置帧的步长，重叠窗处理使用

channal_frame = [];% 通道上，分帧后，帧的数据，存储每个通道上重叠窗的帧信息

MAV = [];% 存储每个通道上的绝对均值
MAV_fn = 0;% 存储一个通道上的绝对均值，用于累加

RMS = [];% 存储每个通道上的均方根值
RMS_fn = 0;% 存储一个通道上的均方根值，用于累加

VAR = [];% 存储每个通道上的方差
VAR_fn = 0;% 存储一个通道上的方差，用于累加

meandate = 0;% 一个通道上当前窗口的绝对均值

ZC = [];% 存储每个通道上的过零点数
ZC_fn = 0;% 存储一个通道上的过零点数，用于累加
th_zc = 0;% 计算过零点数用到的阈值，如何设置这个阈值？%%%%%%%%%

WAMP = [];% 存储每个通道上的willision幅值
WAMP_fn = 0;% 存储一个通道上的willision幅值，用于累加
th_wamp = 0;% 计算willision幅值用到的阈值，如何设置这个阈值？%%%%%%%%%

%% 信号进行重叠窗处理 %% 
for ch_i = 1:channal_sum
    channal_frame(ch_i,:,:) = enframe(action_signal(1:end-1,ch_i), FrameLen, FrameInc);% 进行分帧处理
end

FN = fix((action_signal_L-FrameLen+FrameInc)/FrameInc);% 计算总帧数，向0取进位

%% 绝对均值(MAV);反映sEMG信号的平均强度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的绝对均值的和
        MAV_fn = MAV_fn+sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
    end
    MAV(ch_i) = MAV_fn/FN;% 计算一个通道的绝对均值，并赋值
    MAV_fn = 0;% 累计变量清零，用于下一次累计计算
end

%% 均方根值(RMS);反映sEMG的有效值，表示个肌肉群在动作过程中的贡献力度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的均方根值的和
        RMS_fn = RMS_fn+sqrt(sum((channal_frame(ch_i, frame_i, :)).^2)/(FrameLen-1));
    end
    RMS(ch_i) = RMS_fn/FN;% 计算一个通道的均方根值，并赋值
    RMS_fn = 0;% 累计变量清零，用于下一次累计计算
end

%% 方差(VAR);反映信号的变化的剧烈程度
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算当前通道上当前窗口的绝对均值
        meandate = sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
        % 计算一个通道上的所有窗的方差的和
        VAR_fn = VAR_fn+sum((channal_frame(ch_i, frame_i, :)-meandate).^2)/FrameLen;
    end
    VAR(ch_i) = VAR_fn/FN;% 计算一个通道的均方根值，并赋值
    VAR_fn = 0;% 累计变量清零，用于下一次累计计算
end

%% 过零点数(ZC)
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算当前通道上当前窗口的绝对均值
        meandate = sum(abs(channal_frame(ch_i, frame_i, :)))/FrameLen;
        % 计算一个通道上的所有窗的方差的和
        channal_frame_f = channal_frame(ch_i, frame_i, 1:FrameLen-1);
        channal_frame_b = channal_frame(ch_i, frame_i, 2:FrameLen);
        for signal_i = 1:FrameLen-1
            if abs(channal_frame_b(signal_i)-channal_frame_f(signal_i))>th_zc && (channal_frame_f(signal_i)-meandate)*(channal_frame_b(signal_i)-meandate)<0
                ZC_fn = ZC_fn+1;
            else
                ZC_fn = ZC_fn+0;
            end
        end
    end
    ZC(ch_i) = ZC_fn/FN;% 计算一个通道的均方根值，并赋值
    ZC_fn = 0;% 累计变量清零，用于下一次累计计算
end

%% willision幅值(WAMP);反映相应肌肉的收缩水平
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 计算一个通道上的所有窗的方差的和
        channal_frame_f = channal_frame(ch_i, frame_i, 1:FrameLen-1);
        channal_frame_b = channal_frame(ch_i, frame_i, 2:FrameLen);
        for signal_i = 1:FrameLen-1
            if abs(channal_frame_f(signal_i)-channal_frame_b(signal_i))>th_wamp
                WAMP_fn = WAMP_fn+1;
            else
                WAMP_fn = WAMP_fn+0;
            end
        end
    end
    WAMP(ch_i) = WAMP_fn/FN;% 计算一个通道的均方根值，并赋值
    WAMP_fn = 0;% 累计变量清零，用于下一次累计计算
end

%%% 分帧能量法计算动作起点 %%%
function [action_start, energy_mean] = endpoint_detection(signal_8, signal_L)

% 幅值归一化到[-1,1]，为了加快数据收敛
signal_8 = double(signal_8);
for ch_i = 1:8
    signal_8(:,ch_i) = signal_8(:,ch_i) / max(abs(signal_8(:,ch_i)));
end


FrameLen = 64;% 设置帧的长度
FrameInc = 32;% 设置帧的步长

energy =7;% 这里是需要人为设置的阈值参数,当计算静息状态能量平均值时，设置成很大值

% 进行分帧处理
channal1 = enframe(signal_8(1:end-1,1), FrameLen, FrameInc);
channal2 = enframe(signal_8(1:end-1,2), FrameLen, FrameInc);
channal3 = enframe(signal_8(1:end-1,3), FrameLen, FrameInc);
channal4 = enframe(signal_8(1:end-1,4), FrameLen, FrameInc);
channal5 = enframe(signal_8(1:end-1,5), FrameLen, FrameInc);
channal6 = enframe(signal_8(1:end-1,6), FrameLen, FrameInc);
channal7 = enframe(signal_8(1:end-1,7), FrameLen, FrameInc);
channal8 = enframe(signal_8(1:end-1,8), FrameLen, FrameInc);

FN = fix((signal_L-FrameLen+FrameInc)/FrameInc);% 计算总帧数，向0取进位

channal_energy = zeros(1,8);% 存储每个通道的能量
energy_frame = 0;% 存储一帧的总能量
frame_num = 0;% 帧的计数标志，分帧能量法中的计数帧
energy_frame_sum = 0;% 整段信号的能量总值，在计算静息状态能量平均值用的
energy_frame_num = 0;% 能量不为0的帧的总数，在计算静息状态能量平均值用的
energy_mean = 0;% 能量的平均值，在计算静息状态能量平均值时使用

for FN_i = 10:FN-1
    
    % 求当前帧下每个通道的能量
    channal_energy(1) = sum(channal1(FN_i,:).^2);
    channal_energy(2) = sum(channal2(FN_i,:).^2);
    channal_energy(3) = sum(channal3(FN_i,:).^2);
    channal_energy(4) = sum(channal4(FN_i,:).^2);
    channal_energy(5) = sum(channal5(FN_i,:).^2);
    channal_energy(6) = sum(channal6(FN_i,:).^2);
    channal_energy(7) = sum(channal7(FN_i,:).^2);
    channal_energy(8) = sum(channal8(FN_i,:).^2);
    
    energy_frame = sum(channal_energy);% 计算当前帧下8个通道的能量总和

%     fprintf('energy_frame：\t%f\n', energy_frame);
    
    % 计算静息状态下能量的平均值
    if energy_frame > 0
        energy_frame_sum = energy_frame_sum+energy_frame;
        energy_frame_num = energy_frame_num+1;
    end
    
    if energy_frame > energy        % 当该帧的能量大于设定的阈值时，计数加一
        frame_num = frame_num+1;
    else                            % 当存在一个帧的嫩能量小于，设定的阈值时，清零标志位
        frame_num = 0;
    end
    
    if frame_num == 3               % 如果计数标志位等于3，则检测到标志位，将当前的帧数返回
        action_start = FN_i-2;
        energy_mean = energy_frame_sum/energy_frame_num;% 计算静息状态下能量的平均值
        break;
    end
    
    action_start = 1;
end

energy_mean = energy_frame_sum/energy_frame_num;% 计算静息状态下能量的平均值