how to remove drift from noisy acceleration data

조회 수: 6 (최근 30일)
JUNY
JUNY 2023년 11월 20일
답변: Sulaymon Eshkabilov 2023년 11월 20일
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

답변 (1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023년 11월 20일
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

카테고리

Help CenterFile Exchange에서 Simulation, Tuning, and Visualization에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by