process_participant.m

Go to the documentation of this file.

% =========================================================================== %
%> @brief preprocess and filter the data of a single participant
%>
%> this function preprocesses the data of a single participant. is is written
%> in specific for the SonEMS study participants. the following steps are
%> executed:
%> - load the VICON data
%> - load the Zurich Move data
%> - load the gait summary data
%> - create ground truth positions from the VICON data, i.e. estimate the
%>   position of the JUMP sensor based on the VICON markers
%> - determine the walking axes
%> - determine the walking speed
%> - accumulate the distance based on the walking speed
%> - align the IMU data with the VICON data
%> - find the yaw of the IMU data
%> - filter the data
%>
%> @param basepath: (string) the basepath of the data, e.g. "D:/ETH_Data/"
%> @param ptp: (string) the participant to process, e.g. "SonE_01"
%> @param options the options for the processing:
%> - GT_POS: (logical) create ground truth positions from the VICON data
%> - FIND_WALKING_AXES: (logical) determine the walking axes
%> - FIND_WALKING_SPEED: (logical) determine the walking speed
%> - ACCUMULATE_DIST: (logical) accumulate the distance based on the walking speed
%> - ALIGN_IMU_VICON: (logical) align the IMU data with the VICON data
%> - FIND_YAW: (logical) find the yaw of the IMU data
%> - FILTER: (logical) filter the data
%> - ekf_tuning_file: (string) the file with the tuning parameters for the EKF,
%>   see @ref EKF1::importTuning()
%>
%> @retval none
%>
%> @note the tasks to do which can be set in the options are used to ensure
%> that a task is redone, even if it has been done before, e.g. if filter
%> parameters were changed. if a task is not done yet, it will be executed
%> regardless of the options.
%>
%> @copyright see the file @ref LICENSE in the root directory of the repository
% =========================================================================== %
function process_participant(basepath, ptp, options)
 
arguments
    basepath string
    ptp string
    options.GT_POS logical = false
    options.FIND_WALKING_AXES logical = false
    options.FIND_WALKING_SPEED logical = false
    options.ACCUMULATE_DIST logical = false
    options.ALIGN_IMU_VICON logical = false
    options.FIND_YAW logical = false
    options.FILTER logical = false
    options.ekf_tuning_file string {mustBeFile}
end
 
% convert options
GT_POS             = options.GT_POS;
FIND_WALKING_AXES  = options.FIND_WALKING_AXES;
FIND_WALKING_SPEED = options.FIND_WALKING_SPEED;
ACCUMULATE_DIST    = options.ACCUMULATE_DIST;
ALIGN_IMU_VICON    = options.ALIGN_IMU_VICON;
FIND_YAW           = options.FIND_YAW;
FILTER             = options.FILTER;
 
% preprocess options to ensure dependencies
if GT_POS
    FIND_WALKING_AXES = true;
    FIND_WALKING_SPEED = true;
    ACCUMULATE_DIST = true;
end
 
if FIND_WALKING_AXES
    FIND_WALKING_SPEED = true;
    ACCUMULATE_DIST = true;
end
 
if FIND_WALKING_SPEED
    ACCUMULATE_DIST = true;
    FILTER = true;
end
 
if ALIGN_IMU_VICON
    FIND_YAW = true;
    FILTER = true;
end
 
if FIND_YAW
    FILTER = true;
end
 
%% path definitions
ptp_folder = fullfile(basepath, ptp);
diary(fullfile(ptp_folder, ptp + ".log"))
ptp_num = regexp(ptp, "SonE_(\d+)|Pilot_(\d+)", "tokens");
ptp_num = str2double(ptp_num{1});
vicon_path = fullfile(ptp_folder, "Processed/Vicon/");
zm_raw_path = fullfile(ptp_folder, "RawData/ZM/*.BIN");
summary_path = fullfile(ptp_folder,"Result/GaitSummaryIMU_SonE_" + ptp_num + ".mat");
temp_path = fullfile(ptp_folder, "temp/");
 
log_entry("Start processing participant", ptp_num)
 
if ~isfolder(temp_path)
    % create the temp folder if it does not exist
    mkdir(temp_path)
    log_entry("Created temp folder", ptp_num, 5)
end
 
% log the tasks that need to be done
log_entry("Tasks to do:", ptp_num, 5)
log_entry("GT_POS: " + string(GT_POS), ptp_num, 5)
log_entry("FIND_WALKING_AXES: " + string(FIND_WALKING_AXES), ptp_num, 5)
log_entry("FIND_WALKING_SPEED: " + string(FIND_WALKING_SPEED), ptp_num, 5)
log_entry("ACCUMULATE_DIST: " + string(ACCUMULATE_DIST), ptp_num, 5)
log_entry("ALIGN_IMU_VICON: " + string(ALIGN_IMU_VICON), ptp_num, 5)
log_entry("FIND_YAW: " + string(FIND_YAW), ptp_num, 5)
log_entry("FILTER: " + string(FILTER), ptp_num, 5)
 
% check if the participants JSON already exists
if isfile(fullfile(ptp_folder, "proc_info.json"))
    % load the JSON file
    proc_info = readstruct(fullfile(ptp_folder, "proc_info.json"));
    log_entry("JSON file loaded", ptp_num, 5)
    
    % set status to not done
    proc_info.filtering_done = false;
    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
else
    % create the JSON file
    proc_info = struct();
    log_entry("JSON file not found, created new struct", ptp_num, 4)
    
    % set status to not done
    proc_info.filtering_done = false;
    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
end
 
%% load files (only if not a Pilot participant)
 
if contains(ptp, "SonE")
    
    % vicon data position data
    if ~isfile(temp_path + "vicon_1.mat")
        
        load(fullfile(vicon_path,"GaitSummary_YK","SonE_" + ptp_num + ".mat"), "GaitSummary");
        vicon_1 = GaitSummary.("SonE_" + ptp_num);
        save(temp_path + "vicon_1.mat", "vicon_1");
        
        log_entry("VICON data loaded", ptp_num, 5)
        
        % remove unnecessary variables (only large ones)
        clear GaitSummary
    end
    
    % zurich move, data ankles
    if ~isfile(temp_path + "zm_l_raw.mat") || ~isfile(temp_path + "zm_r_raw.mat")
        
        zm_raw_files = dir(zm_raw_path);
        for i = 1:numel(zm_raw_files)
            this_header = jumpReadData(fullfile(zm_raw_files(i).folder, zm_raw_files(i).name), header_only=true);
            
            if strcmp(this_header.name, "ankle_L")
                zm_l_raw = jumpReadData(fullfile(zm_raw_files(i).folder, zm_raw_files(i).name));
                log_entry(sprintf("Left ankle data loaded from %s", zm_raw_files(i).name), ptp_num, 5);
            elseif strcmp(this_header.name, "ankle_R")
                zm_r_raw = jumpReadData(fullfile(zm_raw_files(i).folder, zm_raw_files(i).name));
                log_entry(sprintf("Right ankle data loaded from %s", zm_raw_files(i).name), ptp_num, 5);
            end
            
        end
        
        save(temp_path + "zm_l_raw.mat", "zm_l_raw");
        save(temp_path + "zm_r_raw.mat", "zm_r_raw");
        
        log_entry("Zurich Move data loaded", ptp_num, 5)
        
    end
    
    % gait summary
    if ~isfile(temp_path + "sum_zm_l.mat") || ~isfile(temp_path + "sum_zm_r.mat")
        
        load(summary_path, "GaitSummaryIMU");
        sum_zm_l = GaitSummaryIMU.ModulesData.ankle_L;
        sum_zm_r = GaitSummaryIMU.ModulesData.ankle_R;
        
        save(temp_path + "sum_zm_l.mat", "sum_zm_l");
        save(temp_path + "sum_zm_r.mat", "sum_zm_r");
        
        log_entry("Gait summary data loaded", ptp_num, 5)
        
        % remove unnecessary variables (only large ones)
        clear GaitSummaryIMU
    end
    
