6
$\begingroup$

I'm given a stream of acceleration data with timestamps. The sampling is non-uniform.

Apart from Euler, is there a way to integrate the acceleration into velocity? Something more accurate or of higher order?

I can store some past acceleration data, but I'd like to output the velocity as quickly as possible, as new data continually arrives.

Improve this question
$\endgroup$
2
  • $\begingroup$ Migrating per OP request. $\endgroup$
    – Willie Wong
    Commented May 27, 2014 at 11:16
  • $\begingroup$ Warning: If you have data from accelerometers, it will be given in the body frame. Unless your platform is entirely non-rotating, you will also need to take into account changes in rotation. This is the well-known INS Mechanisation Equations, of which there are a large number of sub-optimal solutions. $\endgroup$
    – Damien
    Commented May 30, 2014 at 8:57

4 Answers 4

Reset to default
7
$\begingroup$

I think it may be more useful to think of this is as numerical integration of a series of data points rather than as the solution of an ODE. Adams-Bashforth could work as suggested by @Omnomnomnom, but I think there are better methods.

It seems that your acceleration values are "given" meaning that you have no control over the time step or the times at which the acceleration is known. In contrast to solving ODEs, this also means that you do not need to project forward in time to get future velocities to enable the next computation, you can simply wait until you have the new acceleration value to perform the next computation.

In this context, there are several basic integration methods that would work, but I think the simplest will be based on directly integrating piecewise polynomial approximations of your acceleration values. The simplest example of this is the well-known Trapezoid Rule which is just integration of a piecewise linear function and gives second order accuracy. If you want higher-order approximations of the data, you can use higher-order polynomials.

If it's possible to delay computing $v(t_i)$ by one timestep, I would suggest waiting until you know $\{a(t_{i-2}),a(t_{i-1}),a(t_{i}),a(t_{i+1})\}$ and then using a cubic polynomial on these four data points to integrate from $t_{i-1}$ to $t_i$. Using Lagrange interpolating polynomials, this can be done as follows: $$ v(t_i) = v(t_{i-1}) + \int_{t_{i-1}}^{t_i}a(t)dt $$ where you can approximate $a(t)$ as \begin{align} a(t) = & a(t_{i-2})\dfrac{(t-t_{i-1})(t-t_{i})(t-t_{i+1})}{(t_{i-2}-t_{i-1})(t_{i-2}-t_{i})(t_{i-2}-t_{i+1})} \\ & + a(t_{i-1})\dfrac{(t-t_{i-2})(t-t_{i})(t-t_{i+1})}{(t_{i-1}-t_{i-2})(t_{i-1}-t_{i})(t_{i-1}-t_{i+1})} \\ & + a(t_{i-2})\dfrac{(t-t_{i-2})(t-t_{i-1})(t-t_{i+1})}{(t_{i}-t_{i-2})(t_{i}-t_{i-1})(t_{i}-t_{i+1})} \\ & + a(t_{i+1})\dfrac{(t-t_{i-2})(t-t_{i-1})(t-t_{i})}{(t_{i+1}-t_{i-2})(t_{i+1}-t_{i-1})(t_{i+1}-t_{i})} \end{align}

You can then either integrate this function analytically or use a quadrature formula. Two-point Gauss-Legendre quadrature is exact for third order polynomials so that would work.

Edit:
Depending on the language you're working with, there are likely available tools for polynomial fitting, e.g. MATLAB's polyfit that would make this implementation very easy.

  1. Get the polynomial coefficients $y = at^3+bt^22+ct+d$
  2. Use the definite integral, $\int_{t_1}^{t_2}y(t)dt=\left.\left(\frac{a}{4}t^4+\frac{b}{3}t^3+\frac{c}{2}t^2+dt\right)\right|_{t_1}^{t_2}$
