How to program primary submovement duration in matlab (Kinematics analysis for discrete aiming movements)

In a discrete aiming movement the primary sub-movement starts at movement onset till the second zero-crossing in the acceleration profile (i.e., the first negative-to-positive point) after peak velocity.

I am able to program the peak velocity I am wondering if someone can help me based on my progress to identify tprimary movement


%%%%%%velocity
Fs = 200;  % Sampling Frequency
sec = 1/Fs;
N  = 4;  % Order
Fc = 16;  % Cutoff Frequency

n = 100;
data = zeros(n,2);
data(:,1)= randi(996,n,1)+4;
data(:,2)= randi(446,n,1)

[b, a] = butter(N, Fc/(Fs/2));
xres = 0.001;
yres = 0.001; 


xdata = data(1:n,1) * xres;
ydata = data(1:n,2) * yres;

xdata = filtfilt(b, a, xdata);

ydata = filtfilt(b, a, ydata);

xvel = [0; ((xdata(3:end)-xdata(1:end-2))./(2*sec)); 0];

yvel = [0; (ydata(3:end)-ydata(1:end-2))./(2*sec); 0];

res_vel = sqrt(xvel.^2 + yvel.^2);

vel_thres = (5/100)*max(res_vel);

start_time = find(res_vel>vel_thres, 1 );

end_time = find(res_vel<vel_thres, 1, ‘last’ );

Invalid text character. Check for unsupported symbol, invisible character, or pasting of non-ASCII characters.

Error in connector.internal.fevalMatlab

Error in connector.internal.fevalJSON

if end_time<find(max(res_vel))
    
    end_time = length(xdata);
    
end

%%%%%acceleration

xaccel1 = [((xvel1(3:end)-xvel1(1:end-2))./(2*sec))];

yaccel1 = [((yvel1(3:end)-yvel1(1:end-2))./(2*sec))];

xaccel = [xaccel1(20:end-19)];

yaccel = [yaccel1(20:end-19)];

res_accel = sqrt(xaccel.^2 + yaccel.^2);

ave_accel = mean(res_accel);


Movement_duration = end_time - start_time ; % This is the movement time for the whole movement, I want to get the movement time for the primary submovement 

[Walter Roberson]

Please expand on the question ?

[Walter Roberson]

%%%%%%velocity
Fs = 200;  % Sampling Frequency
sec = 1/Fs;
N  = 4;  % Order
Fc = 16;  % Cutoff Frequency

n = 100;
data = zeros(n,2);
data(:,1)= randi(996,n,1)+4;
data(:,2)= randi(7,n,1)

data = 100×2

686 2 91 3 242 2 735 4 454 6 801 5 741 1 902 7 735 2 68 3

[b, a] = butter(N, Fc/(Fs/2));
xres = 0.001;
yres = 0.001; 


xdata = data(1:n,1) * xres;
ydata = data(1:n,2) * yres;

xdata = filtfilt(b, a, xdata);

ydata = filtfilt(b, a, ydata);

xvel = [0; ((xdata(3:end)-xdata(1:end-2))./(2*sec)); 0];

yvel = [0; (ydata(3:end)-ydata(1:end-2))./(2*sec); 0];

res_vel = sqrt(xvel.^2 + yvel.^2);

vel_thres = (5/100)*max(res_vel);

start_time = find(res_vel>vel_thres, 1 );

end_time = find(res_vel<vel_thres, 1, 'last' );

if end_time<find(max(res_vel))
    
    end_time = length(xdata);
    
end

%%%%%acceleration

xaccel1 = [((xvel(3:end)-xvel(1:end-2))./(2*sec))];

yaccel1 = [((yvel(3:end)-yvel(1:end-2))./(2*sec))];

xaccel = [xaccel1(20:end-19)];

yaccel = [yaccel1(20:end-19)];

res_accel = sqrt(xaccel.^2 + yaccel.^2);

ave_accel = mean(res_accel);

plot(xaccel, yaccel)