Integrating discrete values to determine the absolute value of displacement

I have a velocity dataset [stored in 700 x 1 double] that i want to get the total displacement of (eg: I don't want net displacement. For total displacement, any negative displacement I would flip it to positive values). Normally, I would use cumtrapz again to gain the displacement, however any positive or negative displacements would cancel each other out. I don't want that.

I know that the area under the velocity curve is supposed to yield the total displacement. What would be the best way to code it such that the negative area is captured as positive area? Please note that I do not have points that are exactly zero (so point 1 may be +0.7 and point 2 may be -0.5, but it does not have a point in the discrete data which shows 0). I cannot simply just make the y outputs as all absolute value, because integrating between +0.7 and +0.5 yields a different area than +0.7 and -0.5. Any help would be greatly appreciated.

Answer 1

[John D'Errico]

1. Find any consecutive pair of points that change the sign of their velocity.
2. Between such a pair, use (linear) interpolation to locate the time where the velocity was zero.
3. Insert new points at those zeros.
4. Now take the absolute value of all velocities in the new expanded vector of velocities.
5. Use trapz or cumtrapz (as desired) to integrate velocity.

Break "large" problems down into small, managable peices. Then solve each one, and put it all together. A large problem is defined as any problem that is too difficult for you to solve easily.

[Coral]

Hi John,

Thanks for the speedy reply. Would you be able to share the appropriate functions/codes to perform the tasks mentioned above? I'm not exactly familiar with MATLAB and the codes associated with step 1 and 2 specifically.

[John D'Errico]

There are many ways to do this. Take it as a good opportunity to learn MATLAB. Consider the vector V:

V = rand(1,10)*2 -1

V = 1×10

0.2399 -0.1239 0.0360 -0.2833 -0.0479 -0.0424 -0.0588 0.7599 0.4511 -0.5625

So both positives and negatives.

What function returns the sign of each element in a vector? (Hint, what does sign(V) do? Why not try it?

sign(V)

ans = 1×10

1 -1 1 -1 -1 -1 -1 1 1 -1

How can you tell when the sign changes? What does diff do?

diff(sign(V))

ans = 1×9

-2 2 -2 0 0 0 2 0 -2

That is, whenever the sign is constant, diff returns 0. When the vector changes sign, then the differences between consecutive elements of the sign of V will be either +2 or -2. That is, it will be non-zero. And what does find look for?

ind = find(diff(sign(V)))

ind = 1×5

1 2 3 7 9

So there was a change in sign in every case between the indicated element, and the one immediately after. For example...

V(ind(1) + [0 1])

ans = 1×2

0.2399 -0.1239

So different signs. ind and ind+1 locates EVERY such occurrence.

If we know that V(ind) and V(ind+1) always indicate a zero crossing, then how do we locate the intermedaite point where the line segment crossed zero? Again, this is just linear interpolation.

Assume the times for our original vector are in the vector t. I'll assume V was uniformly spaced initially.

t = 0:(numel(V)-1);

I'll just do the interpolation by hand here. two points define a line segment. You will need to accept this simple formula, or you can derive it, as I did on paper offline.

tinterp = (V(ind).*t(ind+1) - V(ind+1).*t(ind))./(V(ind) - V(ind+1));
Vinterp = zeros(size(ind));

So, for example, in between those first two points, we have:

[t(ind(1)), tinterp(1), t(ind(1)+1)]

ans = 1×3

0 0.6595 1.0000

With corresponding velocities at those points:

[V(ind(1)), 0, V(ind(1)+1)]

ans = 1×3

0.2399 0 -0.1239

So it looks as if my quick and dirty linear interpolation worked correctly.

Now, the easiest way to insert those new points in the vector of velocities uses a sort.

Vnew = [V,Vinterp];
[tnew,tags] = sort([t,tinterp],'ascend');
Vnew = Vnew(tags);

Did it work? Yes. Of course.

plot(t,V,'r*-',tinterp,Vinterp,'bs')
yline(0)

So at every point where the curve crossed zero, a new point has been created that lies intermediate, and at the zero crossing.

But what you want is the curve of absolute velocities.

plot(t,V,'r*-',tinterp,Vinterp,'bs',tnew,abs(Vnew),'go-')
yline(0)

We can now trivivally integrate, as

trapz(tnew,abs(Vnew))

ans = 1.7579

The green curve there is what trapz integrated. As you can see, if is the same as the original curve when V was positive, but it mirrors the negative values for V above the axis.

Are there easier ways to have done this? Well, probably. Even certainly so. For example, I could just have pointed you to the function intersections.m, as written by Doug Schwarz, and found on the file exchange. That would have done most of the work for you.

[Coral]

Thank you so much, John. I went ahead and tested out the functions you mentioned, following the steps and trying to understand how they work. This was really informative and I really appreciate the time you've put in to not only put together the functions, but also giving a step-by-step explanation, which is much better than testing and trying random bits of code on MATLAB to try and achieve the desired results. Conceptually, it made sense but I was at a loss as to what functions to employ for this. It's all coming together now that you've explained the coding piece.

I've also just tried to apply the same code on to my dataset and I'm happy to say that it works!

Again, thank you for your help!

Edit: I'm sure there are indeed faster ways, like that function you suggested, but I'm a little old fashioned and prefer to understand what's going on before heading to the shortcuts :)

[John D'Errico]

I'll concede the use of sort there may not have been obvious. :) But once you see it done, it should make sense.

The other line of code that may not be completely obvious is the interpolation step. That comes from a line segment connecting two known points, where you know the line segment crosses y==0. Then find the location on the x-axis where the segment must have crossed the axis. As I did, just draw what is happening on paper, write the corresponding equation for the line, then do some algebra. That formula falls out. In fact, the formula applies for any two endpoints on a line, as long as the line is not horizontal. They need not be on opposite sides of the axis, but then the point of intersection will not be between the endpoints on the segment.

I guess I could have derived the formula, given two points (x1,y1), (x2,y2). The equation of that line goes back to some ancient algebra class I must have taken:

syms x1 y1 x2 y2 x y
lineeq = (y - y1) == (y2 - y1)/(x2-x1)*(x - x1)

lineeq = 

xsol = solve(subs(lineeq,y,0),x,'returnconditions',true)

xsol = struct with fields:
             x: [1×1 sym]
    parameters: [1×0 sym]
    conditions: [1×1 sym]

xsol.x

ans = 

xsol.conditions

ans = 

In that solution, we see the formula I gave, as well as the requirement that y1~=y2, otherwise the line is horizontal and never crosses the axis at all.