Drift in cumtrapz forceplate data
14 次查看(过去 30 天)
显示 更早的评论
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;
-----
0 个评论
采纳的回答
Mathieu NOE
2022-5-19
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]);
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!