MATLAB Answers

Central difference method or Diff for trajectories?

11 views (last 30 days)
Mohamad Mossad
Mohamad Mossad on 10 Sep 2021
Commented: Jan on 10 Sep 2021
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.
  7 Comments
Mohamad Mossad
Mohamad Mossad on 10 Sep 2021
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.

Sign in to comment.

Accepted Answer

Jan
Jan on 10 Sep 2021
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.
  3 Comments
Jan
Jan on 10 Sep 2021
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by