Plotting 1st derivative and 2nd derivative graph from a set of values

I have a set of raw data collected from a displacement time graph and i wish to convert it into a velocity time graph and acceleration time graph. Is there anyway to do that? My data is the time at each interval of displacement so I do not have enough points to plot using the (y1-y2)/(x1-x2) method to extrapolate the gradient. Is there anyway to get the derivatives from the line of best fit? I cannot just differentiate line of best fit polynomial as it becomes a straight line graph after 1.5s so the best method is to find gradient of this graph at many points and plot from there.

Data points:

Answer 1

[Mathieu NOE]

hello

below a useful code for first and second derivative - to apply directlty on your raw data or on fitted / interpolated data (if you want something smoother)

%%%%%%%%%%%%%%%%%%%%%%
function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.

% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n = length (x);
dy = zeros;
ddy = zeros;

% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.

% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.

dy(1) = (-3*y(1) + 4*y(2) - y(3)) / 2*(x(2) - x(1)); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

for i = 2:n-1
   
dy(i) = (y(i+1) - y(i-1)) / 2*(x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;

end

dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / 2*(x(n) - x(n-1));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;

end

[Zhi Xiang Teo]

Hello Mathieu, thank you for the help, here is the results that I got after tabulating the raw data into the code you've provided. Can I check if this is the only set of data I should get?

Also, this is the graph that I got after plotting the values in excel. However, it does not seem to be the graph that I am looking for. Can I ask what do the values represent?

[Feng Du]

Hello Mathieu, thank you for your sharing, your codes is very useful.

I don't quite understand derivatives, can you tell me why these two formulas are used?

dy(1) = (-3*y(1) + 4*y(2) - y(3)) / 2*(x(2) - x(1)); % First derivative

ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

[Mathieu NOE]

hello @feng du

basically it's explained in the comments section - for more info's see the attached paper

all the best

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

[Feng Du]

How nice you are!

Thank you!

[Mathieu NOE]

my pleasure !

Answer 2

[Mathieu NOE]

hello

FYI, you have plotted time vs speed (and not speed vs time) in your excel spreadsheet

in your first plot , also it is not displacement vs time , it's reversed.

otherwise , see below for code implementation of first and second derivative. The curves may not be very smooth if you apply it directly on the raw data. You can do spline or sliding average smoothing if you are in need to smooth the curves.

see below some code for you - hope it helps !!

t = 0:5:85;
s = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];

[velocity, accel] = firstsecondderivatives(t,s);

% raw data + first & second derivatives plot
figure(1),
subplot(3,1,1),plot(t,s,'b');
legend('disp');
subplot(3,1,2),plot(t,velocity,'r');
legend('velocity');
subplot(3,1,3),plot(t,accel,'m');
legend('accel');

% do some interpolation first 
tt = t(1):1:t(end); % refined time axis

% option 1A : spline interpolation
ss = spline(t,s,tt);

% option 1B : sliding avg smooting
N = 20;
ss = myslidingavg(ss, N);

[velocity2, accel2] = firstsecondderivatives(tt,ss);

% interpolated + smoothed data + first & second derivatives plot
figure(2),
subplot(3,1,1),plot(tt,ss,'b');
legend('disp');
subplot(3,1,2),plot(tt,velocity2,'r');
legend('velocity');
subplot(3,1,3),plot(tt,accel2,'m');
legend('accel');


function out = myslidingavg(in, N)
%   OUTPUT_ARRAY = MYSLIDINGAVG(INPUT_ARRAY, N)
%
%  The function 'slidingavg' implements a one-dimensional filtering, applying a sliding window to a sequence. Such filtering replaces the center value in
%  the window with the average value of all the points within the window. When the sliding window is exceeding the lower or upper boundaries of the input
%  vector INPUT_ARRAY, the average is computed among the available points. Indicating with nx the length of the the input sequence, we note that for values
%  of N larger or equal to 2*(nx - 1), each value of the output data array are identical and equal to mean(in).
%
%  *  The input argument INPUT_ARRAY is the numerical data array to be processed.
%  *  The input argument N  is the number of neighboring data points to average over for each point of IN.
%
%  *  The output argument OUTPUT_ARRAY is the output data array.

if (isempty(in)) | (N<=0)                                              % If the input array is empty or N is non-positive,
 disp(sprintf('SlidingAvg: (Error) empty input data or N null.'));     % an error is reported to the standard output and the
 return;                                                               % execution of the routine is stopped.
end % if

if (N==1)                                                              % If the number of neighbouring points over which the sliding 
 out = in;                                                             % average will be performed is '1', then no average actually occur and
 return;                                                               % OUTPUT_ARRAY will be the copy of INPUT_ARRAY and the execution of the routine
end % if                                                               % is stopped.

nx   = length(in);             % The length of the input data structure is acquired to later evaluate the 'mean' over the appropriate boundaries.

if (N>=(2*(nx-1)))                                                     % If the number of neighbouring points over which the sliding 
 out = mean(in)*ones(size(in));                                        % average will be performed is large enough, then the average actually covers all the points
 return;                                                               % of INPUT_ARRAY, for each index of OUTPUT_ARRAY and some CPU time can be gained by such an approach.
end % if                                                               % The execution of the routine is stopped.

out = zeros(size(in));         % In all the other situations, the initialization of the output data structure is performed.

