Drift in cumtrapz forceplate data

Dear comunity,

For a project I am analysing gait paterns with a force plate. Now I am using the function Cumtrapz to integrate the signal of the forceplate to get the velocity of the center of mass (COM) in the vertical direction of the person. The Fbw is the body weight what is measured by the mean of the signal of the forceplate when the person is standing still. Fres is the resulting force when to signal is compensated for the bodyweight. After that I calculate the acceleration of the COM, velocity and position of the COM. dt is 1/samplefrequency. Is there a posibility to adjust this script to the drift of the signal?

----

Fbw = mean(FPdata(1:100));

Fres = FPdata - Fbw;

ACOM = Fres/(Fbw/g);

VCOM = cumtrapz(ACOM)*dt;

PCOM = cumtrapz(VCOM)*dt;

-----

Answer 1

[Mathieu NOE]

hello

yet another example of accelerometer / signal integration drift using cumtrapz.

I see you removed the mean value of the signal wich is quite similar to what detrend does , but detrend can also remove a constant drift (optional) which is quite interesting for baseline correction;

now if you are only interested in the dynamic content of the signal and the DC constants are not the key element, you need to use more detrend and maybe add some additionnal high pass filtering too

see example below : (data file attached)

data = readmatrix('Earthquake.xlsx');
t = data(:,1);
dt = mean(diff(t));       % Average dt
fs = 1/dt;                  % Frequency [Hz] or sampling rate

% Acceleration data
acc = data(:,2)*9.81;         % Add gravity factor ( g to m/s²)
acc = detrend(acc);     % baseline correction

% some additionnal high pass filtering     % baseline correction
N = 2;
fc = 0.25; % Hz
[B,A] = butter(N,2*fc/fs,'high');
acc = filter(B,A,acc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
velocity = cumtrapz(dt,acc);
velocity = detrend(velocity);    % baseline correction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp     = cumtrapz(dt,velocity);
disp = detrend(disp);    % baseline correction

%
figure(1)
subplot(311),plot(t,acc);
title('accel'); ylabel('m/s²');xlabel('time (s)');
subplot(312),plot(t,velocity);
title('velocity'); ylabel('m/s');xlabel('time (s)');
subplot(313),plot(t,disp);
title('disp'); ylabel('m');xlabel('time (s)');


%%fft
nfft = length(acc);
fft_spectrum = abs(fft(disp,nfft));

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*fs/nfft;

figure(2)
plot(freq_vector,fft_spectrum);
title('disp'); ylabel('m');xlabel('Frequency (Hz)');

xlim([0 10]);

[Twan Zwetsloot]

It was the high pass filter I needded, thanks!

[Mathieu NOE]

My pleasure !