肌电信号积分肌电值+均方根值+中值频率+平均功率频率.zip

s= load('09121052HF.CSV'); fs=3000; x1=s(1:16384,1); x2=s(1:16384,2); x3=s(1:16384,3); x4=s(1:16384,4); x5=s(1:16384,5); x6=s(1:16384,6); iemg1=sum(abs(x1))/length(x1) iemg2=sum(abs(x2))/length(x2) iemg3=sum(abs(x3))/length(x3) iemg4=sum(abs(x4))/length(x4) iemg5=sum(abs(x5))/length(x5) iemg6=sum(abs(x6))/length(x6) %求积分肌电值 rms1=sqrt(sum(x1.^2)/length(x1)) rms2=sqrt(sum(x2.^2)/length(x2)) rms3=sqrt(sum(x3.^2)/length(x3)) rms4=sqrt(sum(x4.^2)/length(x4)) rms5=sqrt(sum(x5.^2)/length(x5)) rms6=sqrt(sum(x6.^2)/length(x6)) %求均方根值 L1=length(x1); cx1=xcorr(x1,'unbiased'); cxk1=fft(cx1,L1); px1=abs(cxk1);%求功率谱密度 pxx1=10*log10(px1); f1=(0:L1-1)*fs/L1; subplot(3,2,1),plot(f1(1:L1/2),pxx1(1:L1/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df1=fs/L1; p1=(sum(px1(1:L1/2-1))+sum(px1(1:L1/2)))/2.*df1; pf1=(sum(px1(1:L1/2-1).*[1:L1/2-1]'.*df1)+sum(px1(1:L1/2).*[1:L1/2]'.*df1))/2*df1; MPF1=pf1/p1 %求平均功率频率 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=df1*n_1 %求中值频率 L2=length(x2); cx2=xcorr(x2,'unbiased'); cxk2=fft(cx2,L2); px2=abs(cxk2);%求功率谱密度 pxx2=10*log10(px2); f2=(0:L2-1)*fs/L2; subplot(3,2,2),plot(f2(1:L2/2),pxx2(1:L2/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df2=fs/L2; p2=(sum(px2(1:L2/2-1))+sum(px2(1:L2/2)))/2.*df2; pf2=(sum(px2(1:L2/2-1).*[1:L2/2-1]'.*df2)+sum(px2(1:L2/2).*[1:L2/2]'.*df2))/2*df2; MPF2=pf2/p2 %求平均功率频率 N2=1;pp2=0; while abs(pp2-p2/2)>(px2(N2)+px2(N2+1))/2*df2 pp2=pp2+(px2(N2)+px2(N2+1))/2*df2; N2=N2+1; end n_2=(N2+N2+1)/2; MF2=df2*n_2 %求中值频率 L3=length(x3); cx3=xcorr(x3,'unbiased'); cxk3=fft(cx3,L3); px3=abs(cxk3);%求功率谱密度 pxx3=10*log10(px3); f3=(0:L3-1)*fs/L3; subplot(3,2,3),plot(f3(1:L3/2),pxx3(1:L3/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df3=fs/L3; p3=(sum(px3(1:L3/2-1))+sum(px3(1:L3/2)))/2.*df3; pf3=(sum(px3(1:L3/2-1).*[1:L3/2-1]'.*df3)+sum(px3(1:L3/2).*[1:L3/2]'.*df3))/2*df3; MPF3=pf3/p3 %求平均功率频率 N3=1;pp3=0; while abs(pp3-p3/2)>(px3(N3)+px3(N3+1))/2*df3 pp3=pp3+(px3(N3)+px3(N3+1))/2*df3; N3=N3+1; end n_3=(N3+N3+1)/2; MF3=df3*n_3 %求中值频率 L4=length(x4); cx4=xcorr(x4,'unbiased'); cxk4=fft(cx4,L4); px4=abs(cxk4);%求功率谱密度 pxx4=10*log10(px4); f4=(0:L4-1)*fs/L4; subplot(3,2,4),plot(f4(1:L4/2),pxx4(1:L4/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df4=fs/L4; p4=(sum(px4(1:L4/2-1))+sum(px4(1:L4/2)))/2.*df4; pf4=(sum(px4(1:L4/2-1).*[1:L4/2-1]'.*df4)+sum(px4(1:L4/2).*[1:L4/2]'.*df4))/2*df4; MPF4=pf4/p4 %求平均功率频率 N4=1;pp4=0; while abs(pp4-p4/2)>(px4(N4)+px4(N4+1))/2*df4 pp4=pp4+(px4(N4)+px4(N4+1))/2*df4; N4=N4+1; end n_4=(N4+N4+1)/2; MF4=df4*n_4 %求中值频率 L5=length(x5); cx5=xcorr(x5,'unbiased'); cxk5=fft(cx5,L5); px5=abs(cxk5);%求功率谱密度 pxx5=10*log10(px5); f5=(0:L5-1)*fs/L5; subplot(3,2,5),plot(f5(1:L5/2),pxx5(1:L5/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df5=fs/L5; p5=(sum(px5(1:L5/2-1))+sum(px5(1:L5/2)))/2.*df5; pf5=(sum(px5(1:L5/2-1).*[1:L5/2-1]'.*df5)+sum(px5(1:L5/2).*[1:L5/2]'.*df5))/2*df5; MPF5=pf5/p5 %求平均功率频率 N5=1;pp5=0; while abs(pp5-p5/2)>(px5(N5)+px5(N5+1))/2*df5 pp5=pp5+(px5(N5)+px5(N5+1))/2*df5; N5=N5+1; end n_5=(N5+N5+1)/2; MF5=df5*n_5 %求中值频率 L6=length(x6); cx6=xcorr(x6,'unbiased'); cxk6=fft(cx6,L6); px6=abs(cxk6);%求功率谱密度 pxx6=10*log10(px6); f6=(0:L6-1)*fs/L6; subplot(3,2,6),plot(f6(1:L6/2),pxx6(1:L6/2)) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 df6=fs/L6; p6=(sum(px6(1:L6/2-1))+sum(px6(1:L6/2)))/2.*df6; pf6=(sum(px6(1:L6/2-1).*[1:L6/2-1]'.*df6)+sum(px6(1:L6/2).*[1:L6/2]'.*df6))/2*df6; MPF6=pf6/p6 %求平均功率频率 N6=1;pp6=0; while abs(pp6-p6/2)>(px6(N6)+px6(N6+1))/2*df6 pp6=pp6+(px6(N6)+px6(N6+1))/2*df6; N6=N6+1; end n_6=(N6+N6+1)/2; MF6=df6*n_6 %求中值频率 IEMG=[iemg1,iemg2,iemg3,iemg4,iemg5,iemg6] RMS=[rms1,rms2,rms3,rms4,rms5,rms6] MPF=[MPF1,MPF2,MPF3,MPF4,MPF5,MPF6] MF=[MF1,MF2,MF3,MF4,MF5,MF6]