Velocity from X,Y,Z position coordinates

Hello,

I have tried to compute velcocity data from X,Y,Z position coordinate data. I additionally have the time each X,Y,Z sample was recorded at, so was able to work out the time step. I am hoping someone can tell me if what I have done is correct and also I am looking for advice on what to do next. I believe I should apply a butterworth filter but I am unsure how to do this.

data is F1 (3900,3), where F1(:,1) = X, F1(:,2) = Y , F1(:,3) = Z . Time between each z,yz recording is = dt

F1=Data(:,3:5); %select data
dt = Data(2,2); % time 

for i = 1:length(F1)-1
    f1_length(i) = sqrt((F1(i+1,1) - F1(i,1))^2 + (F1(i+1,2) - F1(i,2))^2 + ...
                        (F1(i+1,3) - F1(i,3))^2);
end
 
F1_velocity = f1_length/dt;

the above code works but not sure how to do filter. I have seen this online, but not sure why f_c = 4 . I think frequency would be my dt(?) but should I also divide by 2?

% low pass filter to extract envelope
f_c = 5; % cut-off frequency
[b, a] = butter(2, f_c./(freq/2), 'low'); % Butterworth LP filter
vel = filtfilt(b, a, vel); % zero-phase filtering (effective order = 4)

finally, I would like to calculate the "area" but not sure how to do this?

[Mathieu NOE]

hello Melissa

• point 1 : yes your way of calculating the velocity is correct , but in the case you would use other (refined) ways of processing derivatives of the position data, you would calculate first the 3 velocity vectors (vx, vy, vz) and after only do the rms averaging (v_avg = sqrt(vx² + vy² + vz²))
• see attached function and demo code for first and second derivatives computations (demo_deriv.m, firstsecondderivatives.m)
• point 2 : the code is correct; if you test it with different values for f_c , you see that the lower f_c the "smoother" the output will be , but you may have to do a compromise between how much smoother it should be without compromising transients that may be of some interest (unless you are only interested in the mean value, so either very low f_c or use directly mean)
• sampling freq = 1/dt
• in case you would do fft / spectrogram analysis , see code attached (simple_demo_wav.m) as example

hope it helps

[melissa sandison]

Hello Mathieu,

Firstly, thank you so much for taking the time to help with my code. I really appreciate it.

I wen't through your code and I understand now how you have obtained velocity profile using forward/central/backward difference. This is most helpful!

Thank you also for the simple demo wav example, I know understand (and I read up on nyquist freq and also allaising).

Honestly, cant thank you enough! all the best!

[Mathieu NOE]

Glad it helped !