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
