Problem using cumtrapz with acceleration data

Hello :)

I am facing problems determing velocity and position from my acceleration data. Therefor I am using cumtrapz, but the results are not what I expect them to be.

Here is what I am doing: I have recorded acceleration data for 3 minutes (without moving the sensor, so the data is just showing the sensor noise). I am using a matlab script, to remove the offset of the recorded data (using mean function). The next step is to determine velocity and position. When I look at the velocity, the course of the function goes back to zero at the end (see first picture). The position resulting from this velocity is about -20 m, which is quite big (I need to implement filtering next).

The problem I am facing with cumtrapz is the following: When I use only the first 120 seconds of the given signal (180 seconds represent the full signal), the signal ends at zero again (see picture 2). The value should be the same as before (looking at the full signal, 180 seconds, upper part in picture 2).

I did a quick check using a sine wave (see picture 3). When I integrate the sine wave using cumtrapz, everything is the way I expect it to be (full signal = 20 seconds, cut off @10 seconds): When comparing the full signal with only 10 seconds of the full signal, the integrated value is the same @10 seconds.

Does anybody have a suggestion what I am doing wrong here?

And here is the matlab script I am using (sorry, forgot that...):

%%init
clear;                          % clear workspace
close all;                      % close opened figures

%%defines
samplefreq_kx  = 0.00125;       % Accel data has been recorded with 1600 hz -> 1600Hz/2=800hz

%%Step 1: Load profile
addpath('C:...');  % Define path with profiles

% Define profile to load here:
load('Profil_20_2.mat');          % Profil 20_2 = Sensor KX022 in hold position, 3 minutes.

t_180_kx = 0:samplefreq_kx:180;   % define time array for 180 seconds
t_120_kx = 0:samplefreq_kx:120;   % define time array for 120 seconds

% Use acceleration data of KX022 @180 seconds 
% Signal has to be converted, because sensor gives out g!
% Use only every second entry in KX022_Signals
% Entry 8001 = Start at 5 seconds (first 5 seconds of signal is not
% important!)
accel_z_kx_180  = KX022_Signals(4,8001:2:end) * 9.80665;    % convert from g to m/s2

% Use acceleration data of KX022 @120 seconds
accel_z_kx_120  = KX022_Signals(4,8001:2:200001) * 9.80665;

%%Remove offset of given signals (for 180 & 120 seconds)
% mean calculates the given offset
% acceleration value has to be zero without offset!
accel_z_kx_180_offsetfree  = accel_z_kx_180  - mean(accel_z_kx_180);
accel_z_kx_120_offsetfree  = accel_z_kx_120  - mean(accel_z_kx_120);

%%Plot offset free signals (180s & 120s)

figure('Name','Acceleration data','NumberTitle','off')
ax1 = subplot(2,1,1);
plot(t_180_kx,accel_z_kx_180_offsetfree)
grid on;
title('KX022 z acceleration (@180s)');
xlabel('time [s]');
ylabel('accel [m/s^2]');

ax2 = subplot(2,1,2);
plot(t_120_kx,accel_z_kx_120_offsetfree)
grid on;
title('KX022 z acceleration (@120s)');
xlabel('time [s]');
ylabel('accel [m/s^2]');

% link y-axis of subplot ax1 and ax2
% First argument in linkaxes defines which subplot is used as reference
% (ax2 has bigger noise)
linkaxes([ax2,ax1],'y');

%%Calculate velocity for 180 & 120 seconds
v_z_kx_180  = cumtrapz(t_180_kx,accel_z_kx_180_offsetfree);
v_z_kx_120  = cumtrapz(t_120_kx,accel_z_kx_120_offsetfree);

figure('Name','Calculated velocity','NumberTitle','off')
ax3 = subplot(2,1,1);
plot(t_180_kx,v_z_kx_180)
grid on;
title('KX022 calculated velocity @180s');
xlabel('time [s]');
ylabel('velocity [m/s]');

ax4 = subplot(2,1,2);
plot(t_120_kx,v_z_kx_120)
grid on;
title('KX022 calculated velocity @120s');
xlabel('time [s]');
ylabel('velocity [m/s]');

linkaxes([ax3,ax4],'y');
linkaxes([ax3,ax4],'x');

%%Calculate position for 180 & 120 seconds
s_z_kx_180  = cumtrapz(t_180_kx,v_z_kx_180);
s_z_kx_120  = cumtrapz(t_120_kx,v_z_kx_120);

figure('Name','Calculated position','NumberTitle','off')
ax5 = subplot(2,1,1);
plot(t_180_kx,s_z_kx_180)
grid on;
title('KX022 calculated position @180s');
xlabel('time [s]');
ylabel('position [m]');

ax6 = subplot(2,1,2);
plot(t_120_kx,s_z_kx_120)
grid on;
title('KX022 calculated position @120s');
xlabel('time [s]');
ylabel('position [m]');

linkaxes([ax5,ax6],'y');
linkaxes([ax5,ax6],'x');

Best regards,

Niko

