hello
beside detrend, you may have to add some highpass filtering that detrend cannot cope with. Of course the choice of the cut off frequency (fc) is a bit touchy as it depends what the disp signal should look like (zero at the end ? ) . If you know the physical displacement by a separate (calibrated) sensor or you know what is the amount of displacement imposed on your sensor , then you can tweak the filter's characteristic. Too much high pass filtering will reduce the "good" signal output , not enough will leave too much dc drift / low frequency noise that will corrupt the integrals.
try this :
clear;
filename = 'calibration experiment #2 11 dec 22.txt';
warning off;
[acc, dt, T] = ReadAccelerometer(filename); % help function added at the bottom
samples = length(T);
% time = dt:dt:length(T)*dt';
time = (0:samples-1)*dt;
Fs = 1/dt;
% Butterworth HP filter
fc = 0.5;
[B,A] = butter(1,2*fc/Fs,'high');
% Derive velocity & displacements
acc.x = detrend(acc.x*9.81);
acc.x = filter(B,A,acc.x);
vel.x = cumtrapz(dt,acc.x);
vel.x = detrend(vel.x);
vel.x = filter(B,A,vel.x);
dis.x = cumtrapz(dt,vel.x);
acc.y = detrend(acc.y/10*9.81);
acc.y = filter(B,A,acc.y);
vel.y = cumtrapz(dt,acc.y);
vel.y = detrend(vel.y);
vel.y = filter(B,A,vel.y);
dis.y = cumtrapz(dt,vel.y);
acc.z = detrend(acc.z/10*9.81);
acc.z = filter(B,A,acc.z);
vel.z = cumtrapz(dt,acc.z);
vel.z = detrend(vel.z);
vel.z = filter(B,A,vel.z);
dis.z = cumtrapz(dt,vel.z);
% Plot results
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaulttextinterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
figure
subplot(3,3,1)
set(gcf,'color','w');
plot(time(1:end),acc.x(1:end)); hold on;
ylabel('Acceleration x (g)')
grid on; grid minor;
subplot(3,3,2)
plot(time(1:end),acc.y(1:end)); hold on;
ylabel('Acceleration y (g)')
grid on; grid minor;
subplot(3,3,3)
plot(time(1:end),acc.z(1:end)); hold on;
ylabel('Acceleration z (g)')
xlabel('Time (s)');
grid on; grid minor;
subplot(3,3,4)
set(gcf,'color','w');
plot(time(1:end),vel.x)
ylabel('Velocity x (m)')
grid on; grid minor;
subplot(3,3,5)
plot(time(1:end),vel.y)
ylabel('Velocity y (m)')
grid on; grid minor;
subplot(3,3,6)
plot(time(1:end),vel.z)
ylabel('Velocity z (m)')
xlabel('Time (s)');
grid on; grid minor;
subplot(3,3,7)
set(gcf,'color','w');
plot(time(1:end),dis.x)
ylabel('Displacement x (m)')
grid on; grid minor;
subplot(3,3,8)
plot(time(1:end),dis.y)
ylabel('Displacement y (m)')
grid on; grid minor;
subplot(3,3,9)
plot(time(1:end),dis.z)
ylabel('Displacement z (m)')
xlabel('Time (s)');
grid on; grid minor;
function [acc, dt, T] = ReadAccelerometer(filename)
dat = readtable(filename);
dt1 = char(dat.ChipTime(1));
dt2 = char(dat.ChipTime(2));
dt = (str2double(dt2(end-2:end)) - str2double(dt1(end-2:end)))/1000;
fns = fieldnames(dat);
fns = fieldnames(dat);
for i = 1:11
dat.(fns{i+3}) = str2double(strrep(dat.(fns{i+3}), ',', '.'));
end
T = dat.ChipTime;
acc.x = dat.ax_g_;
acc.y = dat.ay_g_;
acc.z = dat.az_g_;
end