if rem(N,2)~=1                 % When N is even, then we proceed in taking the half of it:
 m = N/2;                      % m = N     /  2.
else                           % Otherwise (N >= 3, N odd), N-1 is even ( N-1 >= 2) and we proceed taking the half of it:
 m = (N-1)/2;                  % m = (N-1) /  2.
end % if

for i=1:nx,                                                 % For each element (i-th) contained in the input numerical array, a check must be performed:
  dist2start = i-1;         %  index distance from current index to start index (1)
  dist2end = nx-i;          %  index distance from current index to end index (nx)
  if dist2start<nx | dist2end<nx        % if we are close to start / end of data, reduce the mean calculation on centered data vector reduced to available samples
      dd = min(dist2start,dist2end);    % min of the two distance (start or end)
      out(i) = mean(in(i-dd:i+dd));     % mean of centered data , reduced to available samples
 end % if
end % for i

[Ng Joshua]

Hi thank you! its a very helpful code. However, I just realised that I posted my values wrongly the time and displacement column were swapped. If I were to use this code could I just change the data at the start like this:

s = 0:5:85;

t = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];

and run it accordingly?

[Mathieu NOE]

hello

I believe only the display axes are swapped. in excel you should put the time column first.

The code I submitted to you assume that t is time and s is your data. See the functions arguments (t = time , s = data).

so I don't believe that doing this (below) will give the expected results (IMHO)

s = 0:5:85;

t = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];

[Ng Joshua]

Hi i have been trying your code out but I am getting an error message:

Unrecognized function or variable 'firstsecondderivatives'.

So I was wondering if it was an imported function?

[Mathieu NOE]

hello

yes here in attachement

sorry I forgot that !!

all the best

[Zhi Xiang Teo]

[Zhi Xiang Teo]

Hi Mathieu,

We swapped the t and s matrix and the x and y in the function firstsecondderivatives but met an error in line 21 of the function. The graph of displacement and velocity looks great, any idea what went wrong? Also, mad appreciation for your help!

Xav

[Mathieu NOE]

Hi Xav

I think the error came after you swappd x, y in the function - this is 1/ not necessary and 2/ will of course generate an error as the function was written with x, y as the right way to enter the data (and not y , x).

And x is supposed to be the x axis data , so for y = y data (obviously). If you change that order in the function call it cannot work anymore.

so keep firstsecondderivatives(x,y); (and do no modification inside the function code)

and just make sure what correspond to x and y in your data. x = t or s ?

regards

[Zhi Xiang Teo]

Yeah I've tried that before. I edit the input as such:

s = 0:5:85;

t = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];

[velocity, accel] = firstsecondderivatives(t,s);

and the error code appear:

Index exceeds the number of array elements (3).

Error in firstsecondderivatives (line 21)

ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative

[Mathieu NOE]

hi

I have no error on my side when I run the same code as you

I got first and econd derivative computed ok :

Name Size Bytes Class Attributes

accel 1x18 144 double

s 1x18 144 double

t 1x18 144 double

velocity 1x18 144 double

Have you changed anything in the firstsecondderivatives function ?

I ask because your line 21 is not my line 21 - see below (my version)

function [dy, ddy] = firstsecondderivatives(x,y)

% The function calculates the first & second derivative of a function that is given by a set
% of points. The first derivatives at the first and last points are calculated by
% the 3 point forward and 3 point backward finite difference scheme respectively.
% The first derivatives at all the other points are calculated by the 2 point
% central approach.
% The second derivatives at the first and last points are calculated by
% the 4 point forward and 4 point backward finite difference scheme respectively.
% The second derivatives at all the other points are calculated by the 3 point
% central approach.

n = length (x);
dy = zeros;
ddy = zeros;
% Input variables:
% x: vector with the x the data points.
% y: vector with the f(x) data points.
% Output variable:
% dy: Vector with first derivative at each point.
% ddy: Vector with second derivative at each point.
dy(1) = (-3*y(1) + 4*y(2) - y(3)) / 2*(x(2) - x(1)); % First derivative
ddy(1) = (2*y(1) - 5*y(2) + 4*y(3) - y(4)) / (x(2) - x(1))^2; % Second derivative
for i = 2:n-1
   
dy(i) = (y(i+1) - y(i-1)) / 2*(x(i+1) - x(i-1));
ddy(i) = (y(i-1) - 2*y(i) + y(i+1)) / (x(i-1) - x(i))^2;
end
dy(n) = (y(n-2) - 4*y(n-1) + 3*y(n)) / 2*(x(n) - x(n-1));
ddy(n) = (-y(n-3) + 4*y(n-2) - 5*y(n-1) + 2*y(n)) / (x(n) - x(n-1))^2;
end

[Zhi Xiang Teo]

I didn't change anything in in the firstsecondderivatives function. The same problem arise when i make

s = 0:5:85;

t = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];

even when i copied the function you posted above.

[Mathieu NOE]

hello

this is very strange

I tested it again in both cases , and it works by me fine :

s = 0:5:85;
t = [0 0.5 0.72 0.89 1 1.2 1.35 1.41 1.52 1.63 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5];
[velocity, accel] = firstsecondderivatives(t,s);        % it works !!
[velocity2, accel2] = firstsecondderivatives(s,t);      % it works !!

now on your side, did it never work ? it seems that the first time you used it , it was fine , or ?

so what hapenned ? it's a mystery for me ... it never failed on my matlab (tested even diffferent releases)