else
    
    % load the pilot data
    if ~isfile(temp_path + "vicon_1.mat")
        
        
        load(fullfile(vicon_path,sprintf("P%02d.mat",ptp_num)), "GaitSummary");
        vicon_1 = GaitSummary.("P01");
        save(temp_path + "vicon_1.mat", "vicon_1");
        
        log_entry("VICON data loaded", ptp_num, 5)
        
        % remove unnecessary variables (only large ones)
        clear GaitSummary
    end
    
    % zurich move, data ankles
    if ~isfile(temp_path + "zm_l_raw.mat") || ~isfile(temp_path + "zm_r_raw.mat")
        
        zm_l_raw = load(fullfile(ptp_folder, "RawData", "calibrated_ankle_L.mat"));
        zm_r_raw = load(fullfile(ptp_folder, "RawData", "calibrated_ankle_R.mat"));
        
        save(temp_path + "zm_l_raw.mat", "zm_l_raw");
        save(temp_path + "zm_r_raw.mat", "zm_r_raw");
        
        log_entry("Zurich Move data loaded", ptp_num, 5)
        
    end
    
    % gait summary is not needed for pilot participants
    
end
 
%% create VICON ground truth from markers
% there are 4 markers on each foot:
% - heel (Calcaneus - RHEE, LHEE)
% - Toe 1 (Caput metatarsale I - RTO1, LTO1)
% - Toe 3 (Caput metatarsale III - RTO3, LTO3)
% - Toe 5 (Caput metatarsale V - RTO5, LTO5)
%
% the JumpSensor was placed on the metatarsal bones on the proximal side of the
% toe markers. The exact Position was not recorded. Therefore, we will use the
% mean position of the toe markers and transform this position with a fixed
% offset towards the heel marker. An offset of 5cm is used.
%
% note: the vicon data is in mm
 
if ~isfile(fullfile(temp_path,"vicon_2.mat")) || GT_POS
    
    % load the vicon data, if not already done
    if ~exist("vicon_1", "var"); load(temp_path + "vicon_1.mat", "vicon_1"); end
    
    % loop over all trials
    trials = string(fieldnames(vicon_1))';
    for trial = trials
        
        % mean toe left
        vicon_1.(trial).KinematicData.Marker.LTO = (vicon_1.(trial).KinematicData.Marker.LTO1 + vicon_1.(trial).KinematicData.Marker.LTO3 + vicon_1.(trial).KinematicData.Marker.LTO5) / 3;
        % mean toe right
        vicon_1.(trial).KinematicData.Marker.RTO = (vicon_1.(trial).KinematicData.Marker.RTO1 + vicon_1.(trial).KinematicData.Marker.RTO3 + vicon_1.(trial).KinematicData.Marker.RTO5) / 3;
        
        % offset from toe towards heel
        offset = 65; % 6.5 cm
        
        % left foot
        left_vec = vicon_1.(trial).KinematicData.Marker.LHEE - vicon_1.(trial).KinematicData.Marker.LTO;
        left_vec = left_vec ./ vecnorm(left_vec,2,2);
        vicon_1.(trial).KinematicData.Marker.LF = vicon_1.(trial).KinematicData.Marker.LTO + offset * left_vec;
        
        % right foot
        right_vec = vicon_1.(trial).KinematicData.Marker.RHEE - vicon_1.(trial).KinematicData.Marker.RTO;
        right_vec = right_vec ./ vecnorm(right_vec,2,2);
        vicon_1.(trial).KinematicData.Marker.RF = vicon_1.(trial).KinematicData.Marker.RTO + offset * right_vec;
        
        % add the trial name to the struct
        proc_info.trials.(trial).name = trial;
        writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
        
    end
    
    % rename vicon_1 to vicon_2
    vicon_2 = vicon_1;
    
    % save the data
    save(temp_path + "vicon_2.mat", "vicon_2");
    
    log_entry("Ground truth positions created", ptp_num)
    
    % remove unnecessary variables (only large ones)
    clear vicon_1 left_vec right_vec
    
end
 
%% determine the walking axis
if ~isfield(proc_info, "walking_axes") || FIND_WALKING_AXES
    
    if ~FIND_WALKING_AXES
        log_entry("The walking axes have to be determined even though it was not marked (not done yet)", ptp_num, 5)
    end
    
    % load the vicon data, if not already done
    if ~exist("vicon_2", "var"); load(temp_path + "vicon_2.mat", "vicon_2"); end
    
    % loop over all trials
    trials = string(fieldnames(vicon_2))';
    
    sag_axis = 0; % sagital axis
    lon_axis = 0; % longitudinal axis
    lat_axis = 0; % lateral axis
    
    for trial = trials
        
        % get the standard deviation of all three axes
        std_axis_left = std(vicon_2.(trial).KinematicData.Marker.LF);
        std_axis_right = std(vicon_2.(trial).KinematicData.Marker.RF);
        
        % find the walking axis (saggiatal axis) -> highest standard deviation
        [~, idx_left] = max(std_axis_left);
        [~, idx_right] = max(std_axis_right);
        
        % check if the axis is the same for both feet
        if idx_left == idx_right
            if sag_axis > 0 && sag_axis ~= idx_left
                log_entry("The sagital axis of the left and right foot is not constant over all trials", ptp_num, 2, trial)
                
                log_entry(sprintf("Detected Standard Deviations: Left: %.2f %.2f %.2f, Right: %.2f %.2f %.2f", ...
                    std_axis_left(1), std_axis_left(2), std_axis_left(3), std_axis_right(1), std_axis_right(2), std_axis_right(3)), ptp_num, 5, trial)
                
                error("The sagital axis of the left and right foot is not constant over all trials")
            else
                sag_axis = idx_left;
            end
        else
            log_entry("The sagital axis of the left and right foot is not the same", ptp_num, 2, trial)
            
            log_entry(sprintf("Detected Standard Deviations: Left: %.2f %.2f %.2f, Right: %.2f %.2f %.2f", ...
                std_axis_left(1), std_axis_left(2), std_axis_left(3), std_axis_right(1), std_axis_right(2), std_axis_right(3)), ptp_num, 5, trial)
            
            error("The sagital axis of the left and right foot is not the same")
        end
        
        % find the longitudinal axis -> very skewed distribution because feet
        % cannot be inside the ground
        skew_left = abs(skewness(vicon_2.(trial).KinematicData.Marker.LF, 1, 1));
        skew_right = abs(skewness(vicon_2.(trial).KinematicData.Marker.RF, 1, 1));
        skew_comb = skew_left + skew_right;
        
        % find the axis with the highest skewness
        [~, idx_long] = max(skew_comb);
        
        % check if the axis is the same over trials
        if lon_axis > 0 && lon_axis ~= idx_long
            log_entry("The longitudinal axis of the left and right foot is not constant over all trials", ptp_num, 2, trial)
            
            log_entry(sprintf("Combined Skewness: %.2f %.2f %.2f", skew_comb(1), skew_comb(2), skew_comb(3)), ptp_num, 5, trial)
            
            log_entry(sprintf("Detected Skewnesses: Left: %.2f %.2f %.2f, Right: %.2f %.2f %.2f", ...
                skew_left(1), skew_left(2), skew_left(3), skew_right(1), skew_right(2), skew_right(3)), ptp_num, 5, trial)
        else
            lon_axis = idx_long;
        end
        
        % find the lateral axis -> the remaining axis
        lat_axis = setdiff(1:3, [sag_axis, lon_axis]);
        
    end
    
    if sag_axis == 0 || lon_axis == 0 || lat_axis == 0
        log_entry("The walking axes could not be determined", ptp_num, 2)
        
        error("The walking axes could not be determined")
    end
    
    % save the axis
    proc_info.walking_axes.sagittal = num2axis(sag_axis);
    proc_info.walking_axes.longitudinal = num2axis(lon_axis);
    proc_info.walking_axes.lateral = num2axis(lat_axis);
    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
    
    log_entry(sprintf("Walking axes determined: sagittal: %s, longitudinal: %s, lateral: %s", proc_info.walking_axes.sagittal, proc_info.walking_axes.longitudinal, proc_info.walking_axes.lateral), ptp_num, 5)
    