Improve this answer
$\endgroup$
8
  • 2
    $\begingroup$ This is what I thought of as well, however one thing I want to point out is that if there is high frequency noise in the acceleration signal, fitting higher orders may introduce extra errors since the polynomials would try to fit the local noise rather than the underlying function. $\endgroup$
    – Godric Seer
    Commented May 27, 2014 at 18:21
  • $\begingroup$ Agreed. If noise is a concern some form of smoothing spline might be a good choice. $\endgroup$
    – Doug Lipinski
    Commented May 27, 2014 at 18:28
  • $\begingroup$ If the noise is symmetric, aren't the oscillations induced by the interpolant also symmetric, leading to the same integrated values? Or is that not the case in general... $\endgroup$
    – Aurelius
    Commented May 27, 2014 at 18:39
  • $\begingroup$ @Aurelius in a very large majority of cases you are correct. With a cubic fit, the chance of points being just right to cause a sharp peak that throws off the integration is fairly small, but every extra order in the polynomial increases this chance significantly. I don't want to suggest that this solution is bad, just that it has a pitfall that should be watched for. $\endgroup$
    – Godric Seer
    Commented May 27, 2014 at 18:47
  • $\begingroup$ Interesting, I was just thinking out loud, I'm not terribly well-informed on the precise form of the spurious oscillations that can get induced by interpolation - it just intuitively seemed that those too might also be symmetric leading to the same evaluation of the integral. Integration has a nice way of smoothing out data. $\endgroup$
    – Aurelius
    Commented May 27, 2014 at 18:52
2
$\begingroup$

Since you have all the data points already, I agree with @Omnomnomnom that you might want to try numerical quadrature here - a simple trapezoidal rule will be second-order accurate as compared to Euler, which is first-order. (This is equivalent to using Euler for integration, but using the average acceleration between two times as your acceleration value for the step, which should also be 2nd order.)

The easiest higher order approach that comes to mind is to fit some curve to it (e.g. a cubic spline) and that gives you an approximating polynomial that can be integrated "exactly" by quadrature or other means, or you can use any kind of ODE solver like an adaptive Runge-Kutta method. Doing this should be 2 or 3 lines of matlab or numpy so it's a good first effort to see if it's sufficiently accurate, or if you need to use higher order interpolations.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Note that standard cubic splines depend on all data points so adding new data means re-computing everything. You could get around this by using Catmull-Rom cubic splines or just forming a cubic polynomial based on the nearest 4 points. $\endgroup$
    – Doug Lipinski
    Commented May 27, 2014 at 17:57
  • $\begingroup$ Yeah, I had written this before your post; I like your answer much better! My mind immediately went to the "how would I do this in 3 lines of matlab" mode. $\endgroup$
    – Aurelius
    Commented May 27, 2014 at 18:07
  • $\begingroup$ There's definitely something to be said for ease of implementation, especially for testing. I guess polyfit could be used directly in MATLAB to implementing my answer as well. $\endgroup$
    – Doug Lipinski
    Commented May 27, 2014 at 18:11
  • $\begingroup$ Yeah, though polyfit is fixed order for set points, so you'd have to loop/iterate through regions. I was thinking this: pp=csapi(time,accel); [t,y] = ode45(@(t,y) fnval(pp,y),[t0 tf],y0); gets it done quick and dirty in 2 lines (or use integrate instead of ode45). $\endgroup$
    – Aurelius
    Commented May 27, 2014 at 19:22
1
$\begingroup$

The Adams-Bashforth methods are higher-order extensions of Euler's method, which seems to fit what you're looking for.

Since your situation is a particularly simple ODE, however, you may want to consider other methods, including various methods of numerical quadrature.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Is there any specific method I should concentrate on? Skimming the pages, I couldn't find anything that relates to non-uniform, unknown in advance, time intervals. $\endgroup$
    – user1071136
    Commented May 26, 2014 at 21:41
  • $\begingroup$ I thought the adaptive algorithms section might have been useful but I think that's the solution for another problem... I suppose it would be best to try going forward with one of the Adams-Bashforth methods $\endgroup$
    – Omnomnomnom
    Commented May 26, 2014 at 21:51
0
$\begingroup$

You can use the trapezium rule as said by Doug in this way (in MATLAB):

cumtrapz(time,data)

which is the cumulative trapezoidal integration. The 2 inputs are vectors having the same lenght. In this way you can for example obtain the velocity integrating the acceleration. The output is a vector with the same length of the inputs.

Or you can use the Simpson's rule as follows: use quadl() to integrate your data but first you need to create a function in which you interpolate them.

function f = int_fun(x,xdata,ydata)
f = interp1(xdata,ydata,x);

And then feed it to the quadl() function:

integral = quadl(@int_fun,A,B,[],[],x,y) % syntax to pass extra arguments
                                         % to the function
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.