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

16 views (last 30 days)
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.

Accepted Answer

Mathieu NOE
Mathieu NOE on 27 Apr 2022
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]);

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!