end
 
%% determine the walking speed
if ~isfield(proc_info, "walking_speed") || FIND_WALKING_SPEED
    
    if ~FIND_WALKING_SPEED
        log_entry("The walking speed has to be determined even though it was not marked (not done yet)", ptp_num, 5)
    end
    
    % load the vicon data, if not already done
    if ~exist("vicon_2", "var"); load(temp_path + "vicon_2.mat", "vicon_2"); end
    
    % get the walking axis
    try
        sag_axis = axis2num(proc_info.walking_axes.sagittal);
        lon_axis = axis2num(proc_info.walking_axes.longitudinal);
    catch ME
        log_entry("The walking axes are not defined", ptp_num, 2)
        
        % log the error
        log_entry(ME.message, ptp_num, 2)
        for i = 1:numel(ME.stack)
            log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
        end
        
        error("The walking axes are not defined")
    end
    
    % loop over all trials
    trials = string(fieldnames(vicon_2))';
    walking_direction_positions = nan(60 * 60 * 200 * 2,1); % preallocate 1h at 200Hz for both feet
    walking_direction_idx = 1;
    up_down_positions = nan(60 * 60 * 200 * 2,1); % preallocate 1h at 200Hz for both feet
    up_down_idx = 1;
    
    for trial = trials
        
        % store all positions of both feet in the walking direction and the
        % up-down direction
        numel_left = height(vicon_2.(trial).KinematicData.Marker.LF);
        numel_right = height(vicon_2.(trial).KinematicData.Marker.RF);
        walking_direction_positions(walking_direction_idx:walking_direction_idx+numel_left-1) = vicon_2.(trial).KinematicData.Marker.LF(:,sag_axis);
        walking_direction_idx = walking_direction_idx + numel_left;
        walking_direction_positions(walking_direction_idx:walking_direction_idx+numel_right-1) = vicon_2.(trial).KinematicData.Marker.RF(:,sag_axis);
        walking_direction_idx = walking_direction_idx + numel_right;
        
        up_down_positions(up_down_idx:up_down_idx+numel_left-1) = vicon_2.(trial).KinematicData.Marker.LF(:,lon_axis);
        up_down_idx = up_down_idx + numel_left;
        up_down_positions(up_down_idx:up_down_idx+numel_right-1) = vicon_2.(trial).KinematicData.Marker.RF(:,lon_axis);
        up_down_idx = up_down_idx + numel_right;
        
    end
    
    % remove nan values
    walking_direction_positions = walking_direction_positions(~isnan(walking_direction_positions));
    up_down_positions = up_down_positions(~isnan(up_down_positions));
    
    % get upwards direction
    % since the foot is always above the ground, we expect a positive or
    % negative skewness. if the skewness is positive, the upwards direction
    % is positive, otherwise negative
    if skewness(up_down_positions) > 0.5
        proc_info.walking_axes.upwards = num2direction(1);
    elseif skewness(up_down_positions) < -0.5
        proc_info.walking_axes.upwards = num2direction(-1);
    else
        log_entry("The upwards direction could not be determined", ptp_num, 2)
        
        error("The upwards direction could not be determined")
    end
    
    log_entry(sprintf("Upwards direction determined: %s %s-axis", proc_info.walking_axes.upwards, proc_info.walking_axes.longitudinal), ptp_num, 5)
    
    % calculate derivative -> speed
    speed = diff(walking_direction_positions) / 1000 * 200; % mm/5ms -> m/s
    speed = filloutliers(speed, "previous", "quartiles", "ThresholdFactor", 1.8); % 1.8 because it worked well between 1.2 and 2.4
    
    % get walking direction
    % because the foot has to accelerate and slow down, the maximum
    % velocities are higher in the walking direction. in the negative
    % walking direction, they are smaller, because there the constant
    % speed of the treadmill is measured. this results in a left- or
    % right-skewed distribution. if the distribution is left-skewed
    % (negative skew), the walking direction is in the negative
    % direction, otherwise in the positive
    if skewness(speed) < -0.5
        proc_info.walking_axes.forwards = num2direction(-1);
    elseif skewness(speed) > 0.5
        proc_info.walking_axes.forwards = num2direction(1);
    else
        log_entry("The walking direction could not be determined", ptp_num, 2)
        
        error("The walking direction could not be determined")
    end
    
    log_entry(sprintf("Walking direction determined: %s %s-axis", proc_info.walking_axes.forwards, proc_info.walking_axes.sagittal), ptp_num, 5)
    
    % find the sinister direction of the lateral axis
    % we know that a right hand system is used. the lateral axis is the
    % remaining axis.
    %
    % we form a right hand system following forward, sinister, upwards
    %
    % utilizing the cross product of the upwards and forwards direction
    % we can determine the sinister direction
    sinister_vec = cross(...
        axisdir2vec(proc_info.walking_axes.sagittal, proc_info.walking_axes.forwards), ...
        axisdir2vec(proc_info.walking_axes.longitudinal, proc_info.walking_axes.upwards)...
        ) * -1;
    
    [axis, direction] = vec2axisdir(sinister_vec);
    try
        assert(axis == lat_axis, "The sinister direction is not on the lateral axis")
    catch ME
        log_entry("The sinister direction is not the lateral axis", ptp_num, 2)
        
        % log the error
        log_entry(ME.message, ptp_num, 2)
        for i = 1:numel(ME.stack)
            log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
        end
        
        error("The sinister direction is not the lateral axis")
    end
    
    proc_info.walking_axes.sinister = num2direction(direction);
    
    log_entry(sprintf("sinister direction determined: %s %s-axis", proc_info.walking_axes.sinister, proc_info.walking_axes.lateral), ptp_num, 5)
    
    % find actual speed (maximum of the distribution)
    [f,xi] = ksdensity(speed, linspace(min(speed), max(speed), 1000));
    proc_info.walking_speed = abs(xi(f == max(f))); % m/s
    
    % save the struct
    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
    
    log_msg = sprintf("Walking speed determined: %.2f m/s, %.1f km/h", proc_info.walking_speed, proc_info.walking_speed*3.6);
    log_entry(log_msg, ptp_num, 5)
    
    % remove unnecessary variables (only large ones)
    clear walking_direction_positions up_down_positions speed f xi
    
end
 
