Central difference method or Diff for trajectories?

26 visualizzazioni (ultimi 30 giorni)
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 Commenti
Bjorn Gustavsson
Bjorn Gustavsson il 10 Set 2021
(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
Mohamad Mossad il 10 Set 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.

Accedi per commentare.

Risposta accettata

Jan
Jan il 10 Set 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 Commenti
Mohamad Mossad
Mohamad Mossad il 10 Set 2021
Modificato: Mohamad Mossad il 10 Set 2021
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
Jan il 10 Set 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.

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Physics in Help Center e File Exchange

Prodotti


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by