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

2019-04-13 22:01发布

站内文章 / 模拟电子
11154 0 1015
提取的特征值包括： 时域——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');  

热门文章