Velocity integrating to get displacement issues

Hello, I am recording the accelerometer data and I want to get an displacement from acceleration. This means, that I have to integrate acceleration twice. When I integrate acceleration first time I get velocity and I am satisfied with the result (see attached). But strange thing happens when I integrate velocity to get a displacement. The graph is just not right. It should be pretty much the same as velocity. I guess the cumtrapz should integrate velocity the same way as it integrates acecleration. Do you have any idea what am I doing wrong? Thanks in advance

taxis=xlsread('test.xls', 'A:A');
yaxis=xlsread('test.xls','C:C');

length=30;
mask = taxis > length;
taxis(mask) = [];
yaxis(mask) = [];

figure(1)
plot(taxis, yaxis)
xlabel('Time, s')
ylabel('Acceleration in Y axis, m/s^2')
grid on

vel=cumtrapz(taxis,yaxis);
detrend_vel=detrend(vel);

figure(2)
plot(taxis, vel)
hold on
[pks,locs] = findpeaks(detrend_vel, 'minpeakdistance',5, 'minpeakheight', 0.2);
plot(taxis, detrend_vel, 'r',taxis(locs), pks, 'or')
grid on
xlabel('Time, s')
ylabel('Velocity, m/s')
legend('Raw','Detrended')

displacement=cumtrapz(taxis,vel);
detrend_displacement=detrend(displacement);
figure(3)
plot(taxis, displacement)
hold on
plot(taxis, detrend_displacement)
grid on
xlabel('Time, s')
ylabel('Displacement, m')
legend('Raw','Detrended')

Answer 1

[Ameer Hamza]

Try to cumtrapz the detrended velocity

displacement=cumtrapz(taxis,detrend_vel);
detrend_displacement=detrend(displacement);

[Karolis Naujalis]

Hello Ameer, thanks for the quick reply.

I have already tried it - it is not working as well. The real displacement values are ~5-6 cm (its chest compressions during CPR)

[Ameer Hamza]

In that case, you will need to use more complex filtering techniques like Kalman filter. Directly integrating the acceleration value twice to get displacement is a well-known problem. It never works. Unless you have a really accurate accelerometer. Usually, a model-based filtering algorithm, e.g., Kalman filter, is used to get accurate results.

Also, you mentioned that the correct displacement is in the range of 5-6 cm. do you think that the small ripples are the correct signal? And the big sine wave component is the sensor issue? Although I doubt that chest compression can be done at such fast speed.

[Karolis Naujalis]

Ameer,

I used the same sensor and the same code before and it worked for me, but now it doesn't.

Chest compressions are 100-120 times per minute (~ 2 times in one second)

Do you have an idea, why it doesn't work any more for me? Previously I got similar displacement plot as with velocity

[Ameer Hamza]

I guess the small ripples are correct signals, and somehow a lower-frequency component (of frequency 30 seconds) is added due to some external factor. Was the platform itself moving when you recorded the data? Can you attach a small sample of the raw data?

[Karolis Naujalis]

I use the same data as before, see attached.

Platform was not moving, because the CPR manikin was laying flat on the floor.

I'm feeling really stuck now with this

EDIT:

What if I decide to measure the width of peaks (pulse witdth=time) and then I would use simple formula: displacement = peak velocity * pulse width in seconds? that should also give me a displacement

[Ameer Hamza]

As the error signal has a very low frequency, so you can use a high-pass filter to block out that component of your final signal. I set the stop-band frequency to 0.5 Hz, because the required signal is 2 Hz, and the unwanted signal is approximately 1/30 Hz. Try this code

taxis=xlsread('test.xlsx', 'A:A');
yaxis=xlsread('test.xlsx', 'C:C');

length=30;
mask = taxis > length;
taxis(mask) = [];
yaxis(mask) = [];

figure(1)
plot(taxis, yaxis)
xlabel('Time, s')
ylabel('Acceleration in Y axis, m/s^2')
grid on

vel=cumtrapz(taxis,yaxis);
detrend_vel=detrend(vel);

figure(2)
plot(taxis, vel)
hold on
[pks,locs] = findpeaks(detrend_vel, 'minpeakdistance',5, 'minpeakheight', 0.2);
plot(taxis, detrend_vel, 'r',taxis(locs), pks, 'or')
grid on
xlabel('Time, s')
ylabel('Velocity, m/s')
legend('Raw','Detrended')

displacement=cumtrapz(taxis,detrend_vel);
ts = mean(diff(taxis));
fs = 1./ts;
filtered_displacement = highpass(displacement, 0.5, fs);

figure(3)
plot(taxis, displacement)
hold on
plot(taxis, filtered_displacement)
grid on
xlabel('Time, s')
ylabel('Displacement, m')
legend('Raw','Filtered')

[Karolis Naujalis]

Thanks so much for your help, Ameer, that is what I wanted! The filter was all I need, you are amazing! I hoped, that I will get average displacement ~5 cm, but now it's 2.83 cm. Do think it's because of filter or of signal itself?

[Ameer Hamza]

It appears to be part of the original signal. If you look at the ripples of unfiltered, they also don't have a peak-to-peak distance of 10cm (5cm amplitude). If you are sure that it must be 5 and there might be some issue in processing the data, then you can multiply it with a constant factor to make its amplitude equal to 5.

[Karolis Naujalis]

Oh my, Ameer, you are MVP, huge thanks for your effort!

[Ameer Hamza]

I am glad to be of help.