devbox
Gausst松鼠会
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

sEMG的频域特征_semg信号常用的频谱特征提取方法

作者：Gausst松鼠会 | 2024-03-27 20:27:09

赞

踩

semg信号常用的频谱特征提取方法

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

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;

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53

中值频率（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

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8

平均功率频率（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

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22

特征提取总程序

%%% 频域特征提取 %%%
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

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
  • 76
  • 77
  • 78
  • 79
  • 80
  • 81
  • 82
  • 83
  • 84
  • 85
  • 86
  • 87
  • 88
  • 89
  • 90
  • 91
  • 92
  • 93
  • 94
  • 95
  • 96
  • 97
  • 98
  • 99
  • 100
  • 101
  • 102
  • 103
  • 104
  • 105
  • 106
  • 107
  • 108
  • 109
  • 110
  • 111
  • 112
  • 113
  • 114
  • 115
  • 116
  • 117
  • 118
  • 119
  • 120
  • 121
  • 122
  • 123
  • 124
  • 125
  • 126
  • 127
  • 128
  • 129
  • 130
  • 131
  • 132
  • 133
  • 134
  • 135
  • 136
  • 137
  • 138
  • 139
  • 140
  • 141
  • 142
  • 143
  • 144
  • 145
  • 146
  • 147
  • 148
  • 149
  • 150
  • 151
  • 152
  • 153
  • 154
声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/Gausst松鼠会/article/detail/326729
推荐阅读