how to remove drift from noisy acceleration data

7 ビュー (過去 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 Center および File ExchangeSimulation, Tuning, and Visualization についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by