[Edit: correct the sign in one element of one of the rotInit matrices, below.]
You have (linear) acceleration (in the body frame) and gyro data (angular velocity in the body frame), at M instants. You will do
fuse = imufilter('SampleRate',Fs);
and
q = fuse(accelerometerReadings,gyroscopeReadings);
to get q, a Mx1 vector of quaternions.
Or you can do
rotMat= fuse(accelerometerReadings,gyroscopeReadings,'orientation','Rotation matrix');
to obtain a Mx3x3 array of M rotation matrices.
The man page says "imufilter fusion correctly estimates the change in orientation from an assumed north-facing initial orientation." This means the output, q or rotMat, represents the orientations of the body, relative to its initial orientation, at M instants. The manual then says "However, the device's x-axis was pointing southward when recorded. To correctly estimate the orientation relative to the true initial orientation or relative to NED, use ahrsfilter." But the man page for ahrsfilter() indicates that ahrsfilter() expects magnetometer data, which you do not have.
If you use the default reference frame, the NED frame, then you are assuming the global x,y,z axes point north, east, down respectivey. The acceleromter and gyro data come from a sensor with a right-handed x,y,z set of axes. Suppose you know the initial orientation of the sensor relative to global frame. Then you can translate the rotMat matrices from "initial orientation of the body" frame to the global frame by pre-multiplying each 3x3 matrix in rotMat by rotInit, which is a 3x3 matrix representing the initial orientaiton of the body in the global frame.
If the initial orientation of the body is aligned with the NED axes (i.e. body x,y,z axes point north, east, and down respectively), then
rotInit=eye(3);
If the body initially is oriented with X pointing east, Y pointing south, and Z pointing down, then
x0=[0;1;0]; y0=[-1;0;0]; z0=[0;0;1];
rotInit=[x0,y0,z0];
If the body initially is oriented with X pointing south, Y pointing west, and Z pointing down, then
x0=[-1;0;0]; y0=[0;-1;0]; z0=[0;0;1];
rotInit=[x0,y0,z0];
If the body initially is oriented with X pointing northeast and 45 deg up, Y pointing southeast and horizontal, and Z pointing wherever, then
x0=[0.5;0.5;-sqrt(2)]; y0=[-sqrt(2);sqrt(2);0]; z0=cross(x0,y0);
I have not used the Sensor Fusion and Tracking toolbox, so my code above is not guaranteed. Do lots of checks with simulated data to make sure the routines behave as expected when you use relatively simple simulated input data. I have used Matlab and Labview to fuse linear and rotational acceleration data, to estimate head accelerations (rotational and linear) and orientation during soccer heading, and the results were published in several papers. So I do know something about this topic, but I am not familiar with these routines. I don't think this toolbox existed when I did that work. I would have used it.
The interesting challenge with the soccer data is that the sensor is on the edge of the head (strapped to the occiput with an elastic headband). We wanted to estimate linear and rotational accelerations in the middle of the brain. That means we must translate the measured data from the edge to the center. A purely rotational acceleration at the edge will generate linear plus rotational acceleration in the center, and vice versa. The equations are complicated and are best done with quaternions, which were not well-supported in Labview at the time.
