%% Gait Variability Analysis CLBP

% Gait Variability Analysis 
% Script created for MAP 2020-2021
% adapted from Claudine Lamoth and Iris Hagoort
% version1 October 2020

% Input: needs mat file which contains all raw accelerometer data
% Input: needs excel file containing the participant information including
% leg length.
%% Clear and close;

clear;
close all;
%% Load data;
% Select 1 trial. 
% For loop to import all data will be used at a later stage

[FNaam,FilePad] = uigetfile('*.xls','Load phyphox data...');
filename =[FilePad FNaam];
PhyphoxData = xlsread(filename)

%load('Phyphoxdata.mat'); % loads accelerometer data, is stored in struct with name AccData
%load('ExcelInfo.mat');
%Participants = fields(AccData);
%% Settings;
%adapted from GaitOutcomesTrunkAccFuncIH

LegLength = 98;                             % LegLength info not available!
FS = 100;                                   % Sample frequency

Gr = 9.81;                                  % Gravity acceleration, multiplication factor for accelerations
StrideFreqEstimate = 1.00;                  % Used to set search for stride frequency from 0.5*StrideFreqEstimate until 2*StrideFreqEstimate
StrideTimeRange = [0.2 4.0];                % Range to search for stride time (seconds)
IgnoreMinMaxStrides = 0.10;                 % Number or percentage of highest&lowest values ignored for improved variability estimation
N_Harm = 12;                                % Number of harmonics used for harmonic ratio, index of harmonicity and phase fluctuation
LowFrequentPowerThresholds = ...
    [0.7 1.4];                              % Threshold frequencies for estimation of low-frequent power percentages
Lyap_m = 7;                                 % Embedding dimension (used in Lyapunov estimations)
Lyap_FitWinLen = round(60/100*FS);          % Fitting window length (used in Lyapunov estimations Rosenstein's method)
Sen_m = 5;                                  % Dimension, the length of the subseries to be matched (used in sample entropy estimation)
Sen_r = 0.3;                                % Tolerance, the maximum distance between two samples to qualify as match, relative to std of DataIn (used in sample entropy estimation)
NStartEnd = [100];
M  = 5;                                     % Maximum template length
ResultStruct = struct();                    % Empty struct


inputData = (PhyphoxData(:,[1 2 3 4]));      % matrix with accelerometer data
ApplyRealignment = true;
ApplyRemoveSteps = true;
WindowLen = FS*10;


%% Filter and Realign Accdata
% dataAcc depends on ApplyRealignment = true/false
% dataAcc_filt (low pass Butterworth Filter + zerophase filtering

[dataAcc, dataAcc_filt] = FilterandRealignFunc(inputData,FS,ApplyRealignment);

%% Step dectection
% Determines the number of steps in the signal so that the first 30 and last 30 steps in the signal can be removed
% StrideTimeSamples is needed for calculation stride parameters!  

[dataAccCut,dataAccCut_filt,StrideTimeSamples] = StepDetectionFunc(FS,ApplyRemoveSteps,dataAcc,dataAcc_filt,StrideTimeRange);

%% Calculate Stride Parameters
% Outcomeparameters: StrideRegularity,RelativeStrideVariability,StrideTimeSamples,StrideTime

[ResultStruct] = CalculateStrideParametersFunc(dataAccCut_filt,FS,ApplyRemoveSteps,dataAcc,dataAcc_filt,StrideTimeRange);

%% Calculate spatiotemporal stride parameters
% Measures from height variation by double integration of VT accelerations and high-pass filtering
% StepLengthMean; Distance;  WalkingSpeedMean;  StrideTimeVariability; StrideSpeedVariability; 
% StrideLengthVariability; StrideTimeVariabilityOmitOutlier; StrideSpeedVariabilityOmitOutlier; StrideLengthVariabilityOmitOutlier;

[ResultStruct] = SpatioTemporalGaitParameters(dataAccCut_filt,StrideTimeSamples,ApplyRealignment,LegLength,FS,IgnoreMinMaxStrides,ResultStruct);

%% Measures derived from spectral analysis
% IndexHarmonicity_V/ML/AP/ALL ; HarmonicRatio_V/ML/AP ; HarmonicRatioP_V/ML/AP ; FrequencyVariability_V/ML/AP;  Stride Frequency

AccVectorLen = sqrt(sum(dataAccCut_filt(:,1:3).^2,2));
[ResultStruct] = SpectralAnalysisGaitfunc(dataAccCut_filt,WindowLen,FS,N_Harm,LowFrequentPowerThresholds,AccVectorLen,ResultStruct);

%% calculate non-linear parameters
% Outcomeparameters: Sample Entropy, Lyapunov exponents

[ResultStruct] = CalculateNonLinearParametersFunc(ResultStruct,dataAccCut,WindowLen,FS,Lyap_m,Lyap_FitWinLen,Sen_m,Sen_r);

% Save struct as .mat file
% save('GaitVarOutcomesParticipantX.mat', 'OutcomesAcc');







%% AggregateFunction (seperate analysis per minute);
% see AggregateEpisodeValues.m
% 
%