Drift removal of velocity obtained by integrating acceleration data

Hi, can you help me?

I'm taking acceleration data from the IMU sensor and getting a velocity graph through integration and then looking at the velocity of motion by vector sum.

By the way, even though the drift has occurred and the movement is over, it remains at the speed of a point in time, not zero, is there any way to solve this?

I adjusted the offset because there was no movement before 0.47 seconds, and the blue line is the same, but based on the red line, there is no movement after 2.176 seconds (movement ends almost simultaneously with left and right)..

I used "detrend" to correct the trend, but as I did that, the graph tilted and transformed into a motionless section

Do you know how to solve this problem? Any help is appreciated

%% Data preprocessing

% --- left hand ---
T_left = readtable('ALT10.csv', 'VariableNamingRule', 'preserve');
timestamp_left = T_left.uni_time(2:end);
time_left = timestamp_left
time_left = timestamp_left-timestamp_left(1)
xacc_left = T_left.("watch_acc.y")(2:end);
yacc_left = T_left.("watch_acc.y")(2:end);
zacc_left = T_left.("watch_acc.y")(2:end);

% --- right hand ---
T_right = readtable('ART10.csv', 'VariableNamingRule', 'preserve');
timestamp_right = T_right.uni_time(2:end);
time_right = timestamp_right
time_right = timestamp_right-timestamp_right(1)
xacc_right = T_right.("watch_acc.y")(2:end);
yacc_right = T_right.("watch_acc.y")(2:end);
zacc_right = T_right.("watch_acc.y")(2:end);


%% Offset correction and Filtering

stationary_indices = 1:30; % Indices of samples where the device is assumed stationary

% Offset correction for all axes
offset_x_left = mean(xacc_left(stationary_indices));
xacc_left_corrected = xacc_left - offset_x_left;
offset_y_left = mean(yacc_left(stationary_indices));
yacc_left_corrected = yacc_left - offset_y_left;
offset_z_left = mean(zacc_left(stationary_indices));
zacc_left_corrected = zacc_left - offset_z_left;

% Repeat offset correction for the right hand
offset_x_right = mean(xacc_right(stationary_indices));
xacc_right_corrected = xacc_right - offset_x_right;
offset_y_right = mean(yacc_right(stationary_indices));
yacc_right_corrected = yacc_right - offset_y_right;
offset_z_right = mean(zacc_right(stationary_indices));
zacc_right_corrected = zacc_right - offset_z_right;



%% Velocity integration for vector magnitude
xvelocity_left_mag = getvelocity(xacc_left_corrected,time_left);
yvelocity_left_mag = getvelocity(yacc_left_corrected,time_left);
zvelocity_left_mag = getvelocity(zacc_left_corrected,time_left);
xvelocity_right_mag = getvelocity(xacc_right_corrected,time_right);
yvelocity_right_mag = getvelocity(yacc_right_corrected,time_right);
zvelocity_right_mag = getvelocity(zacc_right_corrected,time_right);

%for you information, getvelocity is this..
%function v = getvelocity(acc, time)
  %  v = zeros(1, length(acc));
    
   % for i = 1:length(acc)-1
    %    dt = time(i+1) - time(i);
        % numerical integration
     %   v(i+1) = v(i)+0.5*dt*(acc(i) + acc(i+1));
  %  end
%end


fs = 1000

cutoff_freq = 5; % 보다 낮은 컷오프 주파수 설정
xlp_velocity_left = lowpass(xvelocity_left_mag, cutoff_freq, fs);
ylp_velocity_left = lowpass(yvelocity_left_mag, cutoff_freq, fs);
zlp_velocity_left = lowpass(zvelocity_left_mag, cutoff_freq, fs);

xvelocityR= lowpass(xvelocity_right_mag, cutoff_freq, fs);
yvelocityR = lowpass(yvelocity_right_mag, cutoff_freq, fs);
zvelocityR = lowpass(zvelocity_right_mag, cutoff_freq, fs);

vectorsumR = sqrt(xvelocityR.^2 + yvelocityR.^2 + zvelocityR.^2);
vectorsumL = sqrt(xlp_velocity_left.^2 + ylp_velocity_left.^2 + zlp_velocity_left.^2)

time_left = (timestamp_left - timestamp_left(1)) / 1000; % change to second(S) frome unix_time
time_right = (timestamp_right - timestamp_right(1)) / 1000;

Answer 1

[Walter Roberson]

It turns out that even the best single-mode IMU systems encounter substantial drift

https://en.wikipedia.org/wiki/Inertial_navigation_system#Drift_rate

"Even the best accelerometers, with a standard error of 10 micro-g, would accumulate a 50-meter (164-ft) error within 17 minutes"

Your data only has about 3 significant figures, 4 at best, and your labels talk about mm; I figure that your error is probably on the order of 10 milli-g rather than 10 micro-g, so I estimate that you would probably be about 50 m off in... slightly over 1 second. If your actual error is 1 milli-g then 50 m in 10.2 seconds I estimate, and over your 2.176 seconds... you should still expect over 10 meters positioning error.