sEMG的频域特征

已于 2023-04-16 09:01:50 修改
阅读量3.2k 收藏 56
点赞数 5
分类专栏： 信号处理 文章标签： 信号处理 matlab
于 2021-07-24 15:54:04 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_41990294/article/details/119059249
版权

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

sEMG的频域特征

功率谱最大值（MPS）

请添加图片描述

%% 功率谱最大值(MPS)
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 根据自相关函数计算功率谱
%         Cross_correlation = xcorr(squeeze(channal_frame(ch_i, frame_i, :)),'unbiased'); %计算序列的自相关函数
%         Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
%         Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
        % 根据周期图法计算功率谱
%         window = boxcar(length(channal_frame(ch_i, frame_i, :))); %矩形窗
%         [Psd_fn,f] = periodogram(squeeze(channal_frame(ch_i, frame_i, :)), window); %直接法,周期图功率谱密度估计
        % 现代谱估计法
        Psd_fn = pyulear(reshape(squeeze(channal_frame(ch_i, frame_i, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法
        
        Psd_fn = Psd_fn / (max(Psd_fn)-min(Psd_fn));% 归一化到0~1
        Psd_fn = 10*log10(Psd_fn(1:FrameLen/2)+0.000001);% 换算成分贝，方便观察
        MPS_fn = max(Psd_fn);% 获取当前窗口的最大功率值
        % 获取所有窗口中的最大功率,也即当前通道的最大功率
        if MPS_fn > MPS_ch
            MPS_ch = MPS_fn;
        end
    end
    MPS(ch_i) = MPS_ch;% 计算一个通道的绝对均值，并赋值
    Psd_fn = 0;% 清零，用于下一次循环
    MPS_ch = 0;% 清零，用于下一次循环
end

%%% 三种功率谱的计算结果展示
figure(1);%Create figure window
subplot(2,1,1);
plot(reshape(squeeze(channal_frame(1, 3, :)),1,FrameLen));
title('原始信号');
grid on

subplot(2,1,2);
% 根据自相关函数计算功率谱
% Cross_correlation = xcorr(squeeze(channal_frame(1, 3, :)),'unbiased'); %计算序列的自相关函数
% Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
% Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
% 根据周期图法计算功率谱
% window = boxcar(length(channal_frame(1, 3, :))); %矩形窗
% [Psd_fn,f] = periodogram(squeeze(channal_frame(1, 3, :)), window); %直接法,周期图功率谱密度估计
% 现代谱估计法
Psd_fn = pyulear(reshape(squeeze(channal_frame(1, 3, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法

% Psd_fn = Psd_fn / max(Psd_fn);% 归一化到0~1
Psd_fn = 10*log10(Psd_fn(1:FrameLen/2)+1e-15);% 换算成分贝，方便观察

index=0:round(FrameLen/2-1);
k=index*SampleFre/FrameLen;

plot(k, Psd_fn);
title('AR谱估计');
grid on;

中值频率（MF）

请添加图片描述

%% 中值频率
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        MF_fn = MF_fn+meanfreq(reshape(squeeze(channal_frame(ch_i, frame_i, :)),FrameLen,1),SampleFre);
    end
    MF(ch_i) = MF_fn/FN;% 计算一个通道的中值频率，并赋值
    MF_fn = 0;% 累计变量清零，用于下一次累计计算
end

平均功率频率（MPF）

请添加图片描述

%% 平均功率频率(MPF)
index=0:round(FrameLen/2-1);
f=index*SampleFre/FrameLen;% 计算频率值

for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 根据自相关函数计算功率谱
%         Cross_correlation = xcorr(squeeze(channal_frame(ch_i, frame_i, :)),'unbiased'); %计算序列的自相关函数
%         Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
%         Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
        % 根据周期图法计算功率谱
%         window = boxcar(length(channal_frame(ch_i, frame_i, :))); %矩形窗
%         [Psd_fn, f] = periodogram(squeeze(channal_frame(ch_i, frame_i, :)), window); %直接法,周期图功率谱密度估计
        % 现代谱估计法
        Psd_fn = pyulear(reshape(squeeze(channal_frame(ch_i, frame_i, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法
             
        MPF_fn = MPF_fn+sum(f.*reshape(Psd_fn(1:FrameLen/2), 1, FrameLen/2))/sum(Psd_fn(1:FrameLen/2));
    end
    MPF(ch_i) = MPF_fn/FN;% 计算一个通道的绝对均值，并赋值
    MPF_fn = 0;% 变量清零，用于下一次累计计算
    Psd_fn = 0;% 变量清零，用于下一次累计计算
end

特征提取总程序

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

%% 读取数据 %% 
% filename = '/home/al007/Matlab/聪姐论文程序/聪数据/放松/放松/放松数据/放松_Plot_and_Store_Rep_2.1.csv';
filename = '/home/al007/Matlab/聪姐论文程序/聪数据/肩屈曲数据/肩屈曲_Plot_and_Store_Rep_5.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 = [];% 通道上，分帧后，帧的数据，存储每个通道上重叠窗的帧信息

SampleFre = 2000;%采样频率

Cross_correlation = 0;% 自相关函数
Psd_fn = 0;% 当前帧的功率谱密度函数
MPS_fn = 0;% 当前帧的最大功率
MPS_ch = 0;% 当前通道的最大功率
MPS = [];% 存储每个通道上的频率普最大值

MF = [];% 存储每个通道上的中值频率
MF_fn = 0;% % 存储一个通道上的中值频率，用于累加

MPF = [];% 存储每个通道上的平均功率频率
MPF_fn = 0;% % 存储一个通道上的单边平均功率频率，用于积分运算

%% 信号进行重叠窗处理 %% 
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取进位

%% 功率谱最大值(MPS)
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 根据自相关函数计算功率谱
%         Cross_correlation = xcorr(squeeze(channal_frame(ch_i, frame_i, :)),'unbiased'); %计算序列的自相关函数
%         Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
%         Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
        % 根据周期图法计算功率谱
%         window = boxcar(length(channal_frame(ch_i, frame_i, :))); %矩形窗
%         [Psd_fn,f] = periodogram(squeeze(channal_frame(ch_i, frame_i, :)), window); %直接法,周期图功率谱密度估计
        % 现代谱估计法
        Psd_fn = pyulear(reshape(squeeze(channal_frame(ch_i, frame_i, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法
        
        Psd_fn = Psd_fn / max(Psd_fn);% 归一化到0~1
        Psd_fn = 10*log10(Psd_fn(1:FrameLen/2)+0.000001);% 换算成分贝，方便观察
        MPS_fn = max(Psd_fn);% 获取当前窗口的最大功率值
        % 获取所有窗口中的最大功率,也即当前通道的最大功率
        if MPS_fn > MPS_ch
            MPS_ch = MPS_fn;
        end
    end
    MPS(ch_i) = MPS_ch;% 计算一个通道的绝对均值，并赋值
    Psd_fn = 0;% 清零，用于下一次循环
    MPS_ch = 0;% 清零，用于下一次循环
end

%%% 三种功率谱的计算结果展示
figure(1);%Create figure window
subplot(2,1,1);
plot(reshape(squeeze(channal_frame(1, 3, :)),1,FrameLen));
title('原始信号');
grid on

subplot(2,1,2);
% 根据自相关函数计算功率谱
% Cross_correlation = xcorr(squeeze(channal_frame(1, 3, :)),'unbiased'); %计算序列的自相关函数
% Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
% Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
% 根据周期图法计算功率谱
% window = boxcar(length(channal_frame(1, 3, :))); %矩形窗
% [Psd_fn,f] = periodogram(squeeze(channal_frame(1, 3, :)), window); %直接法,周期图功率谱密度估计
% 现代谱估计法
Psd_fn = pyulear(reshape(squeeze(channal_frame(1, 3, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法

% Psd_fn = Psd_fn / max(Psd_fn);% 归一化到0~1
Psd_fn = 10*log10(Psd_fn(1:FrameLen/2)+0.000001);% 换算成分贝，方便观察

index=0:round(FrameLen/2-1);
k=index*SampleFre/FrameLen;

plot(k, Psd_fn);
title('AR谱估计');
grid on;

%% 中值频率
for ch_i = 1:channal_sum
    for frame_i = 1:FN
        MF_fn = MF_fn+meanfreq(reshape(squeeze(channal_frame(ch_i, frame_i, :)),FrameLen,1),SampleFre);
    end
    MF(ch_i) = MF_fn/FN;% 计算一个通道的中值频率，并赋值
    MF_fn = 0;% 累计变量清零，用于下一次累计计算
end

%% 平均功率频率(MPF)
index=0:round(FrameLen/2-1);
f=index*SampleFre/FrameLen;% 计算频率值

for ch_i = 1:channal_sum
    for frame_i = 1:FN
        % 根据自相关函数计算功率谱
%         Cross_correlation = xcorr(squeeze(channal_frame(ch_i, frame_i, :)),'unbiased'); %计算序列的自相关函数
%         Cross_correlation_fft = fft(Cross_correlation);% 对自相关函数做傅里叶变换
%         Psd_fn = 2 * abs(Cross_correlation_fft(1:fix(length(Cross_correlation_fft)/2))) / length(Cross_correlation_fft);% 得到功率频谱
        % 根据周期图法计算功率谱
%         window = boxcar(length(channal_frame(ch_i, frame_i, :))); %矩形窗
%         [Psd_fn, f] = periodogram(squeeze(channal_frame(ch_i, frame_i, :)), window); %直接法,周期图功率谱密度估计
        % 现代谱估计法
        Psd_fn = pyulear(reshape(squeeze(channal_frame(ch_i, frame_i, :)),1,FrameLen),FrameLen,FrameLen);% 自回归功率谱密度估计-Yule-Walker方法
             
        MPF_fn = MPF_fn+sum(f.*reshape(Psd_fn(1:FrameLen/2), 1, FrameLen/2))/sum(Psd_fn(1:FrameLen/2));
    end
    MPF(ch_i) = MPF_fn/FN;% 计算一个通道的绝对均值，并赋值
    MPF_fn = 0;% 变量清零，用于下一次累计计算
    Psd_fn = 0;% 变量清零，用于下一次累计计算
end
梁小憨憨
关注
  • 5
    点赞
  • 踩
  • 56
    收藏
    觉得还不错? 一键收藏
  • 5
    评论
专栏目录
信号处理频域特征
weixin_39107270的博客
05-18 2505
倒频谱能较好地检测出功率谱上的周期成分，通常在功率谱上无法对边频的总体水平作出定量估计，而倒频谱对边频成分具有“概括”能力，能较明显地显示出功率谱上的周期成分，将原来谱上成族的边频带谱线简化为单根谱线，便于观察，而齿轮发生故障时的振动频谱具有的边频带一般都具有等间隔（故障频率）的结构，利用倒频谱这个优点，可以检测出功率谱中难以辨识的周期性信号。2.通常我们关心的都是正频率区间的结果，有两种截取方法，一种是在fftshift的结果中截后半段，一种是在fft的结果中截前半段，其结果是一样的。
时域、频域特征提取
06-28
提取数据的时域、频域指标，得到相应的时域和频域特征
表面肌电sEMG特征提取的Matlab程序
04-30
提取的特征值包括：时域——RMS，MAV，ZC，VAR频域——平均功率频率MPF，中值频率MF
sEMG项目总结（5）sEMG分析
malele4th
02-24 2万+
sEMG分析 目录 sEMG分析 目录 1肌电信号 2肌电信号采集 3预处理滤波去噪 4特征提取与选择 5分类 1肌电信号 肌电信号的产生 表面肌电信号是由多个运动单元发放的动作电位序列，在皮肤表面呈现的时间和空间上综合叠加的结果 肌电信号特点 （1）sEMG是一维时间动作电位序列 （2）sEMG是交流信号，幅值一般和肌肉运 动力度成正比 （...
频域特征指标及其MATLAB代码实现（重心频率、均方频率、均方根频率、频率方差、频率标准差）
fengzhuqiaoqiu的博客
10-31 2万+
以前我们梳理过常用的有量纲时域特征指标和无量纲时域特征指标，其实特征提取方法中，一些频谱特征指标也是十分有用的。 下述频域特征指标是根据信号的功率谱求得的。在故障诊断或者其他应用场景下，可以与时域指标以及其他指标一同作为特征向量的一部分，丰富特征量种类，提升诊断正确率。 1.理论介绍 1.0.功率谱 功率谱的相关介绍在之前讲到过： Mr.看海：信号频域分析方法的理解（频谱、能量谱、功率谱、倒频谱、小波分析） 下述频域特征都是在功率谱基础上计算求得的，但是为什么不用频谱呢？个人认为是因为功率谱
频域统计特征介绍
zailushag的博客
01-14 1万+
频域特征介绍 频域是通过傅里叶变换将时域进行变换，傅里叶变换具体请看https://blog.csdn.net/l494926429/article/details/51818012/，这是一位大牛，解释的很清除并且易懂。
频域特征提取
热门推荐
qq_37240982的博客
07-18 3万+
频域特征 （1）提取的频域特征 频域分析可按频率观察信号特征，一般情况下，时域的分析更加直观，而频域的表示更加简洁，在频域上观察信号使得问题的分析更加深刻和便捷。目前来说，从时域到频域已成为信号分析的趋势。但是，这两种分析手段是相互联系、相辅相成和互有优势的。常规的频谱分析是指对信号进行傅立叶变换以进行分析。频谱分析包括幅度频谱和相位频谱且幅度频谱是最常用的。当减速器关键部件的健康状态发生改变时，样本数据信号频谱中的频率分量会相应地改变。故可以先通过分析振动信号的频域特征准确地表征信号频谱信息，然后获悉在不
sEMG 信号 预处理：去噪，分割，特征提取_sEMG去噪_
10-04
sEMG 信号 预处理：去噪，分割，特征提取
sEMG肌电信号.zip
09-03
四通道的表面肌电信号，共5组数据，分别为拇指，食指，中指，无名指，小指的运动状态
MATLAB频域特征提取.zip
12-03
MATLAB频域特征提取.
时域频域29个特征提取.zip
09-08
提取信号时域，频域特征共29维。可以用来进行信号分析
声音信号的时频域特征分析
11-22
针对三种声音分析了时频域特征，总结各特征的表达能力。
对肌电信号计算积分肌电值+均方根值+中值频率+平均功率频率
05-02
matlab程序，可对肌电信号进行处理。
sEMG基于肌电信号的深度学习数据集
04-04
delsys设备采集的表面肌电信号，16个手势动作，每个动作维持6s，休息4s，进行6次循环。类别标签通过最大面积法修正，效果较好。
SEMG检测电极的设计
01-20
SEMG（Surface electromyogram）即表面肌电信号，是指当骨骼肌收缩时，肌纤维产生的微弱动作电位在皮肤表面募集，并通过表面电极所采集到的信号。SEMG检测电极是用来检测表面肌电信号的有源电极，它是一种无创检测...
一维信号的频域特征分析python
abc1234abcdefg的博客
03-14 1万+
1.讲解了功率谱的实现 2.利用功率谱为参考的频域特征求取
时域、频域、时频特征提取技术
u011661076的博客
04-21 2万+
文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。 一、时域特征提取技术 当我们想要研究一个生物信号时，有时从原始信号中获取的信息并不明显。因此，需要提取一些特征来表示信号。可以研究代表信号波形的特征，也可以把信号的样本看作是来自随机源，因此可以得到一个分布，从中提取信息。通过分析它们在这个分布中的统计特性，或者在某些情况下的混沌行为，可以从混沌分形理论中进行测量。下面列出了一些提取原始生物电信号(时域)特征的最流行的方法。 1、从波形中提取特征：提取的最直接的特征是可以直接从信号中
sEMG动作分类前馈环节设计
最新发布
06-09
sEMG动作分类前馈环节的设计主要涉及以下几个方面： 1. 信号预处理：sEMG信号本身存在着干扰，如电源干扰、肌肉震颤等，需要进行滤波和去噪处理，以提高信号质量。 2. 特征提取：特征提取是将原始sEMG信号转化为能够描述运动特征的特征向量的过程。常用的特征包括时域特征、频域特征、时频域特征等。 3. 特征选择：特征选择是从提取的特征中选取最具代表性的特征。常用的方法包括相关系数、卡方检验、互信息等。 4. 分类器设计：分类器是将输入的特征向量映射到对应的动作标签的模型。常用的分类器包括支持向量机、决策树、神经网络等。 在设计前馈环节时，需要根据具体的应用场景和数据特点，选择合适的预处理方法、特征提取方法、特征选择方法和分类器。同时需要进行充分的实验验证，不断优化和调整设计方案，以提高分类性能和实时性。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值