表面肌电sEMG特征提取的Matlab程序

提取的特征值包括：

时域——RMS，MAV，ZC，VAR

频域——平均功率频率MPF，中值频率MF

数据为：thisdata，带入你的数据

 

% 时域特征
% 绝对平均值(Mean Absolute Value,MAV)
%思路:截取一段时间（定义滑动窗口）length_t,滑动大小delta_t,单位ms
length_t=1000;
delta_t=50;
j=1;
for i=1:delta_t:size(thisdata,1)
    if i+length_t>size(thisdata,1)
        break;
    end
    meandata(j)=(sum(abs(thisdata(i:i+length_t))))/length_t;
    j=j+1;
end
figure(2)
plot(meandata);
xlabel('取样点数/1');
ylabel('肌电电压值/uV');

% 均方根值(Root Mean Square,RMS)
%思路:截取一段时间（定义滑动窗口）length_t,滑动大小delta_t,单位ms
length_t=1000;
delta_t=50;
j=1;
for i=1:delta_t:size(thisdata,1)
    if i+length_t>size(thisdata,1)
        break;
    end
    rms_data(j)=sqrt((sum((thisdata(i:i+length_t)).^2))/length_t);
    j=j+1;
end
figure(2)
plot(rms_data);
xlabel('取样点数/1');
ylabel('肌电电压值/uV');

% 方差(Variance,VAR)

%思路:截取一段时间（定义滑动窗口）length_t,滑动大小delta_t,单位ms
length_t=1000;
delta_t=50;
j=1;
for i=1:delta_t:size(thisdata,1)
    if i+length_t>size(thisdata,1)
        break;
    end
    mean_thisdata=mean(thisdata(i:i+length_t));
    rms_data(j)=(sum((thisdata(i:i+length_t)-mean_thisdata).^2))/length_t;
    j=j+1;
end
figure(2)
plot(rms_data);
xlabel('取样点数/1');
ylabel('肌电电压值/uV^2');

% 过零点数(Zero Crossing,ZC)
%思路:截取一段时间（定义滑动窗口）length_t,滑动大小delta_t,单位ms
%阈值：相邻值大于一定值，threshold
length_t=300;
delta_t=50;
ZC_threshold=50;
j=0;
for i=1:delta_t:size(thisdata,1)
    if i+length_t>size(thisdata,1)
        break;
    end
    ZC_data(j+1)=0;
    meanvalue_thisdata=mean(thisdata(i:i+length_t-1));
    for k=1:length_t-1        
        if (thisdata(i+k-1)-meanvalue_thisdata)*(thisdata(i+k)-meanvalue_thisdata)<0
               if abs(thisdata(i+k-1)-thisdata(i+k))>ZC_threshold
                    ZC_data(j+1)=ZC_data(j+1)+1;
               end           
        end
    end
    j=j+1;
end
figure(2)
plot(ZC_data);
xlabel('取样点数/1');
ylabel('过零点数/1');



% Willison幅值(Willison Amplitude,WA)
信号幅值的变化次数

% 频域特征
%% 平均功率频率(Mean Power Frequency,MPF)
%思路:用窗函数截取数据，然后分析
%阈值：相邻值大于一定值，threshold
delta_t=50;
fs=1000;%采样频率
length_t=1024;%窗函数的数据长度,最好为2的次方，否则不足的数据会被补0
%L1=length(thisdata);%数据总长度，和数据点数一样
%w=hamming(length_t);%窗长度为   的汉明窗,解决栅栏效应
w=hamming(length_t);

j=1;
k=1;
for i=1:floor(length(thisdata)/delta_t)
    if delta_t*(i-1)+1+(length_t-1)>length(thisdata)
        break;
    end
    getdata=w.*thisdata(delta_t*(i-1)+1:delta_t*(i-1)+1+(length_t-1))';
    mean_getdata=mean(getdata);%求均值
    cx1=xcorr(getdata-mean_getdata,'unbiased');%求自相关函数，无偏估计

    cxk1=fft(cx1,length_t);%傅里叶变换
    px1=abs(cxk1);%求功率谱密度
    pxx1=10*log10(px1);

    f1=(0:length_t-1)*fs/length_t;
%     plot(f1(1:length_t/2),pxx1(1:length_t/2))
%     xlabel('频率/Hz');ylabel('功率谱/dB');
%     %title('平均功率谱图');
%     grid on  %做功率谱图
    
    df1=fs/length_t;
    p1=(sum(px1(1:length_t/2-1))+sum(px1(1:length_t/2)))/2.*df1;
    pf1=(sum(px1(1:length_t/2-1).*[1:length_t/2-1]'.*df1)+sum(px1(1:length_t/2).*[1:length_t/2]'.*df1))/2*df1;
    MPF1(j)=pf1/p1; %求平均功率频率
    j=j+1;
    
    N1=1;pp1=0;
    while abs(pp1-p1/2)>(px1(N1)+px1(N1+1))/2*df1
         pp1=pp1+(px1(N1)+px1(N1+1))/2*df1;
         N1=N1+1;
    end
    n_1=(N1+N1+1)/2;
    MF1(k)=df1*n_1; %求中值频率
    k=k+1;
end
figure(2)
plot(MPF1);
xlabel('取样点数/1');
ylabel('平均功率频率/Hz');
figure(3)
plot(MF1);
xlabel('取样点数/1');
ylabel('中值频率/Hz');