# How to calculate position from 3 dimensional acceleration data?

I discussed this problem yesterday with Pickett here. We concluded that the best solution may be just Interpolate and then NIntegrate to get the position vector. I also suggested smoothing but as Pickett commented you can lose important details with it.

How to calculate the position from acceleration data in Mathematica?

Example where moving to one direction and other directions about the same

Import["http://pastebin.com/raw.php?i=jZ57mqZT"]


or

time_tick,acc_X_value,acc_Y_value,acc_Z_value
0.008387,-7.051625,-0.432767,-6.701011
0.041984,-7.308113,-0.712300,-6.199558
0.074841,-6.989672,-1.105712,-6.235771
0.108211,-7.580313,-0.931228,-5.587518
0.141834,-7.547990,-0.979114,-5.437576
0.174273,-7.075867,-1.138783,-6.130123
0.208107,-7.554275,-0.835906,-5.692567
0.240980,-7.329661,-0.960558,-5.710225
0.274663,-7.546344,-0.690752,-5.827994
0.308129,-6.949119,-0.860447,-6.631278
0.341972,-6.716275,-0.842939,-6.727797
0.374906,-7.340585,-0.757642,-6.366709
0.408325,-6.646092,-0.905340,-6.730790
0.440989,-6.814441,-1.069348,-6.686945
0.474559,-7.092926,-0.681474,-6.726151
0.508442,-6.090618,-0.887832,-7.290455
0.541162,-7.464041,-0.712600,-6.281712
0.574674,-6.314932,-0.937214,-6.787804
0.608281,-6.961689,-0.820642,-6.581596
0.641125,-7.042347,-0.777395,-6.264054
0.674738,-6.658662,-0.941104,-6.425518
0.708979,-6.956152,-0.764526,-6.639957
0.741439,-6.618408,-0.791462,-6.837934
0.774397,-6.924129,-0.730108,-6.721961

• Its not so much that smoothing will lose details its that you don't know what the acceleration was between timestamps. You have to assume some kind of interpolation (even linear interpolation is an assumption) and each different interpolation scheme will give you a different position; that is your "model risk". – Ymareth Oct 15 '14 at 14:19
• Have you checked that? :) – Öskå Oct 15 '14 at 14:29
• Following the link ↑↑↑ you have that youtube video. – Öskå Oct 15 '14 at 14:36