% accumulate the distance based on the speed
if ~isfile(temp_path + "vicon_3.mat") || ACCUMULATE_DIST
    
    if ~ACCUMULATE_DIST
        log_entry("The distance has to be accumulated even though it was not marked (not done yet)", ptp_num, 5)
    end
    
    % load the vicon data, if not already done
    if ~exist("vicon_2", "var"); load(temp_path + "vicon_2.mat", "vicon_2"); end
    
    % get the walking axis (sagital axis)
    try
        axis = axis2num(proc_info.walking_axes.sagittal);
    catch ME
        log_entry("The walking axis is not defined", ptp_num, 2)
        
        % log the error
        log_entry(ME.message, ptp_num, 2)
        for i = 1:numel(ME.stack)
            log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
        end
        
        error("The walking axis is not defined")
    end
    
    % get the walking direction
    try
        walking_direction = direction2num(proc_info.walking_axes.forwards);
    catch ME
        log_entry("The walking direction is not defined", ptp_num, 2)
        
        % log the error
        log_entry(ME.message, ptp_num, 2)
        for i = 1:numel(ME.stack)
            log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
        end
        
        error("The walking direction is not defined")
    end
    
    % get the walking speed
    try
        speed = proc_info.walking_speed;
    catch ME
        log_entry("The walking speed is not defined", ptp_num, 2)
        
        % log the error
        log_entry(ME.message, ptp_num, 2)
        for i = 1:numel(ME.stack)
            log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
        end
        
        error("The walking speed is not defined")
    end
    
    % create a new struct
    vicon_3 = struct();
    
    % loop over all trials
    trials = string(fieldnames(vicon_2))';
    for trial = trials
        
        % detect if the treadmill has been started or stopped at the start or
        % the end. if so, cut the data at the start and the end
        % the thresholds have been experimentally determined analyzing the
        % fastest and slowest walking speeds
        
        % create a detector based on the standard deviation of the feet in the
        % sagital axis
        movstd_result_l = movstd(vicon_2.(trial).KinematicData.Marker.LF(:,axis), 500);
        movstd_result_r = movstd(vicon_2.(trial).KinematicData.Marker.RF(:,axis), 500);
        movstd_result = max([movstd_result_l, movstd_result_r], [], 2);
        
        % check if there's a phase where the treadmill is not running (std < 100)
        if any(movstd_result < 100)
            
            log_entry("Treadmill start or stop detected", ptp_num, 3, trial)
            
            % not walking detector (remove double detections)
            se = ones(200,1); % 1s kernel
            not_walking = movstd_result < 100;
            not_walking = imdilate(not_walking, se);
            not_walking = imerode(not_walking, se);
            not_walking = diff(not_walking); % -1: start of movement, 1: end of movement
            
            % check if there are multiple starts and stops
            if sum(not_walking == -1) > 1 || sum(not_walking == 1) > 1
                log_entry("Multiple starts or stops detected", ptp_num, 3, trial)
                
                % keep only last start in the first half and the first stop in
                % after the first start
                
                % find the last start in the first half
                length_half = floor(length(not_walking) / 2);
                start_idx = find(not_walking(1:length_half) == -1, 1, "last");
                if isempty(start_idx)
                    start_idx = 1;
                end
                
                % find the first stop after the first start
                stop_idx = find(not_walking(start_idx:end) == 1, 1, "first");
                if isempty(stop_idx)
                    stop_idx = length(not_walking);
                else
                    stop_idx = stop_idx + start_idx - 1;
                end
                
                % set all other starts and stops to 0
                not_walking(1:start_idx-1) = 0;
                not_walking(stop_idx+1:end) = 0;
            end
            
            % define the start
            if any(not_walking == -1)
                idx_start = find(not_walking == -1, 1, "first") + 1000; % 5s after the start
                
                log_entry("Beginning of the trial must be cut", ptp_num, 5, trial)
            else
                idx_start = 1;
            end
            
            % define the end
            if any(not_walking == 1)
                idx_end = find(not_walking == 1, 1, "last") - 1000; % 5s before the end
                
                log_entry("End of the trial must be cut", ptp_num, 5, trial)
            else
                idx_end = height(vicon_2.(trial).KinematicData.Marker.LF);
            end
            
            % plot the start and stop detection
            start_stop_figure = figure(Visible="off");
            plot(movstd_result)
            xline(idx_start, "g", "Start", "LineWidth", 2)
            xline(idx_end, "r", "Stop", "LineWidth", 2)
            title("Start and stop detection")
            legend("Standard deviation", "Start/Stop", "Location", "best")
            grid on;
            filename = fullfile(temp_path, string(datetime("now", "Format", "yy-MM-dd_HH-mm-ss")) + "_" + trial + "_start_stop_detection.png");
            saveas(gcf, filename);
            close(start_stop_figure)
            
            log_entry("Start and stop detection plot saved", ptp_num, 5, trial)
            
            % save the indices
            vicon_2.(trial).idx_start = idx_start;
            vicon_2.(trial).idx_end = idx_end;
            
            % clip the vicon data
            vicon_2.(trial).KinematicData.Marker.LF = vicon_2.(trial).KinematicData.Marker.LF(idx_start:idx_end,:);
            vicon_2.(trial).KinematicData.Marker.RF = vicon_2.(trial).KinematicData.Marker.RF(idx_start:idx_end,:);
            vicon_2.(trial).KinematicData.Marker.LHEE = vicon_2.(trial).KinematicData.Marker.LHEE(idx_start:idx_end,:);
            vicon_2.(trial).KinematicData.Marker.RHEE = vicon_2.(trial).KinematicData.Marker.RHEE(idx_start:idx_end,:);
            
            % clip the gait events
            vicon_2.(trial).GaitEvents.HSleftlocs(vicon_2.(trial).GaitEvents.HSleftlocs < idx_start) = [];
            vicon_2.(trial).GaitEvents.HSleftlocs(vicon_2.(trial).GaitEvents.HSleftlocs > idx_end) = [];
            vicon_2.(trial).GaitEvents.HSleftlocs = vicon_2.(trial).GaitEvents.HSleftlocs - idx_start + 1;
            vicon_2.(trial).GaitEvents.TOleftlocs(vicon_2.(trial).GaitEvents.TOleftlocs < idx_start) = [];
            vicon_2.(trial).GaitEvents.TOleftlocs(vicon_2.(trial).GaitEvents.TOleftlocs > idx_end) = [];
            vicon_2.(trial).GaitEvents.TOleftlocs = vicon_2.(trial).GaitEvents.TOleftlocs - idx_start + 1;
            vicon_2.(trial).GaitEvents.HSrightlocs(vicon_2.(trial).GaitEvents.HSrightlocs < idx_start) = [];
            vicon_2.(trial).GaitEvents.HSrightlocs(vicon_2.(trial).GaitEvents.HSrightlocs > idx_end) = [];
            vicon_2.(trial).GaitEvents.HSrightlocs = vicon_2.(trial).GaitEvents.HSrightlocs - idx_start + 1;
            vicon_2.(trial).GaitEvents.TOrightlocs(vicon_2.(trial).GaitEvents.TOrightlocs < idx_start) = [];
            vicon_2.(trial).GaitEvents.TOrightlocs(vicon_2.(trial).GaitEvents.TOrightlocs > idx_end) = [];
            vicon_2.(trial).GaitEvents.TOrightlocs = vicon_2.(trial).GaitEvents.TOrightlocs - idx_start + 1;
            
            % clip the imu data in the gait summary
            % note: this is done further below, because we don't want to load the
            % gait summary data here.
            
            % add flag to the struct
            proc_info.trials.(trial).start_stop_clipped = true;
            writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
            
        end
        
        % copy necessary data
        vicon_3.(trial).KinematicData.Marker.LF = vicon_2.(trial).KinematicData.Marker.LF / 1000; % mm -> m
        vicon_3.(trial).KinematicData.Marker.RF = vicon_2.(trial).KinematicData.Marker.RF / 1000; % mm -> m
        vicon_3.(trial).KinematicData.Marker.LHEE = vicon_2.(trial).KinematicData.Marker.LHEE / 1000; % mm -> m
        vicon_3.(trial).KinematicData.Marker.RHEE = vicon_2.(trial).KinematicData.Marker.RHEE / 1000; % mm -> m
        vicon_3.(trial).GaitEvents = vicon_2.(trial).GaitEvents;
        
        % calculate the distance based on the speed
        dist = cumsum(speed * 0.005 * ones(height(vicon_2.(trial).KinematicData.Marker.LF),1)); %distance in meters
        
        % add to the position data
        vicon_3.(trial).KinematicData.Marker.LF(:,axis) = vicon_3.(trial).KinematicData.Marker.LF(:,axis) + (dist * walking_direction);
        vicon_3.(trial).KinematicData.Marker.RF(:,axis) = vicon_3.(trial).KinematicData.Marker.RF(:,axis) + (dist * walking_direction);
        
    end
    
    % save the struct
    save(temp_path + "vicon_3.mat", "vicon_3")
    
    log_entry("Distance accumulated", ptp_num)
    
    % remove unnecessary variables
    clear vicon_2 dist axis walking_direction
    
