newData = importdata('Data.txt', '\t', 2);
data=newData.data;
A=9.84*[data(:,2) data(:,3) data(:,4)];
AA=A;
W=2*pi*[data(:,5) data(:,6) data(:,7)];
Angle=[data(:,8) data(:,9) data(:,10)]/(180*pi);
F=1;
Q=[1 0 0 0];
deltaT=0.01;
G=[0 0 9.84];
Ts=[1 1 1];
Te=[1 1 1];
winSize=20;
for i=1:length(A)
% %解算四元数
% preQ=Q;
% Q(1)=Q(1)+0.5*deltaT*(-preQ(2)*W(i,1)-preQ(3)*W(i,2)-preQ(4)*W(i,3));
% Q(2)=Q(2)+0.5*deltaT*(preQ(1)*W(i,1)+preQ( 3)*W(i,3)-preQ(4)*W(i,2));
% Q(3)=Q(3)+0.5*deltaT*(preQ(1)*W(i,2)-preQ(2)*W(i,3)+preQ(4)*W(i,1));
% Q(4)=Q(4)+0.5*deltaT*(preQ(1)*W(i,3)+preQ(2)*W(i,2)-preQ(3)*W(i,1));
% %根据姿态角求得旋转矩阵,并却出重力分量
% Cbe=[Q(1)*Q(1)+Q(2)*Q(2)-Q(3)*Q(3)-Q(4)*Q(4) 2*(Q(2)*Q(3)-Q(1)*Q(4)) 2*(Q(1)*Q(3)+Q(2)*Q(4));
% 2*(Q(1)*Q(4)) Q(1)*Q(1)-Q(2)*Q(2)+Q(3)*Q(3)-Q(4)*Q(4) 2*(Q(3)*Q(4)-Q(1)*Q(2));
% 2*(Q(2)*Q(4)-Q(1)*Q(3)) 2*(Q(1)*Q(2)+Q(3)*Q(4)) Q(1)*Q(1)-Q(2)*Q(2)-Q(3)*Q(3)+Q(4)*Q(4)];
%计算旋转矩阵
CbeZYX=eul2rotm([Angle(i,3),Angle(i,2),Angle(i,1)]);
AA(i,:)=AA(i,:)-G;
%去除重力加速度分量
A(i,:)=(CbeZYX*A(i,:).').'-G;
end
% figure(F)
% F=F+1;
% subplot(2,2,1);
% plot(1:length(A),A(:,1));
% legend('A');
% hold on;
% plot(1:length(AA),AA(:,1));
% legend('AA');
% grid on;
for j=1:3
zero=mean(A(1:20,j));
A(:,j)=A(:,j)-zero;
end
figure(F)
F=F+1;
subplot(2,2,1);
plot(A(:,1));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度x/g');
title('滤除重力加速度');
subplot(2,2,2);
plot(A(:,2));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度y/g');
subplot(2,2,3);
plot(A(:,3));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度z/g');
for j=1:3
winS=2;
winE=winS+winSize-1;
isMotion=0;
zero=mean(A(1:10,j));
Eamax=max(abs(A(1:10,j)-zero));
threshold=2*Eamax;
while 1
[winMax,index]=max(abs(A(winS:winE,j)-zero));
%记录运动开始和运动结束时刻
% std(A(winS:winE,j))
% abs(winMax)
if isMotion==0
if winMax>threshold%静止状态下,大于阈值,运动开始
isMotion=1;
Ts(j)=index+winS-1;
As=A(Ts(j),j);
A(Te(j):Ts(j),j)=0;%运动开始前加速度为0
elseif winMax>Eamax%大于幅值,刷新阈值
Eamax=winMax;
threshold=2*winMax;
end
else winMax<threshold %运动状态下,连续小于阈值,运动结束
isMotion=0;
Te(j)=index+winS;
Ae=A(Ts(j),j);
A(Te(j):length(A),j)=0;%结束运动后加速度为0
end
if winE==length(A) %所有数据均处理完毕,跳出循环
break
end
winS=winE+1;%滑动到下一个窗口
if length(A)-winE<winSize
winE=length(A);
else
winE=winS+winSize-1;
end
end
if Ts(j)==Te(j)%无运动过程,一直处于静止状态
Ts(j)=2;
Te(j)=length(A);
A(:,j)=0;
As=0;
Ae=0;
end
%加速度补偿
A(Ts(j):Te(j),j)=A(Ts(j):Te(j),j)-As;
end
figure(F)
F=F+1;
subplot(2,2,1);
plot(A(:,1));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度x/g');
title('运动加速度补偿');
subplot(2,2,2);
plot(A(:,2));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度y/g');
subplot(2,2,3);
plot(A(:,3));
grid on;
legend('A');
xlabel('时间/s');
ylabel('加速度z/g');
%计算速度
V=A;
for j=1:3
for i=Ts(j):Te(j)
V(i,j)=V(i-1,j)+A(i-1,j)*deltaT;
end
end
%速度补偿
for j=1:3
V(Ts(j):Te(j),j)=V(Ts(j):Te(j),j)-V(Ts(j),j);
for i=Ts(j):Te(j)
V(i,j)=V(i,j)-V(Te(j),j)*(i-Ts(j))/(Te(j)-Ts(j));
end
end
%计算位移
S=V;
for j=1:3
for i=Ts(j):Te(j)
S(i,j)=S(i-1,j)+V(i-1,j)*deltaT+0.5*deltaT*deltaT*A(i-1,j);
end
end
%位移补偿
for j=1:3
S(Ts(j):Te(j),j)=S(Ts(j):Te(j),j)-S(Ts(j),j);
for i=Ts(j):Te(j)
S(i,j)=S(i,j)-S(Te(j),j)*(i-Ts(j))/(Te(j)-Ts(j));
end
end
figure(F);
F=F+1;
subplot(2,2,1);
plot(1:length(S),S(:,1));
grid on;
legend('s');
ylabel('位移x/g');
title('位移');
subplot(2,2,2);
plot(1:length(S),S(:,2));
grid on;
legend('s');
ylabel('位移x/g');
subplot(2,2,3);
plot(1:length(S),S(:,3));
grid on;
legend('s');
ylabel('位移x/g');
subplot(2,2,4);
plot(S(:,1),S(:,2));
grid on;
legend('s');
title('y-x')
% dataCount=dataCount+1;
% if dataCount==winSize || length(A)-winE<winSize
% dataCount=0;
% winE=i;
% %零状态补偿
% for j=1:3
% if std(A(winS:winE,j))<0.1 && abs(mean(A(winS:winE,j)))<0.5 %静止状态
% Ts=winE;
% A(winS:winE,j)=0;
% else
% Te=winE;
% end
% end
%
% %计算速度
% for k=winS:winE
% V=[V;
% V(k-1,1)+A(k-1,1)*deltaT , V(k-1,2)+A(k-1,2)*deltaT ,V(k-1,3)+A(k-1,3)*deltaT];
% end
% %速度的累积误差补偿
% for k=winS:winE
% for j=1:3
% V(k,j)=V(k,j)-V(winE,j)*(k-winS)/(winE-winS);
% end
% end
% %计算位移
% for k=winS:winE
% S=[S;
% S(k-1,1)+0.5*deltaT*(V(k,1)+V(k-1,1)) , S(k-1,1)+0.5*deltaT*(V(k,2)+V(k-1,2)) , S(k-1,1)+0.5*deltaT*(V(k,3)+V(k-1,3))];
% end
% %位移的累积误差补偿
% for k=winS:winE
% for j=1:3
% S(k,j)=S(k,j)-S(winE,j)*(k-winS)/(winE-winS);
% end
% end
% winS=winE+1;
% end
figure(F);
F=F+1;
subplot(2,2,1);
plot(V(:,1));
grid on;
hold on;
plot(A(:,1));
legend('V','A');
xlabel('时间/s');
ylabel('速度x/g');
subplot(2,2,2);
plot(V(:,2));
grid on;
hold on;
plot(A(:,2));
legend('V','A');
xlabel('时间/s');
ylabel('速度y/g');
subplot(2,2,3);
plot(V(:,3));
grid on;
hold on;
plot(A(:,3));
legend('V','A');
xlabel('时间/s');
ylabel('速度z/g');
figure(F);
F=F+1;
stem3(S(:,1),S(:,2),S(:,3));
grid on;
xlabel('x/s');
ylabel('y/g');
zlabel('z/g');
figure(F);
F=F+1;
plot3(S(:,1),S(:,2),S(:,3));
grid on;
xlabel('x/s');
ylabel('y/g');
zlabel('z/g');
评论22