Assuming that all the initial positions (x, y and z) and the initial velocities (x', y' and z') are equal to 0 you can do:

adat = Rest@Import["http://pastebin.com/raw.php?i=jZ57mqZT"];
{ax, ay, az} = Interpolation /@ (adat[[All, {1, #}]] & /@ {2, 3, 4});
{xt, yt, zt} = (x /. Quiet@First@NDSolve[{
x == 0, x' == 0,
x''[t] == #[t]
}, x, {t, First@adat[[All, 1]], Last@adat[[All, 1]]}] & /@ {ax, ay, az});

PlotLabel -> "Acceleration", Frame -> True,
FrameLabel -> {"Time", "Aceleration"}] Plot[{xt[t], yt[t], zt[t]}, {t, First@adat[[All, 1]], Last@adat[[All, 1]]},
PlotLabel -> "Positions", Frame -> True,
FrameLabel -> {"Time", "Position"}] ParametricPlot3D[{xt[t], yt[t], zt[t]}, {t, First@adat[[All, 1]], Last@adat[[All, 1]]},
PlotRange -> {{0, -2}, {-1, 1}, {0, -2}},
AxesLabel -> {"X", "Y", "Z"}] • If you switch between InterpolationOrders and Methods do you end up in (very) different final positions? – Ymareth Oct 15 '14 at 15:14

A different interpretation would be that the acceleration is always in the same direction, but that the sensor actually rotates. In that case:

{pt, px, py, pz} = Transpose@Rest@Import["http://pastebin.com/raw.php?i=jZ57mqZT"];
{pa, pθ, pϕ} = Transpose@CoordinateTransform["Cartesian" -> "Spherical", Transpose@{px, py, pz}];

θi = Interpolation[Transpose@{pt, pθ}];
ϕi = Interpolation[Transpose@{pt, pϕ}];
ai = Interpolation[Transpose@{pt, pa}];


Using NDSolve and assuming that the acceleration data is the acceleration observed in an object (as opposed to measured within the object)

hi = (h /.
First@Quiet@
NDSolve[{h == 0, h' == 0, h''[t] == -ai[t]},
h, {t, Min@pt, Max@pt}]);


So we can plot the elevation and azimuthal angle

ListLinePlot[{pθ, pϕ}, PlotRange -> {-2 π, 2 π}, FrameTicks ->{Automatic, π Range[-2, 2]}, Frame -> True, PlotLabel -> "Orientation", FrameLabel -> {"Time", "Angles"}] The acceleration

Plot[hi[t], {t, Min@pt, Max@pt}, Frame -> True,  PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"}] And the trajectory of an object accelerating as in the data

Plot[hi[t], {t, Min@pt, Max@pt}, Frame -> True,  PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"}] To visualize the orientation over time

Animate[
Graphics3D[
{
Opacity
, Red
, Arrow[{{0, 0, 0}, {1, 0, 0}}]
, Green
, Arrow[{{0, 0, 0}, {0, 1, 0}}]
, Blue
, Arrow[{{0, 0, 0}, {0, 0, 1}}]
, Opacity[0.2]
, GeometricTransformation[
Cuboid[{-(1/2), -(1/2), -(1/2)}, {1/2, 1/2, 1/2}]
, Composition @@ {RotationTransform[ϕi[t], {0, 0, 1}],
RotationTransform[θi[t], {0, 1, 0}]}
]
}, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}], {t, Min@pt, Max@pt}] Now that I figured out how to do Composition transformations just for eye candy we could mix the trajectory with the rotation angle.

WARNING I'm mixing pears and oranges here, this is misleading and the reality here depends on the details on how and where the acceleration was measured. The data can NOT come from an accelerometer on a falling cube with initial velocity zero (as in the animation) as the acceleration measured in free fall would be minimal until terminal velocity is reached. If the acceleration if measured from 'outside' then we could not measure the angles.

The angles were estimated assuming that gravity acceleration always points down, and the observed magnitude (Norm) was fairly constant and close to $\vec{g}$ assuming units of $m/s^2$. If this comes from an accelerometer, like the ones in a phone, implies that it was not falling, but just hanging.

Therefore, this animation is not a physically accurate picture, just fun with Mathematica.

Animate[
Grid[{{
Graphics3D[
{
Opacity
, Red
, Arrow[{{0, 0, 0}, {1, 0, 0}}]
, Green
, Arrow[{{0, 0, 0}, {0, 1, 0}}]
, Blue
, Arrow[{{0, 0, 0}, {0, 0, 1}}]
, Opacity[0.2]
, GeometricTransformation[
Cuboid[{-(1/2), -(1/2), -(1/2)}, {1/2, 1/2, +(1/2)}]
, Composition @@ {
TranslationTransform[{0, 0, hi[t]}]
, RotationTransform[ϕi[t], {0, 0, 1}]
, RotationTransform[θi[t], {0, 1, 0}]

}]
}
, PlotRange -> {{-2, 2}, {-2, 2}, {-5, 2}}
, ImageSize -> 300]
,
, Plot[hi[tt], {tt, Min@pt, Max@pt}, Frame -> True,
PlotLabel -> "Trajectory", FrameLabel -> {"Time", "Hight"},
ImageSize -> 600,
Epilog -> {PointSize[Medium], Point[{t, hi[t]}]}]
}}], {t, Table[ti[k], {k, 1, 24, 0.1}]}] • +1 Could you explain or link for reference about {{{Cos[θi[t]] Cos[ϕi[t]], -Sin[ϕi[t]], Cos[ϕi[t]] Sin[θi[t]]}, {Cos[θi[t]] Sin[ϕi[t]], Cos[ϕi[t]], Sin[θi[t]] Sin[ϕi[t]]}, {-Sin[θi[t]], 0, Cos[θi[t]]}} – Junho Lee Oct 16 '14 at 0:12
• @JunhoLee It came from composing the two rotations by hand, now I do it more elegantly. – rhermans Oct 16 '14 at 8:15
• Aha thanks a lot :) – Junho Lee Oct 16 '14 at 8:51
• Correct this data is measured by having the phone in hand, well-observed -- nice visualisation +1! – hhh Oct 16 '14 at 9:21