Issue with magnitude of cumtrapz and integration of second order system acceleration data

I am working on a school project in which we analyze a second order system to create a magnitude function of a motor driven spring-mass-dampener system.

I have gotten the code to the point where I like the general shape of the graphs it gives, but some of the magntudes are far off (ie ~4000 meters) by the time it gets to the displacement graph. I'm pretty sure at this point I need to run this through a filter. I've been on forums about filters for the last hour and I have no clue what I'm doing there. I attached my lab data if it makes it easier.

This is my current code:

clearvars

ColVstart = 1;

ColXstart = 2;

%Initializing Global Variables

%Data sheet must have Time in collumn A and Acceleration in collumn C

for h = 1:17

%Allow data to be read on each trial sheet

t=xlsread('310FinalLabDataWorking.xlsx',h,'A2:A5000');

a=xlsread('310FinalLabDataWorking.xlsx',h,'C2:C5000');

%Read Data for h trial

sizeV=size(a);

start = sizeV(1)-1000;

a = a(start:sizeV);

t = t(1:1001);

fa = fit(t,a,'fourier3');

as = fa(t);

v = cumtrapz(as);

%Integration of data for velocity

sizev = size(v);

bestfitv=detrend(v,1);

fv = fit(t,bestfitv,'fourier3');

vs = fv(t);

x = cumtrapz(vs);

bestfitx=detrend(x,1);

speed = 15+5*h;

figure(h);

subplot(1,3,1);

plot(t,a);

title('acceleration v time');

xlabel('time');

ylabel('acceleration');

subplot(1,3,2);

plot(t,bestfitv)

title('velocity v time');

xlabel('time');

ylabel('velocity');

subplot(1,3,3);

plot(t,bestfitx)

title('displacement v time');

xlabel('time');

ylabel('displacement');

% if h < 14 %Puts data on a seperate sheet columns A:Z

% VCol = char([ColVstart+64,50]);

% XCol = char([ColVstart+65,50]);

% VcollumnName = string(sprintf('V%d',speed));

% XcollumnName = string(sprintf('X%d',speed));

% VNam = char([ColVstart+64,49]);

% XNam = char([ColVstart+65,49]);

% xlswrite('310FinalLabDataWorking.xlsx',{VcollumnName},19,VNam);

% xlswrite('310FinalLabDataWorking.xlsx',{XcollumnName},19,XNam);

% xlswrite('310FinalLabDataWorking.xlsx',bestfitv,19,VCol);

% xlswrite('310FinalLabDataWorking.xlsx',bestfitx,19,XCol);

% ColVstart = ColVstart+2;

% ColXstart = ColXstart+2;

% else %Puts data on a sperate sheet columns AA:end

% VCol = char([65,ColVstart+38,50]);

% XCol = char([65,ColVstart+39,50]);

% VcollumnName = string(sprintf('V%d',speed));

% XcollumnName = string(sprintf('X%d',speed));

% VNam = char([65,ColVstart+38,49]);

% XNam = char([65,ColVstart+39,49]);

% xlswrite('310FinalLabDataWorking.xlsx',{VcollumnName},19,VNam);

% xlswrite('310FinalLabDataWorking.xlsx',{XcollumnName},19,XNam);

% xlswrite('310FinalLabDataWorking.xlsx',bestfitv,19,VCol);

% xlswrite('310FinalLabDataWorking.xlsx',bestfitx,19,XCol);

% ColVstart = ColVstart+2;

% ColXstart = ColXstart+2;

% end

end

Any help creating a filter or fix for this problem would be much appreciated.

Ps. sorry if my code is inefficient, I'm trying here.

Answer 1

[Mathieu NOE]

hello

accelerometer's signals double integration can be a bit tricky. We must compromise between being as accurate as possible while in the same time avoid integrating constants (can be sensor offset) or drift (poor quality sensor).

So it's a matter of using detrend and sometimes even add a bit of high pass filtering - as in the code below. The filter cut off frequency must be chhosen to not be "too high" , means , removing valuable info from the data. So a good thing to do is to perform a fft of the accel signal and see below which frequency the energy start to drop enough so we can put our filter frequency at this value without impacting too much the accuracy of the double integration.

the data file of this demo code is attached .

Hope it helps

data = readmatrix('GEE 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]);

[Jared Lee Pratt]

Thank you this helped a lot.

[Mathieu NOE]

My pleasure !