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

sEMG的频域特征

• 功率谱最大值（MPS）
• 中值频率（MF）
• 平均功率频率（MPF）
• 特征提取总程序

功率谱最大值（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);% 归一化到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;

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