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]