Central difference method or Diff for trajectories?

Hello,

I have a dataset (2000 to 5000 sample size) of positions for (30 to 100) trajectories in 3D obtained from tracking particles in images with an exact known time step or frame rate. The problem when calculating the velocities of the particles is that of course if I use the "diff" function, the velocity array would have a smaller size than the positions, which is not what I want since I use both the velocities and positions in the same function to calculate another thermodynamic quanitity.

So I learned about the central difference method/gradient which gives a same-sized velocity array, the question is:

Is that correct to use it here? and which one since there's "gradient", "Dgradient" and "central diff". I don't know if that's relevant but the velocities should be normally distributed. I would also like to know if there's a preferred normality test for such data and sample size because most of the defined functions for the tests usually give different results.

I would really appreciate the help.

Thank you.

[Walter Roberson]

Jan's DGradient in the File Exchange https://www.mathworks.com/matlabcentral/fileexchange/29887-dgradient is optimized and supports more methods. However it is not necessarily more accurate.

The usual test for normal distribution is KS Test https://www.mathworks.com/help/stats/kstest.html

[Mohamad Mossad]

Thank you for the info,

Do you mean it's not necessarily more accurate than the Diff function or the default gradient?

[Walter Roberson]

Not necessarily more accurate than gradient() when you use the (default) 1stOrder option.

I do not know anything about the mathematics of the 2ndOrder option.

[Bjorn Gustavsson]

If you are talking about trajectories for particles that are accelerated by some forces you might consider what your velocities are supposed to represent. Some high-accuracy numerical methods for integrating equations of motion use leapfrog-type schemes where the velocities are estimated at the "+1/2-timestep"s and not explicitly at the same time-steps of the positions. You might be interested in reading up on the Størmer-Verlet method.

[Mohamad Mossad]

Thank you Bjorn,

I do know about the Velocity Verlet algorithm, having used it before in an MD simulation but it hadn't occured to me to use it on this experimental data. I use a random laser force kicks to "heat" or accelerate dust particles in a plasma, that's why I'm interested in studying the velocity distribution of the particles. Do you think the Verlet algorithm fit such data? The particles should maintain a normal "Maxwellian" distribution before the "heating" or acceleration process for an ideal plasma in equilibrium.

[Bjorn Gustavsson]

(Since I work in auroral physics in Norway (which is what Størmer worked on and were from respectively) there is no way I will ever acknowledge the VV-algorithm named as such - but others might find me petty on this point). My idea was mainly that the velocities estimated with diff(r)/dt would naturally fit with an interpretation based on Størmer-Verlet-integration as the average velocity half-way between the different time-steps where you have the positions determined. You "just" have to remember to shift the v-estimates half a time-step forward. Exactly how it works on your experiment I have no idea - but I'm curious to see. Just go ahead and test!

[Mohamad Mossad]

Oh, I had no idea, I've always read, heard and known it as the VV-algorithm, you do have the right to correct me for using it, although it is the commonly used name. I always felt rightful attribution is a problem in our field, similar to the Ibn Sahl and Snellius. Regarding applying it on exp. data, I will consider that and see what I get. Thank you for the hint.

Answer 1

[Jan]

There are different methods for numerical differentiation. Assume you have the positions store in x and the times in t. Then:

1. Forward difference: v(i) = (x(i) - x(i+1)) / (t(i) - t(i+1))
2. Backward difference: v(i) = (x(i) - x(i+1)) / (t(i) - t(i+1))
3. Two sided difference: v(i) = (x(i+1) - x(i-1)) / (t(i+1) - t(i-1))
4. 2nd order 2 sided:

v(i) = ((x(i+1) * (t(i)-t(i-1)) / (t(i+1)-t(i))) -
             (x(i-1) * (t(i+1)-t(i)) / (t(i)-t(i-1)))) / (t(i+1)-t(i-1))
            + x(i) * (1.0 / (t(i)-t(i-1)) - 1.0 / (t(i+1)-t(i)))            

For evenly spaced data, so all t are equidistant, 4. is the same as 3.

The 2-sided difference is more accurate and less sensitive for noise.

5. If you have a model for the position, you can fit the parameters of the model to the data (e.g. a least squares fit) and reduce the noise coming from the measurement of the position. Then you can determine the velocity of the fitted model. Example: Throwing a heavy sphere in the gravity and ignoring the friction, which is valid as long as the velocity is low. Then you can assume that the trajectory is a parabola and the speed grows linearily.

6. Without having an explicit and validated model, you can assume some properties of the derivatives, e.g. that the higher order derivatives vanish, or in other words the position can be described as polynomial with a small number of terms. For objects moving in the air with friction you can assume, that the forces (which change the acceleration) depend on the position and the velocity only, but not on the acceleration or the jerk. This mean, that the position can be modelled locally by a polynomial of a certain order to reduce noise. The Savitzky-Golay-Filter fits a polynomial to the local environment. e.g. 10 neighboring points. Then you can calculate the smoothed derivative at the center of this polynomial. See e.g. https://www.mathworks.com/matlabcentral/fileexchange/3514-savitzky-golay-smoothing-filter

The difference between the methods is the handling of the noise. The two-sided difference is less sensitive for noise, but has a lower resolution. Imagine this x:

x = [1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]

While the one-sided difference finds a high frequency change of the speed, the two-sided difference claims, that the speed is 0.

If you have a model and adjust the parameters by a method for fitting, this process cares for the noise of the measured positions. If you know, that x is the position of a spring pendulum, you have measured the positions with the eigenfrequency, you can easily designe a sine wave and determine the velocity with subpixel resolution based on the given x.

Using the version 6 with fitting a simple general model like a polynomial to the data, you suppress high frequency noise in the velocities. Theís is critical, if the assumptions about the motion is wrong and you remove important information.

There is no general "best" method. It depends on the physical nature of the process and on what you want to measure. Imagine you measure the sea level with a frequency of a second with a laser. If you want to calculate the average height of the waves, you need another smoothing approach that if you want to determine the tides. Calculating the derivatives of noisy data has a strong relation to filtering the data.

[Walter Roberson]

@Jan

Your formula for forward and backward difference are the same

[Mohamad Mossad]

Thank you Jan,

Your explanation is very detailed and informative. I already applied option 6, because I can't model the trajectories of particles. I used the mentioned Savistky Golay filter on the positions and used a window interval of 7 frames for a polynomial of 2nd order, since I was using a normality test (as an indicator for the goodness of the window length) to see how many particles whose velocities are normally distributed (particles are in plasma, so their velocities should follow Maxwell-Boltzmann distribution, which is a normal distribution). Of course such tests aren't always the best choice, and QQ plots can be better when it comes to large data sizes. Having used both diff() (forward) and gradient() (central), they almost give the same number of normally distributed velocities of particles, although gradient() gives fewer, the velocity array has the same length as the positions array which is benefecial for me.

Edit: @Jan I had already applied both the diff() and gradient() after the Savitsky-Golay filter, I didn't know about your attached function to calculate the smoothed derivative on FileExchange. How different is it from the default filter on MATLAB, except the derivative?

[Jan]

Matlab's sgolayfilt determines the polynomial and replaces the midpoint by the value of the polynomial. Afterwards you can use diff or gradient or my DGradient to obtain the derivative. The linked function savitzky-golay-smoothing filter does both in one step. If you have the parameters of the polynomial, it is cheap to determine the value of the derivatives directly. This avoids some rounding errors also.