It's possible that your confusion stems from the fact that there are multiple solutions to the problem. While your accelerometer can tell you which way is up, it cannot distinguish between North and West. I.E. If you rotate the device about the vertical axis, the outputs of the accelerometers won't change.
How can you distinguish between North and West? The best way is probably to add a digital compass to the mix. Alternatively, you may not care to know the difference between real North and real West. You may only want two orthogonal horizontal axes. I'll assume the latter.
Define our frames first. The device's frame is (X, Y, Z). The earth's frame is (V, H1, H2).
Lets assume your accelerometer readings (Ax, Ay, Az) are in the range -1 .. +1, where +1 means straight up. Immediately you know which direction is up: it's simply (Ax, Ay, Az). But how do we obtain a horizontal axis? There's a function called the Cross Product, which takes any two vectors as inputs, and returns another vector at right angles to both. Therefore, we can easily find a vector at right angles to Up. In C:
Vector3D V = (Ax, Ay, Az);
Vector3D H1 = RANDOM_VECTOR;
Vector3D H2 = CrossProduct(V, H1);
H1 = CrossProduct(V, H2);
Normalise(H1);
Normalise(H2);
So, now we have a vertical vector V, and two horizontal vectors H1 and H2. Now we just need to rotate the gyroscope readings into the same frame.
Let's call the gyroscope readings (Gx, Gy, Gz). We're going to convert them into earth frame rotation coordinates (GV, GH1, GH2). All you have to do is to think about the gyro readings like a single 3D vector. Which way is it pointing in the device's frame. Which way is it pointing in the Earth's frame?
The answer is simply:
GV = (Gx*V.x) + (Gy*V.y) + (Gz*V.z);
GH1 = (Gx*H1.x) + (Gy*H1.y) + (Gz*H1.z);
GH2 = (Gx*H2.x) + (Gy*H2.y) + (Gz*H2.z);
(I hope that's right)...
Non-zero rates are normal for MEMS accelerometers and gyros. This is persistent error. It is eliminated by somehow making sure that the device is stationary for a couple of seconds (so the output can stabilize), then getting a reading. Henceforth, subtracting this steady-state error from all future measurements. Look up the datasheet of your sensor - there will be maximum values for this and other types of measurement tolerances.
Now, the much more complex subject of fusing the accelerometer, gyro and compass data. This can get hugely complicated, using Kalman filter, like Apolo once did. It can, however, be quite simple as well.
The general idea is that the magnetic sensor has slow response, low accuracy, but the error does not increase. On the other hand, a gyro's output is velocity, which is integrated to get angular position. The error grows very fast - generally you can't do dead reckoning for more than a minute with only a giro. The accelerometer is worse - it outputs acceleration, which gets integrated twice!
So, a simple fusing filter would be some linear combination of the readings of the accelerometer and compass, with the coefficient in front of the gyro descreasing over tyme.
Here is a discussion by much more knowledgeable people than me on the topic.
Note: What you are trying to do is called dead reckoning.
Best Answer
The IMU should have the same axes for the accelerometer, gyroscope and whatever else there is inside. Just look it up in the datasheet.
My approach to this would be to first make sure that the axes of your IMU are alligned with the axes of your vehicle since the IMU can only measure the velocity, euler angles and so on of itself.
To get the euler angles from your IMU you can use this:
https://robotics.stackexchange.com/questions/4511/euler-angles-from-9dof-imu
assuming your IMU has a Magnetometer. Otherwise there are many tutorials on how to get the roll, pitch and yaw angle from IMUs:
https://engineering.stackexchange.com/questions/3348/calculating-pitch-yaw-and-roll-from-mag-acc-and-gyro-data
After this, you already have the euler angles you need and you can rotate the vectors you want (acceleration and angular velocity) into the earth frame. This can be done with quaternions as you mentioned. For this there are also a lot of tutorials online. The madgwick paper where you got the algorithm from has a good explanation for this problem.
So to clarify: You calculate the orientation of your reference frame relative to your world frame (Euler angles) with your IMU. Then you can rotate whatever vectors your IMU gives you (acceleration, angular velocity, ...) from your reference frame into the earth frame and use them (calculate trajectory, compensate gravity, calculate movement etc.)