Clear Filters
Clear Filters

how to remove drift from noisy acceleration data

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

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 20 Nov 2023
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

Products

Community Treasure Hunt

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

Start Hunting!