how to remove drift from noisy acceleration data

Hello, I need your help.

I used an imu sensor to get acceleration data and integral acceleration data to find the velocity

But there's a drift at the end

(It doesn't move after 2 seconds, but it's not a zero value in the graph.)

So I tried offset adjustment, low-pass filter, detrend, etc., but I couldn't find a solution. I want to reduce the drift at the end, what should I do?

I will attach the code I used and the graph I executed

Please,, help me...Any help is appreciated

%% Data preprocessing

% --- left ---
T_left = readtable('ALT1.csv', 'VariableNamingRule', 'preserve');
timestamp_left = T_left.uni_time(2:end);
time_left = timestamp_left

time_left = 492×1

65713070 65713080 65713090 65713100 65713236 65713240 65713250 65713260 65713270 65713280

time_left = timestamp_left-timestamp_left(1)

time_left = 492×1

0 10 20 30 166 170 180 190 200 210

yacc_left = T_left.("watch_acc.y")(2:end);

% --- right ---
T_right = readtable('ART1.csv', 'VariableNamingRule', 'preserve');
timestamp_right = T_right.uni_time(2:end);
time_right = timestamp_right

time_right = 485×1

65713236 65713246 65713256 65713266 65713276 65713286 65713296 65713306 65713316 65713326

time_right = timestamp_right-timestamp_right(1)

time_right = 485×1

0 10 20 30 40 50 60 70 80 90

yacc_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_y_left = mean(yacc_left(stationary_indices));
yacc_left_corrected = yacc_left - offset_y_left;

% Repeat offset correction for the right hand
offset_y_right = mean(yacc_right(stationary_indices));
yacc_right_corrected = yacc_right - offset_y_right;

%% Velocity integration for vector magnitude %!!! This is the integral equation I used !!!
velocity_left_mag = getvelocity(yacc_left_corrected,time_left);
velocity_right_mag = getvelocity(yacc_right_corrected,time_right);



%use low-pass
fs = 1000

fs = 1000

cutoff_freq = 5; 
lp_velocity_left = lowpass(velocity_left_mag, cutoff_freq, fs);
lp_velocity_right = lowpass(velocity_right_mag, cutoff_freq, fs);

time_left = (timestamp_left - timestamp_left(1)) / 1000; % change to second(s) from unix_tinestanp
time_right = (timestamp_right - timestamp_right(1)) / 1000; 



%% 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

Answer 1

[Sulaymon Eshkabilov]

There are a couple of critical issues in your code. (1) Sampling frequency is 100 Hz and NOT 1000 Hz. and (2) Low-pass filter should be applied to acceleration data and NOT velocity data. And some detrending migh be necessary for velocity.

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

% --- right ---
T_right = readtable('ART1.csv', 'VariableNamingRule', 'preserve');
timestamp_right = T_right.uni_time(2:end);
time_right = timestamp_right;
time_right = timestamp_right-timestamp_right(1);
yacc_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_y_left = mean(yacc_left(stationary_indices));
yacc_left_corrected = yacc_left - offset_y_left;

% Repeat offset correction for the right hand
offset_y_right = mean(yacc_right(stationary_indices));
yacc_right_corrected = yacc_right - offset_y_right;
% Acceleration is filtered:
fs = 100;  % See 1/(time_left(ii+1)-time_left(ii)) and 1/(time_right(ii+1)-time_right(ii))
cutoff_freq = 5; 
yacc_left_corrected = lowpass(yacc_left_corrected, cutoff_freq, fs);
yacc_right_corrected = lowpass(yacc_right_corrected, cutoff_freq, fs);
% Velocity integration for vector magnitude %!!! This is the integral equation I used !!!
velocity_left_mag = getvelocity(yacc_left_corrected,time_left);
velocity_right_mag = getvelocity(yacc_right_corrected,time_right);
velocity_L1 = detrend(velocity_left_mag,'linear');
velocity_R1 = detrend(velocity_right_mag,'linear');
time_left = (timestamp_left - timestamp_left(1)) / 1000; % change to second(s) from unix_tinestanp
time_right = (timestamp_right - timestamp_right(1)) / 1000; 
figure
plot(time_left, velocity_L1, 'r', 'LineWidth', 2); hold on; grid on
plot(time_right, velocity_R1, 'b', 'LineWidth', 2); 
legend('v_{left} from getvelocity', 'v_{right} from getvelocity',...
       'Location', 'best');
xlim([0, 5])

% 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