Assuming the coordinate frames are as follows:
ECI, the world frame
BODY1, the IMU1 body frame
BODY2, the IMU2 body frame
Examining your quaternion arithmetic:
q_start1 = mean(Q1(t0->5s)); % ECI->BODY1_t0
q_start2 = mean(Q2(t0->5s)); % ECI->BODY2_t0
% (NOTE: The above should be normalized)
Now in order to express their frames in world frame:
Q1_w = Q1 * Inverse(q_start1); % ECI->BODY1_t * inverse(ECI->BODY1_t0) = ECI->BODY1_t * BODY1_t0->ECI
Q2_w = Q2 * Inverse(q_start2); % ECI->BODY2_t * inverse(ECI->BODY2_t0) = ECI->BODY2_t * BODY2_t0->ECI
And to obtain the relative rotation between them:
Q_relative_w= Inverse(Q1_w) * Q2_w; % inverse(ECI->BODY1_t * BODY1_t0->ECI) * (ECI->BODY2_t * BODY2_t0->ECI)
Expanding that last one:
(ECI->BODY1_t0 * BODY1_t->ECI) * (ECI->BODY2_t * BODY2_t0->ECI) =
ECI->BODY1_t0 * (BODY1_t->ECI * ECI->BODY2_t) * BODY2_t0->ECI =
ECI->BODY1_t0 * (BODY1_t->BODY2_t) * BODY2_t0->ECI
So, you have what I think you are after, BODY1_t->BODY2_t, wrapped inside of other ECI to BODY stuff. I think what you should be after is just the BODY1_t->BODY2_t stuff. E.g.,
q_t0 = inverse(q_start1) * q_start2; % inverse(ECI->BODY1_t0) * ECI->BODY2_t0 = BODY1_t0->BODY2_t0
Then calculate this downstream:
q_t = inverse(Q1) * Q2; % inverse(ECI->BODY1_t) * ECI->BODY2_t = BODY1_t->BODY2_t
And see if all of the q_t calculations seems to be fairly constant.
For a discussion of MATLAB quaternion conventions (you need to know this in order to do the quaternion arithmetic correctly), see these links:
Best Answer