There are different methods for numerical differentiation. Assume you have the positions store in x and the times in t. Then:
- Forward difference: v(i) = (x(i) - x(i+1)) / (t(i) - t(i+1))
- Backward difference: v(i) = (x(i) - x(i+1)) / (t(i) - t(i+1))
- Two sided difference: v(i) = (x(i+1) - x(i-1)) / (t(i+1) - t(i-1))
- 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.