end
 
%% find IMU data start and end
% in sum_zm_l.synced.trim(i).acc and sum_zm_r.synced.trim(i).acc the already
% trimmed and processed data of the IMU is stored. Since we do not want to use
% preprocessed data, we sync this preprocessed data to the raw data in zm_l_raw and
% zm_r_raw. This can be done with alignsignals.
 
if ~isfile(fullfile(temp_path, "zm_aligned.mat")) || ~isfile(fullfile(temp_path, "vicon_4.mat")) || ALIGN_IMU_VICON
    
    % load the vicon data, if not already done
    if ~exist("vicon_3", "var"); load(temp_path + "vicon_3.mat","vicon_3"); end
    
    if contains(ptp_folder, "SonE")
        
        
        if ~ALIGN_IMU_VICON
            log_entry("The IMU data has to be aligned even though it was not marked (not done yet)", ptp_num, 5)
        end
        
        % load the raw data, if not already done
        if ~exist("zm_l_raw", "var"); load(temp_path + "zm_l_raw.mat", "zm_l_raw"); end
        if ~exist("zm_r_raw", "var"); load(temp_path + "zm_r_raw.mat", "zm_r_raw"); end
        
        % load the gait summary data, if not already done
        if ~exist("sum_zm_l", "var"); load(temp_path + "sum_zm_l.mat", "sum_zm_l"); end
        if ~exist("sum_zm_r", "var"); load(temp_path + "sum_zm_r.mat", "sum_zm_r"); end
        
        % create a new struct
        zm_aligned = struct();
        
        % loop over all trials
        trials = string(fieldnames(vicon_3))';
        for trial = trials
            
            % find the row where the trial is
            trial_row = find(string({sum_zm_l.synced.trim.Trials}) == trial);
            
            % if the data had to be clipped, the gait summary data has to be clipped as well
            % note: only the acc data is clipped, because the other data is not used,
            % however, this (c|sh)ould be done as well
            if isfield(vicon_3.(trial), "idx_start") || isfield(vicon_3.(trial), "idx_end")
                sum_zm_l.synced.trim(trial_row).acc = sum_zm_l.synced.trim(trial_row).acc(vicon_3.(trial).idx_start:vicon_3.(trial).idx_end,:);
                sum_zm_r.synced.trim(trial_row).acc = sum_zm_r.synced.trim(trial_row).acc(vicon_3.(trial).idx_start:vicon_3.(trial).idx_end,:);
            end
            
            % left
            idx_start_l = finddelay(vecnorm(sum_zm_l.synced.trim(trial_row).acc, 2, 2), vecnorm(zm_l_raw.tt.acc, 2, 2)) + 1;
            
            % right
            idx_start_r = finddelay(vecnorm(sum_zm_r.synced.trim(trial_row).acc, 2, 2), vecnorm(zm_r_raw.tt.acc, 2, 2)) + 1;
            
            % clip zm data TODO: the clipping should only be done after the check if the max correlation
            zm_aligned.(trial).l = zm_l_raw;
            zm_aligned.(trial).l.tt = zm_aligned.(trial).l.tt(idx_start_l:idx_start_l+height(sum_zm_l.synced.trim(trial_row).acc)-1,:);
            zm_aligned.(trial).l.t_start = zm_aligned.(trial).l.tt.t(1);
            
            zm_aligned.(trial).r = zm_r_raw;
            zm_aligned.(trial).r.tt = zm_aligned.(trial).r.tt(idx_start_r:idx_start_r+height(sum_zm_r.synced.trim(trial_row).acc)-1,:);
            zm_aligned.(trial).r.t_start = zm_aligned.(trial).r.tt.t(1);
            
            % check if the data is properly aligned
            
            corr_l_x = corr(sum_zm_l.synced.trim(trial_row).acc(:,1), zm_aligned.(trial).l.tt.acc(:,1));
            corr_l_y = corr(sum_zm_l.synced.trim(trial_row).acc(:,2), zm_aligned.(trial).l.tt.acc(:,2));
            corr_l_z = corr(sum_zm_l.synced.trim(trial_row).acc(:,3), zm_aligned.(trial).l.tt.acc(:,3));
            corr_r_x = corr(sum_zm_r.synced.trim(trial_row).acc(:,1), zm_aligned.(trial).r.tt.acc(:,1));
            corr_r_y = corr(sum_zm_r.synced.trim(trial_row).acc(:,2), zm_aligned.(trial).r.tt.acc(:,2));
            corr_r_z = corr(sum_zm_r.synced.trim(trial_row).acc(:,3), zm_aligned.(trial).r.tt.acc(:,3));
            try
                assert(corr_l_x > 0.95, "The x-axis of the left IMU data is not properly aligned")
                assert(corr_l_y > 0.95, "The y-axis of the left IMU data is not properly aligned")
                assert(corr_l_z > 0.95, "The z-axis of the left IMU data is not properly aligned")
                assert(corr_r_x > 0.95, "The x-axis of the right IMU data is not properly aligned")
                assert(corr_r_y > 0.95, "The y-axis of the right IMU data is not properly aligned")
                assert(corr_r_z > 0.95, "The z-axis of the right IMU data is not properly aligned")
            catch ME
                log_entry(ME.message, ptp_num, 2, trial)
                
                % check if the left and right sensor were swapped
                corr_l_x_sw = corr(sum_zm_l.synced.trim(trial_row).acc(:,1), zm_aligned.(trial).r.tt.acc(:,1));
                corr_l_y_sw = corr(sum_zm_l.synced.trim(trial_row).acc(:,2), zm_aligned.(trial).r.tt.acc(:,2));
                corr_l_z_sw = corr(sum_zm_l.synced.trim(trial_row).acc(:,3), zm_aligned.(trial).r.tt.acc(:,3));
                corr_r_x_sw = corr(sum_zm_r.synced.trim(trial_row).acc(:,1), zm_aligned.(trial).l.tt.acc(:,1));
                corr_r_y_sw = corr(sum_zm_r.synced.trim(trial_row).acc(:,2), zm_aligned.(trial).l.tt.acc(:,2));
                corr_r_z_sw = corr(sum_zm_r.synced.trim(trial_row).acc(:,3), zm_aligned.(trial).l.tt.acc(:,3));
                
                if corr_l_x_sw > 0.95 && corr_l_y_sw > 0.95 && corr_l_z_sw > 0.95 && corr_r_x_sw > 0.95 && corr_r_y_sw > 0.95 && corr_r_z_sw > 0.95
                    log_entry("The left and right IMU data were swapped", ptp_num, 3, trial)
                    
                    % add flag to the struct
                    proc_info.trials.(trial).imu_swapped = true;
                    
                else
                    
                    % remove the trial from the structs
                    zm_aligned = rmfield(zm_aligned, trial);
                    vicon_3 = rmfield(vicon_3, trial);
                    
                    log_entry("Trial removed", ptp_num, 3, trial)
                    
                    corr_min = min([corr_l_x, corr_l_y, corr_l_z, corr_r_x, corr_r_y, corr_r_z]);
                    corr_min_sw = min([corr_l_x_sw, corr_l_y_sw, corr_l_z_sw, corr_r_x_sw, corr_r_y_sw, corr_r_z_sw]);
                    
                    log_entry(sprintf("Minimum correlation: %.2f", corr_min), ptp_num, 5, trial)
                    log_entry(sprintf("Minimum correlation swapped: %.2f", corr_min_sw), ptp_num, 5, trial)
                    
                    % add flag to the struct
                    proc_info.trials.(trial).imu_raw_data_missmatch = true;
                    proc_info.trials.(trial).removal_reason = "IMU raw data missmatch";
                    
                    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                    
                    continue
                end
            end
            
            %% check the alignment of the VICON and IMU data
            % the data should be aligned in time. also the number of samples should be the
            % same. here, first the number of samples is checked and then the periods of
            % movement are compared.
            
            % check number of samples
            d_samples_l = height(vicon_3.(trial).KinematicData.Marker.LF) - height(zm_aligned.(trial).l.tt.acc);
            d_samples_r = height(vicon_3.(trial).KinematicData.Marker.RF) - height(zm_aligned.(trial).r.tt.acc);
            try
                assert(d_samples_l <= 1, "The number of samples of the left IMU data and the VICON data is not the same")
                assert(d_samples_r <= 1, "The number of samples of the right IMU data and the VICON data is not the same")
            catch ME
                log_entry(ME.message, ptp_num, 3, trial)
                
                % add flag to the struct
                proc_info.trials.(trial).imu_vicon_num_of_samples_missmatch = true;
                writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                
                if d_samples_l > 5 || d_samples_r > 5
                    % remove the trial from the structs
                    zm_aligned = rmfield(zm_aligned, trial);
                    vicon_3 = rmfield(vicon_3, trial);
                    
                    log_entry("The samples difference is too high. The trial was removed.", ptp_num, 2, trial)
                    
                    log_entry(sprintf("Difference in samples: left %d, right %d", d_samples_l, d_samples_r), ptp_num, 5, trial)
                    
                    % add removal reason
                    proc_info.trials.(trial).removal_reason = "IMU VICON num of samples missmatch";
                    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                    
                    continue
                    
                else
                    log_entry("The samples difference is small. The trial was not removed.", ptp_num, 4, trial)
                    
                    log_entry(sprintf("Difference in samples: left %d, right %d", d_samples_l, d_samples_r), ptp_num, 5, trial)
                    
                end
                
            end
            
            % check periods of movement
            % here the following is done:
            % VICON:
            % - the speed of movement is calculated from the position data
            % - the speed is filtered with a lowpass filter
            % - take the norm of the speed
            % - the periods of movement are detected with a threshold
            %
            % IMU:
            % - the magnitude of the angular velocity is calculated
            % - the magnitude is filtered with a lowpass filter
            % - the norm of the angular velocity is taken
            % - the periods of movement are detected with a threshold
            
            % VICON
            % speed
            v_vicon_l = diff(vicon_3.(trial).KinematicData.Marker.LF) / 0.005; % m/s
            v_vicon_r = diff(vicon_3.(trial).KinematicData.Marker.RF) / 0.005; % m/s
            
            % lowpass filter
            [b, a] = butter(2, 0.1);
            v_vicon_l = filtfilt(b, a, v_vicon_l);
            v_vicon_r = filtfilt(b, a, v_vicon_r);
            
            % norm
            v_vicon_l = vecnorm(v_vicon_l,2,2);
            v_vicon_r = vecnorm(v_vicon_r,2,2);
            
            % threshold
            th_vicon = 0.06; % m/s
            vicon_periods_l = v_vicon_l > th_vicon;
            vicon_periods_r = v_vicon_r > th_vicon;
            
            % IMU
            % angular velocity
            w_l = zm_aligned.(trial).l.tt.gyr - mean(zm_aligned.(trial).l.tt.gyr);
            w_r = zm_aligned.(trial).r.tt.gyr - mean(zm_aligned.(trial).r.tt.gyr);
            
            % lowpass filter
            [b, a] = butter(2, 0.1);
            w_l = filtfilt(b, a, w_l);
            w_r = filtfilt(b, a, w_r);
            
            % norm
            w_l = vecnorm(w_l,2,2);
            w_r = vecnorm(w_r,2,2);
            
            % threshold
            th_imu = 65; % deg/s
            imu_periods_l = w_l > th_imu;
            imu_periods_r = w_r > th_imu;
            
            % check if the periods align
            l_min = min([height(vicon_periods_l), height(imu_periods_l)]);
            r_min = min([height(vicon_periods_r), height(imu_periods_r)]); % should be the same, but just in case
            
            corr_l = corr(vicon_periods_l(1:l_min), imu_periods_l(1:l_min));
            corr_r = corr(vicon_periods_r(1:r_min), imu_periods_r(1:r_min));
            try
                assert(corr_l > 0.75, "The periods of movement of the left IMU data and the VICON data do not align")
                assert(corr_r > 0.75, "The periods of movement of the right IMU data and the VICON data do not align")
            catch ME
                log_entry(ME.message, ptp_num, 3, trial)
                
                % check if left and right are swapped
                corr_l_sw = corr(vicon_periods_l(1:l_min), imu_periods_r(1:l_min));
                corr_r_sw = corr(vicon_periods_r(1:r_min), imu_periods_l(1:r_min));
                
                if corr_l_sw > 0.5 && corr_r_sw > 0.5
                    log_entry("The left and right IMU data were swapped", ptp_num, 3, trial)
                    
                    this_zm_l = zm_aligned.(trial).l;
                    this_zm_r = zm_aligned.(trial).r;
                    
                    zm_aligned.(trial).l = this_zm_r;
                    zm_aligned.(trial).r = this_zm_l;
                    
                    corr_l = corr_l_sw;
                    corr_r = corr_r_sw;
                    
                end
                
                % add flag to the struct
                proc_info.trials.(trial).imu_vicon_periods_missmatch = true;
                writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                
                if corr_l < 0.5 || corr_r < 0.5
                    
                    % remove the trial from the structs
                    zm_aligned = rmfield(zm_aligned, trial);
                    vicon_3 = rmfield(vicon_3, trial);
                    
                    log_entry("The periods of movement do not align. The trial was removed.", ptp_num, 2, trial)
                    
                    log_entry(sprintf("Correlation of periods: left %.2f, right %.2f", corr_l, corr_r), ptp_num, 5, trial)
                    
                    % add removal reason
                    proc_info.trials.(trial).removal_reason = "movement periods missmatch";
                    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                    
                    continue
                    
                else
                    
                    log_entry("The periods of movement do not perfectly align, but we'll try.", ptp_num, 4, trial)
                    
                    log_entry(sprintf("Correlation of periods: left %.2f, right %.2f", corr_l, corr_r), ptp_num, 5, trial)
                    
                end
                
                % save a plot of the periods
                period_figure = figure(Visible="on");
                subplot(2,1,1)
                plot(vicon_periods_l(1:l_min), '.', "MarkerSize", 10)
                hold on
                plot(imu_periods_l(1:l_min), '.', "MarkerSize", 10)
                xlim([length(vicon_periods_l)/2 - 500, length(vicon_periods_l)/2 + 499]);
                title(sprintf("Left foot, correlation: %.2f", corr_l))
                legend("VICON", "IMU", "Location", "best")
                grid on;
                subplot(2,1,2)
                plot(vicon_periods_r(1:r_min), '.', "MarkerSize", 10)
                hold on
                plot(imu_periods_r(1:r_min), '.', "MarkerSize", 10)
                xlim([length(vicon_periods_r)/2 - 500, length(vicon_periods_r)/2 + 499]);
                title(sprintf("Right foot, correlation: %.2f", corr_r))
                legend("VICON", "IMU", "Location", "best")
                grid on;
                filename = fullfile(temp_path, string(datetime("now", "Format", "yy-MM-dd_HH-mm-ss")) + "_" + trial + "_periods.fig");
                saveas(period_figure, filename);
                close(period_figure)
                
                log_entry("Periods plot saved", ptp_num, 5, trial)
            end
            
            log_entry("IMU data aligned", ptp_num, 5, trial)
            
            
        end
        
    else
        
        % load the already aligned data
        zm_ankle_l = load(fullfile(ptp_folder, "Processed", "ZM", "calibrated_ankle_L.mat"));
        zm_ankle_r = load(fullfile(ptp_folder, "Processed", "ZM", "calibrated_ankle_R.mat"));
        
        % create a new struct
        zm_aligned = struct();
        
        % the only trial is "Walk_Norm"
        trial = "Walk_Norm";
        
        % left
        zm_aligned.(trial).l = zm_ankle_l;
        zm_aligned.(trial).r = zm_ankle_r;
        
    end
    
    % save the aligned data
    save(fullfile(temp_path, "zm_aligned.mat"), "zm_aligned")
    vicon_4 = vicon_3;
    save(fullfile(temp_path, "vicon_4.mat"), "vicon_4")
    log_entry("IMU data aligned", ptp_num)
    
    % remove unnecessary variables
    clear vicon_3 sum_zm_l sum_zm_r zm_l_raw zm_r_raw trial_row idx_start_l
    clear idx_start_r corr_l_x corr_l_y corr_l_z corr_r_x corr_r_y corr_r_z
    clear d_samples_l d_samples_r b a v_vicon_l v_vicon_r th_vicon
    clear vicon_periods_l vicon_periods_r w_l w_r th_imu imu_periods_l
    clear imu_periods_r l_min r_min corr_l corr_r
