function [subject]=preprocessing_individual(f, StopCrit,forcetype,subject)
%%[subject]=preprocessing_individual(f,trialtype, StopCriteria,nsub,forcetype,subject,reg,p)
% input arguments:
% f:           the force factor leading to the force perturbing the
%              arm movement by f*velocity; to caluclate the predicted force.
% StopCrit:    either end of state moving (moving) or end of stop moving state (stopmoving);
% forcetype:   'Handle'/'robot' which channeldata to use 'Handle' or 'robot' measured forces;

nTrials = subject.data.Trials;                                              % amount of trials
%%
[b,a]=butter(5,12/500);                                                     % filter parameters                                                                                                 
%% predefined variables to allocate memory 
subject.posx = NaN(nTrials,2000);
subject.posy = NaN(nTrials,2000);
subject.misses = NaN(nTrials,1);

%% all possible targets:
[targets,i] = unique(subject.data.TargetPosition(:,1)); 
targets(:,2) = subject.data.TargetPosition(i,2);
for trial = 1:nTrials
    disp(trial)
    frames = 1:subject.data.FrameData.Frames(trial);                       % index of frames of the current trial
    px = subject.data.FrameData.RobotPosition(trial,frames,1);             % x position data of current trial
    py = subject.data.FrameData.RobotPosition(trial,frames,2);             % y position data of current trial
    pt = subject.data.TargetPosition(trial,1:2);                           % target location
    time = subject.data.FrameData.TrialTime(trial,frames);                 % time vector
    State = subject.data.FrameData.State(trial,frames);                    % States of this trial and frames
    %%
    if strcmp(forcetype,'Handle')                                                                                                               
        ChannelForces = subject.data.FrameData.ChannelHandleForce(trial,frames);% Forces measured at the handle
    elseif strcmp(forcetype,'Robot')                                                                                                            
        ChannelForces = subject.data.FrameData.ChannelForces(trial,frames);% Forces measured at the robot
    end
    
    %%  we first interpolate to true 1000 Hz 
    itime = 0:0.001:max(time);                                             % time to which we interpolate
    ipx = interp1(time,px,itime,'spline','extrap');                        % interpolated x-position
    ipy = interp1(time,py,itime,'spline','extrap');                        % interpolated y-position
    iF = interp1(time,ChannelForces,itime,'linear','extrap');              % interpolated Forces
    iState = interp1(time,State,itime,'nearest','extrap');                 % interpolated States
    %% and filtered the interpolated position data 
    ifpx = filtfilt(b,a,ipx);
    ifpy = filtfilt(b,a,ipy);
    %% velocity and speed
    vx = gradient(ifpx,1e-3);                                               % calculate x velocity from filtered and interpolated position data
    vy = gradient(ifpy,1e-3);                                               % calculate y velocity from filtered and interpolated position data
    speed= sqrt(vx.^2+vy.^2);                                               % calculate speed
    distance_to_target = sqrt((ifpx-pt(1)).^2+(ifpy-pt(2)).^2);             % distance to the target location
    ds_t = gradient(distance_to_target,1e-3);                               % ensuring the participant is moving towards the target location
    distance_to_start = sqrt((ifpx).^2+(ifpy).^2);                          % distance to the target location
    ds_s = gradient(distance_to_start,1e-3);                                % ensuring the participant is moving away from the start location
    speed_for_sample_selection = speed;
    speed_for_sample_selection(sign(-ds_t)>0|sign(ds_s>0)) = speed(sign(-ds_t)>0|sign(ds_s>0));
    speed_for_sample_selection(sign(-ds_t)<=0&sign(ds_s<=0)) = -speed(sign(-ds_t)<=0&sign(ds_s<=0));
    
    %% relevant samples
    if strcmp(StopCrit,'Moving')
        raw_samples = find(iState == 8);
        % If there are blips before the start of the movement, take the next frame
        % after the last blip as the start of index of movement
        if sum(diff(find(iState == 8))>1)>0
            sample_starts = find(iState == 8);
            raw_samples = sample_starts(find(diff(sample_starts)>1,1,'last')+1):raw_samples(end);% if so take the next frame after the last blip as the start of index of movement
        end
    elseif strcmp(StopCrit,'StopMoving')
        raw_samples = find(iState == 8| iState == 9);
        % If there are blips before the start of the movement, take the next frame
        % after the last blip as the start of index of movement
        if sum(diff(find(iState == 8))>1)>0
            sample_starts = find(iState == 8);
            raw_samples = sample_starts(find(diff(sample_starts)>1,1,'last')+1):raw_samples(end);% if so take the next frame after the last blip as the start of index of movement
        end
    end 
    %%  make sure we take the bellshape around peak velocity
    [v,i] = max(speed_for_sample_selection(raw_samples));
    [V,I] = findpeaks(speed_for_sample_selection(raw_samples),'MinPeakHeight',0.6*v); 
    %% find peaks of at least 60% of the maximum peak
     if i~=I(1)
        disp('check this trial')
        i = I(1);
        v =V(1);
        subject.misses(trial) = 4;
     end
    %% finding the start and end of the movement that is larger than 15% of maximum speed
    start = find(speed(raw_samples(1:i))<0.15*v,1,'last');
    finish_speed = find(speed(raw_samples(i:end))<0.15*v,1,'first');
    finish_sign = find(gradient(speed(raw_samples(i:end)))>0 & speed(raw_samples(i:end))<(0.5*v),1,'first');
    finish = min([finish_speed,finish_sign]);
    samples = raw_samples(start:(i+finish-1));
   %% write sampled data to subject.r and subject.posx/posy
    if  isempty(samples)
        disp('no Samples')
        %pause
        AI = NaN;
        subject.r(trial) =AI;
        subject.fpredicted{trial} =NaN;
        subject.f{trial} = NaN;
        subject.misses(trial) = NaN;
        n = length(samples);
        subject.posx(trial,1:n) = NaN;
        subject.posy(trial,1:n) = NaN;
        subject.misses(trial) = 1;
    else
        
        AI = mean(iF(samples)./(f*speed(samples)));
        subject.r(trial) =AI;
        subject.fpredicted{trial} =f*speed(samples);
        subject.f{trial} = iF(samples);
        subject.speed{trial} = speed(samples);
        n = length(samples);
        subject.posx(trial,1:n) = ipx(samples);
        subject.posy(trial,1:n) = ipy(samples);
    end
   

end

