Double Integration of Numeric Data - Acceleration -> Displacement

I have a data set of accelerlation data points over 10ms, sampled at 1MHz. I am trying to take this acceleration data and convert it to displacement data. I have a referenced graph that is attached that I am trying to get my data to match, but it seems like they did not do a cumtrapz, and only did a point by point integral.

When I tried to use Cumtrapz it seems like, A09 it does not have the return back to -5mm after a 2nd order cumtrapz. I have tried doing a point by point integration (A9(1:end-1)+A9(2:end)) * (dt/2); and this did not bring create the same graph.

I am looking for some insight on what the right function it should be/ the reason after doing a "double integral" these two graphs do not match.

Thank you for the insight.

[John D'Errico]

I know what a definite integral means. I know what an indefinite integral means. Yet, despite having spent an entire career as a mathematician, I don't know what a definite point by point integral means. If you invent your own jargon, then expect people to be confused.

What do you mean by the vectors in A9.mat and A2.mat? Yes, apparently these are accelerations, based on your comments. at least that is what you seem to think. I might suppose they are taken at some point on a chest, given the figure title. But you might need to contact the source of the data to know the meaning, and where the failure arises.

To answer what is i think your question, is if you use cumtrapz twice on those data, assuming they would be validly measured accelerations taken from a chest cavity, AND they represent a complete beath cycle, then at some point in the cycle, the result should return to roughly the start point. You don't tell us dt, but that is just a scaling.

Could it be that this is only data from a partial breathing cycle? We don't know.

load A2.mat
plot(Datashort1);

Assuming those are accelerations, taken by an accelerometer at a constant (unknown) time step, the first integral would be velocity.

However, if we look at the figure you show, that shows a negative displacement curve. Yet the accelerations start off with a positive acceleration.

So it might be possible the accelerations are sign swapped, where the accelerometer is set to measure a positive acceleration as chest compression. That may make sense, since the figure (y axis) is labeled compression. Hard to be sure. I would check with the source.

[Kyle]

I am that you are as confused as I was. This is the language my boss used, and it confused the hell out of me to. If I am using acceleration, which these data are, as this is the raw data.

These data points are two different exposure to a blast wave and it is the acceleratoin of the chest as measured by an accelerometer located on the sturnum of a cadaver. In what I was assuming to be the 2nd integral of A9 ( test group A, test number 9) The data never rebounds as shown in the reference plot, which is using the exact same raw data vector as input. And my boss believes this is because it is due to me calculating the cumlative intergral instead of the point by point integral. Which I did not follow his logic and explination.

Edit: I have updated the name of the question to reflect a conventional understanding on integration.

Answer 1

[Torsten]

Assuming x(0) = v(0) = 0 in both cases, I get

load("A2.mat")
load("A9.mat")
a2 = Datashort1;
a9 = Datashort21;
figure
plot(0:1e-3:10,[a2,a9])

v2 = cumtrapz(-a2)*1e-3;
v9 = cumtrapz(-a9)*1e-3;
figure
plot(0:1e-3:10,[v2,v9])

x2 = cumtrapz(v2)*1e-3;
x9 = cumtrapz(v9)*1e-3;
figure
plot(0:1e-3:10,[x2,x9]*1e-2)

If v(0) is not equal to 0, things will change remarkably.

It would help if we knew the physical unit of your acceleration data.

[Kyle]

This aligns with what I got when I performed this, but as you can see in the reference graph that I attached, starting with the same raw data sets, this did not match the outputs that are in the previous graph and I am attempting to remedy this. If both the reference graph and the output that you present above are from the same accerleration data. Then what other methods would be used to calculate these deflections?

[Torsten]

The data you gave don't allow to remedy this - whatever method you choose. Always assuming v(0) = 0 !

[Walter Roberson]

Observe,

cumtrapz([3 3])

ans = 1×2

0 3

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cumtrapz() always calculates the result as-if the first output is zero. It is always taking the area relative to the baseline that is the first input point.

This is particularly important if you cumtrapz() the cumtrapz() result.

It is the difference between computing a*t^2 and a*t^2 + v*t + c . The v*t + c is important for calculating net displacement, but you cannot recover it when using cumtrapz() without taking special steps.

[Torsten]

@Walter Roberson

Isn't this exactly the condition x(0) = v(0) = 0 that allows using cumtrapz two times without correction ?

[Kyle]

I believe the chest v was 0 at t = 0 as it was steady state before the impulse was applied to it. The curx of this question was if I were to take the double derivative of the deflection presented in the png, that should give me the same output as the chest acc input. But for some reason, when I do the cumsum ( again most prominately at A9 the data does not have this recovery back to steady state that is seen in the graph presented in the png, and I was trying to discover why.

[Torsten]

Steady state doesn't mean v(0) = 0 in general. But from the initial plot in the .png where position almost remains at x = 0, it seems that x(0) = v(0) = 0 is adequate to assume.

Did you try to recover the accelation data from the displacement curve you and me deduced from the data ?

Just to check: what is the physical unit of your displacement data ?

If you add the lines

figure
a2_compare = (x2(3:end)-2*x2(2:end-1)+x2(1:end-2))/1e-6;
a9_compare = (x9(3:end)-2*x9(2:end-1)+x9(1:end-2))/1e-6;
plot(1e-3:1e-3:10-1e-3,[-a2_compare,a2(2:end-1),-a9_compare,a9(2:end-1)])

to the code above, you'll see that the second derivatives of the integrated curves perfectly coincide with the given acceleration data.

[Kyle]

Yup. Again we are in agreement here. I think the issue my PI is having with this output is that it is using the cumulative intergral, instead of a point by point. But after running a point by point, which is X+ X+1 that obviously does not even come close to the output of the cumtrapz curves.

In answer to your question the accerleration takes the form of g. But all that will do will change the scaling of the y axis and not affect the shape of the displacement graph. As far as I am conserned at this point, the reference graph is wrong and the data I have ploted here, when using these data as an input are the correct curves.

[Torsten]

In answer to your question the accerleration takes the form of g. But all that will do will change the scaling of the y axis and not affect the shape of the displacement graph.

I only noticed that the unit for g used in your acceleration data is quite strange. If the order of magnitude for the displacement compared to Picture1.jpg should be approximately equal, the unit is mm*10^(-2)/(10^(-3)*s)^2 .

But after running a point by point, which is X+ X+1 that obviously does not even come close to the output of the cumtrapz curves.

You can use a "point-by-point" method to integrate, but then you must add the results - and this is exactly what "cumtrapz" does:

load("A2.mat")
load("A9.mat")
a2 = Datashort1;
a9 = Datashort21;
figure
plot(0:1e-3:10,[a2,a9])

v2(1) = 0;
x2(1) = 0;
for i = 2:numel(a2)
  v2(i) = v2(i-1) - 1e-3*(a2(i-1)+a2(i))/2;
  x2(i) = x2(i-1) + 1e-3*(v2(i-1)+v2(i))/2;
end
v9(1) = 0;
x9(1) = 0;
for i = 2:numel(a9)
  v9(i) = v9(i-1) - 1e-3*(a9(i-1)+a9(i))/2;
  x9(i) = x9(i-1) + 1e-3*(v9(i-1)+v9(i))/2;
end
figure
plot(0:1e-3:10,[x2;x9]*1e-2)