end
 
%% filter the data
 
% load the proc info
proc_info = readstruct(fullfile(ptp_folder, "proc_info.json"));
 
 
if ~isfile(fullfile(temp_path, "processed_data.mat")) || FILTER
    
    if ~FILTER
        log_entry("The data has to be filtered even though it was not marked (not done yet, 3)", ptp_num, 5)
    end
    
    % flag as not done
    writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"));
    
    % load the aligned data, if not already done
    if ~exist("zm_aligned", "var"); load(temp_path + "zm_aligned.mat", "zm_aligned"); end
    
    % create struct to store the filtered data
    filtered = struct();
    
    % check if at least one trial is available
    if ~isempty(fieldnames(zm_aligned))
        
        if ~exist("vicon_4", "var"); load(temp_path + "vicon_4.mat", "vicon_4"); end
        
        % get the walking speed
        try
            speed = proc_info.walking_speed;
        catch ME
            log_entry("The walking speed is not defined", ptp_num, 2)
            
            % log the error
            log_entry(ME.message, ptp_num, 2)
            for i = 1:numel(ME.stack)
                log_entry(ME.stack(i).name + " (Line: " + string(ME.stack(i).line) + ")", ptp_num, 2)
            end
            
            error("The walking speed is not defined")
        end
        
        %% instance of filter
        ekf = EKF1;
        dt = 0.005;
        
        % import tuning parameters, if available
        if isfield(options, "ekf_tuning_file")
            assert(ekf.importTuning(options.ekf_tuning_file), "The tuning parameters could not be imported")
            log_entry("Tuning parameters imported", ptp_num, 5)
        end
        
        % loop over all trials
        trials = string(fieldnames(zm_aligned))';
        
        % check if the yaw angles are already defined
        % the yaw angles are used to calculate the initial position. also,
        % during the filtering process, the filtered angle is nudged
        % towards this angle. is is ASSUMED that the mean yaw angle of each
        % foot is constant over all the trials and during each trial.
        % please note also, that this angle does not necessarly reflect the
        % actual toe out angle since the JUMP sensor might nit be perfectly
        % aligned with the feet.
        %
        % if no yaw angle is determined yet, the walking trial with the
        % longest duration is filtered without x and y position nudging.
        % therefore the trajectory will drift in a directrion. based on
        % this drift, the mean yaw angle can be calculated which is then
        % used for the rest of the trials.
        find_yaw = false;
        if ~isfield(proc_info, "left_yaw") || isnan(proc_info.left_yaw) || ~isfield(proc_info, "right_yaw") || isnan(proc_info.right_yaw) || FIND_YAW
            % find the trial with the longest duration
            durations = zeros(1, numel(trials));
            for i = 1:numel(trials)
                durations(i) = height(zm_aligned.(trials(i)).l.tt);
            end
            [~, idx] = max(durations);
            trials = [trials(idx), trials];
            
            log_msg = sprintf("No yaw angle is defined yet. The trial %s is filtered first to determine the yaw angle", trials(1));
            log_entry(log_msg, ptp_num, 4);
            
            find_yaw = true;
        end
        
        for trial = trials
            
            t_before_filter = tic;
            
            %% left foot
            
            % zero velocity detection
            zv_l = jumpZeroVelocityDetection(zm_aligned.(trial).l, "threshold_acc", 1.2, "threshold_gyr", 15, "windowSize", 30, "plot", false);
            
            % start and end during stance
            idx_l(1) = find(zv_l, 1, "first");
            idx_l(2) = find(zv_l, 1, "last");
            
            % initial values
            % initial gyroscope bias
            w_zv_l = zm_aligned.(trial).l.tt.gyr(zv_l,:);
            ekf.b_w_0 = [mean(w_zv_l(1:20,1)), mean(w_zv_l(1:20,2)), mean(w_zv_l(1:20,3))] * pi/180;
            
            % measurements
            ekf.zv = zv_l(idx_l(1):idx_l(2));
            ekf.z_acc = zm_aligned.(trial).l.tt.acc(idx_l(1):idx_l(2),:)';
            ekf.z_vel = zeros(size(ekf.z_acc));
            ekf.z_gyr = zm_aligned.(trial).l.tt.gyr(idx_l(1):idx_l(2),:)' * pi/180; % convert to rad/s
            ekf.z_pos = zeros(size(ekf.z_acc));
            
            % check if the left foot yaw angle is already defined
            if find_yaw
                proc_info.left_yaw = nan;
                
                yaw = 0; % rad
                
            else
                yaw = proc_info.left_yaw; % rad
            end
            
            % initial orientation
            %> @remark: there's something odd with the way I calculate the initial
            %> orientation. Since the trust in this initial value is very low and
            %> the results are pretty well, I will leave it as it is for the moment.
            %> But do not take this as a good example.
            
            % mean accelerations during stance
            a_zv_l = zm_aligned.(trial).l.tt.acc(zv_l,:);
            a_x = mean(a_zv_l(1:20,1));
            a_y = mean(a_zv_l(1:20,2));
            a_z = mean(a_zv_l(1:20,3));
            % calculate initial orientation
            roll = atan2(a_y, a_z);
            pitch = atan2(-a_x, sqrt(a_y^2 + a_z^2));
            rot_vec = [roll, pitch, yaw];
            ekf.q_0 = quaternion(rot_vec, "rotvec").normalize.compact';
            
            log_entry(sprintf("Initial orientation left: roll %.1f°, pitch %.1f°, yaw %.1f°", rad2deg(roll), rad2deg(pitch), rad2deg(yaw)), ptp_num, 5, trial)
            
            % pseudo position update
            ekf.z_pos(1,92:end) = cumsum(speed * dt * ones(size(ekf.z_pos(1,92:end))));
            if find_yaw
                % disable x and y position pseudo updates if yaw must be found
                ekf.z_pos(1,:) = nan; %this line turns off x-position pseudo updates
                ekf.z_pos(2,:) = nan; %this line turns off y-position pseudo updates
                
                log_entry("Pseudo position updates disabled", ptp_num, 5, trial)
            end
            
            % pseudo yaw update
            ekf.z_yaw = yaw * ones(1, size(ekf.z_acc,2));
            
            % Run left filter
            t_before_filter_l = tic;
            smoothed_l = ekf.run(plot=false);
            log_entry(sprintf("Filtering left foot took %.2f seconds", toc(t_before_filter_l)), ptp_num, 5, trial)
            
            %% right foot
            
            % zero velocity detection
            zv_r = jumpZeroVelocityDetection(zm_aligned.(trial).r, "threshold_acc", 1.2, "threshold_gyr", 15, "windowSize", 30, "plot", false);
            
            % start and end during stance
            idx_r(1) = find(zv_r, 1, "first");
            idx_r(2) = find(zv_r, 1, "last");
            
            % initial values
            % initial gyroscope bias
            w_zv_r = zm_aligned.(trial).r.tt.gyr(zv_r,:);
            ekf.b_w_0 = [mean(w_zv_r(1:20,1)), mean(w_zv_r(1:20,2)), mean(w_zv_r(1:20,3))] * pi/180;
            
            % check if the right foot yaw angle is already defined
            if find_yaw
                proc_info.right_yaw = nan;
                
                yaw = 0; % rad
                
            else
                yaw = proc_info.right_yaw; % rad
            end
            
            % initial orientation
            % mean accelerations during stance
            a_x = mean(zm_aligned.(trial).r.tt.acc(38:104,1));
            a_y = mean(zm_aligned.(trial).r.tt.acc(38:104,2));
            a_z = mean(zm_aligned.(trial).r.tt.acc(38:104,3));
            % calculate initial orientation
            roll = atan2(a_y, a_z);
            pitch = atan2(-a_x, sqrt(a_y^2 + a_z^2));
            rot_vec = [roll, pitch, yaw];
            ekf.q_0 = quaternion(rot_vec, "rotvec").normalize.compact';
            
            log_entry(sprintf("Initial orientation right: roll %.1f°, pitch %.1f°, yaw %.1f°", rad2deg(roll), rad2deg(pitch), rad2deg(yaw)), ptp_num, 5, trial)
            
            % measurements
            ekf.zv = zv_r(idx_r(1):idx_r(2));
            ekf.z_acc = zm_aligned.(trial).r.tt.acc(idx_r(1):idx_r(2),:)';
            ekf.z_vel = zeros(size(ekf.z_acc));
            ekf.z_gyr = zm_aligned.(trial).r.tt.gyr(idx_r(1):idx_r(2),:)' * pi/180; % convert to rad/s
            ekf.z_pos = zeros(size(ekf.z_acc));
            
            % pseudo position update
            ekf.z_pos(1,92:end) = cumsum(speed * dt * ones(size(ekf.z_pos(1,92:end))));
            if find_yaw
                % disable x and y position pseudo updates if yaw must be found
                ekf.z_pos(1,:) = nan; %this line turns off x-position pseudo updates
                ekf.z_pos(2,:) = nan; %this line turns off y-position pseudo updates
                
                log_entry("Pseudo position updates disabled", ptp_num, 5, trial)
            end
            
            % pseudo yaw update
            ekf.z_yaw = yaw * ones(1, size(ekf.z_acc,2));
            
            % Run right filter
            t_before_filter_r = tic;
            smoothed_r = ekf.run(plot=false);
            log_entry(sprintf("Filtering right foot took %.2f seconds", toc(t_before_filter_r)), ptp_num, 5, trial)
            
            % save the yaw angles and dismiss the filtered data
            if find_yaw
                
                % calculate the mean yaw angle
                % left
                end_x = smoothed_l(8, end);
                end_y = smoothed_l(9, end);
                proc_info.left_yaw = -1 * atan2(end_y, end_x);
                
                % right
                end_x = smoothed_r(8, end);
                end_y = smoothed_r(9, end);
                proc_info.right_yaw = -1 * atan2(end_y, end_x);
                
                % save the struct
                writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
                
                find_yaw = false;
                
                log_msg = sprintf("Yaw angles determined: left: %.1f °, right: %.1f °", rad2deg(proc_info.left_yaw), rad2deg(proc_info.right_yaw));
                log_entry(log_msg, ptp_num)
                
                continue
            end
            
            % extract trajectories
            traj_l = smoothed_l(8:10,:)';
            traj_r = smoothed_r(8:10,:)';
            
            % add start and end (match original length)
            filtered.(trial).l = [nan(idx_l(1)-1,3); traj_l; nan(height(zm_aligned.(trial).r.tt)-idx_l(2),3)];
            filtered.(trial).r = [nan(idx_r(1)-1,3); traj_r; nan(height(zm_aligned.(trial).r.tt)-idx_r(2),3)];
            
            log_entry("Data filtered", ptp_num, 4, trial)
            
            log_entry(sprintf("Filtering took %.2f seconds", toc(t_before_filter)), ptp_num, 5, trial)
            
            % add filtered flag to the struct
            proc_info.trials.(trial).filtered = true;
            writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
            
        end
        
        %% save the data
        result = struct();
        result.jump_filtered = filtered;
        result.vicon = vicon_4;
        save(fullfile(temp_path, "processed_data.mat"), "result")
        
        % write the EKF parameters to the log file
        log_entry("EKF parameters", ptp_num, 5)
        log_entry(sprintf("P_0_q: %.1e", ekf.P_0_q), ptp_num, 5)
        log_entry(sprintf("P_0_v: %.1e", ekf.P_0_v), ptp_num, 5)
        log_entry(sprintf("P_0_p: %.1e", ekf.P_0_p), ptp_num, 5)
        log_entry(sprintf("P_0_b_w: %.1e", ekf.P_0_b_w), ptp_num, 5)
        log_entry(sprintf("var_acc_xy: %.1e", ekf.var_acc_xy), ptp_num, 5)
        log_entry(sprintf("var_acc_z: %.1e", ekf.var_acc_z), ptp_num, 5)
        log_entry(sprintf("var_vel_ZUPT: %.1e", ekf.var_vel_ZUPT), ptp_num, 5)
        log_entry(sprintf("var_gyr: %.1e", ekf.var_gyr), ptp_num, 5)
        log_entry(sprintf("var_pos_x: %.1e", ekf.var_pos_x), ptp_num, 5)
        log_entry(sprintf("var_pos_y: %.1e", ekf.var_pos_y), ptp_num, 5)
        log_entry(sprintf("var_pos_z: %.1e", ekf.var_pos_z), ptp_num, 5)
        log_entry(sprintf("var_yaw_pseudo: %.1e", ekf.var_yaw_pseudo), ptp_num, 5)
        log_entry(sprintf("sigma_w_p: %.1e", ekf.sigma_w_p), ptp_num, 5)
        log_entry(sprintf("var_a: %.1e", ekf.var_a), ptp_num, 5)
        log_entry(sprintf("var_b_w: %.1e", ekf.var_b_w), ptp_num, 5)
        
    else
        log_entry("No trials available", ptp_num, 2)
    end
    
end
 
log_entry("Processing completed", ptp_num);
 
% set the status to completed
proc_info.filtering_done = true;
proc_info.filtering_done_time = datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss');
writestruct(proc_info, fullfile(ptp_folder, "proc_info.json"))
 
 
end