Better numerical integration

[atfox98]

I was going through the user's manual for the JS220 and noticed on page 33 that you are using Euler's method for computing numeric integrals.

In python this is simple, for an array x and time spacing dt you can write an accumulator like this:

for i in range(1,len(x)):
    x[i] = x[i-1] + dx[i]*dt

I think I found the dx[i]*dt step in device.py, lines 159/160.

A slightly better (more accurate) method for numerical integration is the 4th Order Runge-Kutta method. You can read more about the RK methods here: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods.

I found this page too that shows the differences in accuracy in an easy to understand format: https://felix11h.github.io/blog/euler-midpoint-rk

It takes a few more cycles of the CPU to compute, but the same accumulator can be written in python like this:

for i in range(1,len(x)):
        K1 = dx[i-1] * dt
        K2 = (dx[i-1] + 0.5*K1) * dt
        K3 = (dx[i-1] + 0.5*K2) * dt
        K4 = (dx[i-1] + K3) * dt
        x[i] = x[i-1] + (K1 + 2*(K2 + K3) + K4)/6

Here's the same accumulator, but with the 3/8 method mentioned in the Wikipedia article (supposed to have smaller error coefficients):

for i in range(1,len(x)):
        K1 = dx[i-1] * dt
        K2 = (dx[i-1] + K1/3) * dt
        K3 = (dx[i-1] + K2-K1/3) * dt
        K4 = (dx[i-1] + K3-K2+K1) * dt
        x[i] = x[i-1] + (K1 + 3*(K2 + K3) + K4)/8

I was going to try this out on the hardware, but I don't have it yet and I'm not sure I'll have the time once I do. Thought opening an issue might be the best way to get this simple improvement in front of you.

[mliberty1]

Hi @atfox98 and thanks for taking the time to think about how to make Joulescopes more accurate!

The JS220 actually performs statistics computation on the instrument in the FPGA. Changing functionality is a little more complicated than a few lines of python, but still possible. Multiplication is fine, but divisions are very expensive, especially since it's for two channels (charge, energy) at 1 Msps. We could implement alternate integration methods on the host using the sample data, too.

I have used other integration methods for inertial navigation products where time alignment is critical. With charge and energy, it's the total that matters, not time alignment. Euler integration on bandwidth-limited signals is usually pretty good if you do not care about time.

If you have the time, I would be very interested in hearing what you find for the accuracy of different methods on real-world signals!

[atfox98]

I am coming from an IMU background (dead reckoning) where Euler just didn't cut it (could be up to 6 inches off after 700 feet and going ~130-150MPH). We tried lots of integration schemes - RK methods, midpoint, simpsons, euler, etc. - and the RK4 3/8 method was the best performer out of the lot of them.

The applications are definitely different, and it may not make much of a difference for Joulescope, but I still do think a change may be worth it for that little extra bit of accuracy and advertising oomph 😊 (given you have the headroom for it on an FPGA going at 1 Msps)

If I can think of a signal where the difference is profound I’ll reply back here. If you want to set up a time to talk more I’m open to that too, just shoot me an email. It’s my github username @gmail.com

[mliberty1]

I am definitely interested in investigating alternate integration methods for Joulescopes. I'll add it to the todo list, but it will likely be some time before I can personally investigate. We need to remain focused on the Joulescope UI and implementing all of the committed, new JS220 features! I am happy to support you if you want to take a close look in the meantime.

I will leave this issue open. Feel free to investigate more and post what you find!