[John D'Errico]

Since the one thing you have not provided is your code, it is difficult to guess. It would be much easier if you simply attached your script as an m-file to a comment, or best is if you can paste it into your question. Please use the "{} Code" button to format the code so it is readable though.

[Student88]

sorry, I totally forgot that :) I will update my question right now ;)

[John D'Errico]

I'll be darned! COMMENTS! :)

Answer 1

[John D'Errico]

Ok, so you have a long sequence of acceleration measurements, sampled at uniform times.

You subtract off the mean, then use cumtrapz to integrate them, yielding velocity. The starting velocity is zero, as one would expect, but you are surprised the end velocity is also zero. (It is good to check what you see. Plot everything! Test what you see against intuition. If it makes sense, then you are probably doing things right. So all of this is good to see.)

Lets look at the reason why the endpoint velocity is coming out zero. I don't have your data, so I can't run the code on your example, but I should have an explanation for you.

When you subtract off the mean of the series, you are performing a numerical integration, of sorts. cumtrapz is one form of numerical integration, that applies a cumulative trapezoidal rule. In fact, that last value from cumtrapz should be exactly what trapz would yield on the same series. Cumtrapz just gives you the intermediate results too.

But look at what you did to the series before you passed it into cumtrapz. You subtracted off the mean of the series. In fact, if you take a series like that, sampled at uniform times, then compute the mean, you are applying rectangle rule to the series. In effect, that mean subtraction is adjusting the series so that the integral will be zero! It does so only approximately, because rectangle rule is not quite as good an approximation as trapezoidal rule. ANY such numerical integration is only an approximation of course. But on a equally sampled series, especially a moderately long one, it turns out that rectangle rule and trapezoidal rule will be very nearly identical in their prediction.

So by mean-subtracting the series before you throw cumtrapz at it, you are essentially enforcing the result that the endpoint integral will be zero. Again, I predict this will not be exact. So if you look carefully at that final value, it should be only approximately zero. Close, but not exactly zero. The difference will be the difference in the two numerical integration estimates for that series over that domain.

This explains why when you changed the length of the series, the endpoint integral was always approximately zero. That was because the mean subtraction was also over the same domain.

[Student88]

That's funny :D You where just a few minutes faster than me :)

Thanks for the reply ;)

Is there a more accurate way to determine the offset than via mean?

[John D'Errico]

I'll make up some random data. Not very interesting data, I'll admit. But random.

x = rand(1,100);

First, see that trapz and the last element from cumtrapz are identical.

xc = cumtrapz(x);
xc(end)
ans =
       49.996

trapz(x)
ans =
       49.996

As well, sum here is also a variant of rectangle rule. Since the step size is 1, it IS a rectangle rule.

sum(x)
ans =
       50.488

Mean is the same thing as sum of course, since the mean of a vector is just the sum divided by the number of terms.

So, what happens when you subtract off the mean, and then apply trapz, or cumtrapz? If your goal is to adjust the function the series represents, so that the FUNCTION has average value of zero over that interval, then it is best to use the same integration method to compute the average value of the function over that interval.

So what is the support of this function? It lives on the interval [1,n], here n=100. So the length of the interval is n-1.

avefun = trapz(x)/(n-1);

Now, lets test the end point result.

trapz(x - avefun)
ans =
     0

xc = cumtrapz(x - avefun);
xc(end)
ans =
   2.2204e-16

Zero (or effectively so) in both cases. Don't get misled by that tiny trash amount there. It is zero, except that some trash in the least significant digits got in the way. Still zero.

Had I used the mean though, I'd have gotten a result that was not quite as good, because the mean uses an implicit rectangle rule integration instead of trapz.

trapz(x - mean(x))
ans =
     0.013365

All of these differences are small. In fact, it probably does not matter a whole lot which adjustment scheme you used there, since for a sufficiently long vector, the difference will be relatively a small number. Compared to the junk in the data, it might be insignificant. More important is to understand where the difference arises and why.

Yes, if this vector were accelerations for a spaceship on the way to Mars, then it might matter, because the life of the astronauts in your care might be in the balance. A big part of numerical computing is knowing when to worry about the approximation you will make, because you are constantly making one approximation after another. Having a good feel for the important ones can be crucial, and a valuable skill to hone.

[Student88]

Thanks a lot :) That was really helpful :)

Yes, the difference between those two ways is small, but it might be a way to get better results when calculating position from the given acceleration :)

Answer 2

[Student88]

I found the error :)

After verifying the process of cumtrapz with excel (I got exactly the same results in excel), I found the error accidentally.

%%Remove offset of given signals (for 180 & 120 seconds)
% mean calculates the given offset
% acceleration value has to be zero without offset!
accel_z_kx_180_offsetfree  = accel_z_kx_180  - mean(accel_z_kx_180);
accel_z_kx_120_offsetfree  = accel_z_kx_120  - mean(accel_z_kx_120);

When you calculate the offset via mean function, there are (of course...) two different offsets for the full signal (180 seconds) and for the shortened signal (120 seconds). This is, why the plots for 180 & 120 seconds look different. The attached image shows the correct plot, with the red line as the shortened signal :)

Best regards,

Niko

[John D'Errico]

But again, that all happens because the mean of a series can be interpreted as application of rectangle rule, so just another integration rule like trapz.