%%     Gait Variability Analysis CLBP
%% Script created for MAP 2020-2021 LD
%% Adapted from Claudine Lamoth and Iris Hagoort
%% Version 4 April 2021
%
% Input: needs .mat file which contains all raw accelerometer data
% Input: needs excel file containing the participant information (including leg length, distance covered)
%
% Attention:
% Pre-processing CLBP and controls data dissimilar (substraction gravity, order, detrend)
% CLBP Participant 6 was initially deleted from analysis (no leglength available), has been added to AccData
% Input data should be at least 10 sec. long! (see if endPos-startPos < 10*FS statement)
% Start and End position depend on the turns detected within the minute!
% if ApplyRemoveSteps = true; first and last step will be deleted. Therefore Windowlen = size(dataAccCut,1) changed in GaitOutcomesTrunkAccFuncIH
% Data can be found in a seperate repository called / Data_LBP_Lars
% Please set path to these .mat files
% P35 was excluded in the previous study! P60 new participant! Check if path = correct!
%
%% Clear and close;

clearvars; close all;

%
%% Load data;

%% Load data
load('AccData_LD.mat');                     % Loads accelerometer data, is stored in struct with name AccData
load('GeneralInformation_LD.mat');

%% Settings;

Participants = fields(AccData);             % Raw Acc Data (not reorderd, detrended, et cetera.)
FS = 100;                                   % Sample frequency
LegLengthsM = GeneralInformation.LegLength; % In meter!
ResultStruct = struct();                    % Empty struct
Duration = 6;                               % Duration of measurement in minutes

plotit = true;                              % If true, plot turn detection

%
%% Check allignment of sensor for data patients
% First check VT direction --> if this is flipped over, invert AP signal as well!
% Expect problems with the ML direction.
% Note: patient 6 / 7 / 8 are measured WITHOUT gravity! --> problems with realignment
% ML : (-) = left (+) = right
%
% For patient 4: VT signal looks odd, although AP direction is correct.
% For Patient 6: AP signal looks odd, although VT direction is correct.
%
%% Realign and Detrend the data
% Apply realignment as described in Rispens S, Pijnappels M, van Schooten K, Beek PJ, Daffertshofer A, van Die?n JH (2014).
% ''Gait characteristics are robust against limited repositioning error of sensors at the lumbar spine, in particular if our offline realignment is applied''
% (1) First realignment, (2) detrend, not vice versa --> because otherwise realignment doesn't work properly!

for i = [ 1 2 3 4 5]
    data = AccData.([char(Participants(i))])(:,[3,2,4]); % reorder data to 1 = V; 2 = ML; 3 = AP
    % Consistency of gait characteristics as determined from acceleration data collected at different trunk locations. Gait Posture 2014;40(1):187-92.
    [RealignedAcc, ~] = RealignSensorSignalHRAmp(data, FS);
    AccData.([char(Participants(i))]) = RealignedAcc;
    AccData.([char(Participants(i))]) = detrend(AccData.([char(Participants(i))])); % Detrend the data  
    clear data RealignedAcc
end

for i = [6 7 8] % Data collected without g, check if realignment works properly!
    AccData.([char(Participants(i))]) = AccData.([char(Participants(i))])(:,[3,2,4]); % reorder data to 1 = V; 2 = ML; 3 = AP
    disp(['Invert VT direction for participant: ' char(Participants(i))])
    AccData.([char(Participants(i))])(:,1) = AccData.([char(Participants(i))])(:,1) *-1;    % Invert VT signal
    AccData.([char(Participants(i))]) = detrend(AccData.([char(Participants(i))]));         % Detrend the data
    % NOT SURE WHAT TO DO WITH ML DIRECTION THOUGH!!
end

%
%% Determine epochs of 1 minute straight walking (for patients and controls seperately in function)
% Determine the location of the turns
% Relate the step frequency to the turns!

for i =  1 : length(Participants)
    inputData = AccData.([char(Participants(i))]);
    Distance = GeneralInformation.Distance(i);
    ApplyRealignment = false;
    [locsTurns,DataFilt,FootContacts] = DetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,Distance,ResultStruct);
    Turns.([char(Participants(i))]) = locsTurns;
    FilteredData.([char(Participants(i))]) = DataFilt;      % This data is realigned and filtered with a low-pass butterworth filter
    StepContacts.([char(Participants(i))]) = FootContacts;
    clear locsTurns DataFilt FootContacts
end

%
% Check if the Footcontacts and peaks in AP signal are correct for all participants
% Signals were corrected for horizontal tilt, and low-pass filtered using a Butterworth filter (2th order; cut-off frequency 0.6 Hz). Footcontacts were determined from the forward accelerations. I.e., the basic pattern of acceleration responds to a pattern predicted by an inverted pendulum; during single support, forward acceleration increaes when, after mid-stace, the body is falling forward and downward. ... '' Consequently, the forward acceleration of the lower trunk changes sign from positive to negative. The peak forward acceleration preceding the change of sign coincides with the instant of foot contact. The AP acceleration signal was low-pass filtered (4th order Butterworth filter at 2 Hz). A median filter was used to exclude outliers in the data due to turning points in the gait assessment. Foot contact data were used to calculate stride times, which were defined as the time interval between two ipsilateral foot contacts. Mean stride time,coefficient of variation (CV) of stride time, and stride frequency were calculated for each patient. Further parameters were calculated to assess variations throughout the gait cycle
% CHECK IF LOW PASS BUTTERWORTH FILTER IS CORRECT --> CUT-OFF FREQUENCY 20/(FS/2)

for k = 1 : length(Participants)
    StepContact = StepContacts.([char(Participants(k))]);
    dataAcc = AccData.([char(Participants(k))]);
    
        figure(k)
        subplot(311), plot(dataAcc(1:1000,1))    % VT
        subplot(312), plot(dataAcc(1:1000,2))    % ML
        subplot(313), plot(dataAcc(1:1000,3))    % AP
        set(findall(gcf,'type','axes'),'Xgrid', 'on','Xtick', StepContact, 'xlim',[1,1000])  % length(dataAcc) instead of 6000
        title(['Check footcontacts AP for participant: ', ([char(Participants(k))])]);
end

%
% Check PP1:

for k = 1
    StepContact = StepContacts.([char(Participants(k))]);
    dataAcc = AccData.([char(Participants(k))]);
    
        figure(k)
        subplot(311), plot(dataAcc(1:36000,1))    % VT
        subplot(312), plot(dataAcc(1:36000,2))    % ML
        subplot(313), plot(dataAcc(1:36000,3))    % AP
        set(findall(gcf,'type','axes'),'Xgrid', 'on','Xtick', StepContact, 'xlim',[1,36000])  % length(dataAcc) instead of 6000
        title(['Check footcontacts AP for participant: ', ([char(Participants(k))])]);
end

%
% Calculate the parameters for each minute seperately

Epochs = 1:60*FS:42000;                                     % Create epochs of 1 minute for 6MWT

for i = 1%:length(Participants)
    tic
    [~, EpochId] = histc(Turns.([char(Participants(i))])(:,1),Epochs);      % Classify the location of turns in the correct minute
    locsTurnsPerMinute = [Turns.([char(Participants(i))]) EpochId];
    LegLength = GeneralInformation.LegLength(i);                            % LegLength will be used for calculation of steplength 
    
    for IndexMinute = 1:Duration                            % Minute 1 till 6
        startPos = Epochs(IndexMinute);                     % Sample [1,6001,12001...]
        index = find(locsTurnsPerMinute(:,3)==IndexMinute)  % Find the turns within timeframe IndexMinute
        
        if isempty(index)                                   % If there is no turn detected, iterate only once (startPos - endPos)
            iteration = 1
        else                                                % Number of iterations depends on number of turns within timeframe IndexMinute
            iteration = size(index,1)
        end
        
        for k = 1:iteration
            
            if ~isempty(index)                              % There are turns detected in IndexMinute
                endPos = locsTurnsPerMinute(index(k),1)-1;  % End position is index of start turn (column 1)
            else
                endPos = Epochs(IndexMinute+1);             % If no turns detected, inputdata should be the complete minute. 
                disp('There was no turn detected within this minute, please verify (i.e. visual check or change parameters turn detection!');
            end
            
            if endPos-startPos < 13*FS                      % Otherwise error in StrideFrequencyFrom3dAcc (AccXYZ must be at least ten seconds long, and ApplyRemovesSteps), 
                disp('EndPos-StartPos < 13*FS = ERROR')
                continue                                    % Pass control to next iteration of for loop
            end
            
            inputData = AccData.([char(Participants(i))])(startPos:endPos,:);
            WindowLen = 10*FS;
            ApplyRealignment = false;                       % Data is already reorderd + detrend.
            ApplyRemoveSteps = true;                        % See Comment prof. Bert Otten! Remove 1 step before and after the turn!
            [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps); % Store main outcome parameters in struct
        
            %if k ==1                                        % only the first time
            Parameters = fieldnames(ResultStruct);
            NrParameters = length(Parameters);
            %end
    
            for j = 1:NrParameters                          % only works if for every bin we get the same outcomes (which is the case in this script)
            DataStraight.([char(Parameters(j))])(k) = ResultStruct.([char(Parameters(j))]);
            end
       
        if ~isempty(index)                                  % Start position of the next iteration should be at the end of the turn (column 2 of locsTurnsPerMinute)      
            startPos = locsTurnsPerMinute(index(k),2)+1;
        else
            continue                                        % If there is no turn detected within a minute, pass control to next iteration of for loop 
        end
        
        end
        
        % Calculate mean for inside straight for every minute (Indexminute) and store for each participant(i) seperately in struct
        for j = 1:NrParameters
            OutcomesAcc.([char(Participants(i))]).([char(Parameters(j))])(IndexMinute) = nanmean(DataStraight.([char(Parameters(j))]));
        end
        clear ResultStruct;
    end
    clear DataStraight;
    
    % Calculate mean for inside straight (6 minutes) for each participant and store in OutcomesAcc struct at column 7 
    for j = 1:NrParameters
            OutcomesAcc.([char(Participants(i))]).([char(Parameters(j))])(7) = nanmean(OutcomesAcc.([char(Participants(i))]).([char(Parameters(j))]));
    end
    toc
end
%uisave % Save workspace variables to MAT-file 

%
%% Create same matrix for controls as CLBP (detrended / reorder the data!)

% Create function!
load('ControlsData.mat');                % Acc Data of healthy participants. Detrended and in order; AP ML V
load('ControlsInformation_LD.mat');

%Concatenate both matrices (2x3 min.) vertically, and store in struct ControlsData
% ControlsData.P04 = cat(1,P04Walk1_IPod,P04Walk2_IPod);
% ControlsData.P12 = cat(1,P12Walk1_IPod,P12Walk2_IPod);
% ControlsData.P29 = cat(1,P29Walk1_IPod,P29Walk2_IPod);
% ControlsData.P48 = cat(1,P48Walk1_IPod,P48Walk2_IPod);
% ControlsData.P56 = cat(1,P56Walk1_IPod,P56Walk2_IPod);
% ControlsData.P60 = cat(1,P60Walk1_IPod,P60Walk2_IPod);
% ControlsData.P61 = cat(1,P61Walk1_IPod,P61Walk2_IPod);
% ControlsData.P63 = cat(1,P63Walk1_IPod,P63Walk2_IPod);

 Controls = fields(ControlsData);

%
%% Realign and Detrend the data (performed in previous stage)
% VT / AP direction inverted!
% 'Flipped over' around the ML direction, no need to invert ML direction at this stage.

for i = 1 : length(Controls)
    ControlsData.([char(Controls(i))]) = ControlsData.([char(Controls(i))])(:,[3,2,1]);        % Reorder data to 1 = V; 2 = ML; 3 = AP, same as in controls!
    ControlsData.([char(Controls(i))])(:,[1]) = ControlsData.([char(Controls(i))])(:,[1])*-1;  % Invert VT direction
    ControlsData.([char(Controls(i))])(:,[3]) = ControlsData.([char(Controls(i))])(:,[3])*-1;  % Invert AP direction
end
disp('PAY ATTENTION HERE: ONLY RUN THIS SECTION ONCE, OTHERWISE THE ORDER WILL BE REVERSED..')

%
% Determine the location of the turns for pwCLBP
% ControlsDetermineLocationTurnsFunc is adapted to controls data! (realignment performed)

ResultStruct = struct();                    % Empty struct

for i = 1 : length(Controls)
    plotit = true; 
    inputData = ControlsData.([char(Controls(i))]);
    ApplyRealignment = false;
    [locsTurns,DataFilt,FootContacts] = ControlsDetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,ResultStruct); % use DetermineLocationTurnsFunc
    Turns.([char(Controls(i))]) = locsTurns;
    FilteredData.([char(Controls(i))]) = DataFilt;      % This data is realigned
    StepContacts.([char(Controls(i))]) = FootContacts;
    clear locsTurns DataFilt FootContacts
end

%
% Check if the Footcontacts and peaks in AP signal are correct for all participants

for k = 1 : length(Controls)
    StepContact = StepContacts.([char(Controls(k))]);
    dataAcc = ControlsData.([char(Controls(k))]);
    
    figure()
    subplot(311), plot(dataAcc(1:3000,1))    % VT
    subplot(312), plot(dataAcc(1:3000,2))    % ML
    subplot(313), plot(dataAcc(1:3000,3))    % AP
    set(findall(gcf,'type','axes'),'Xgrid', 'on','Xtick', StepContact, 'xlim',[1,1000])  % length(dataAcc) instead of 6000
    title(['Check footcontacts AP for participant: ', ([char(Controls(k))])]);
end

%

Epochs = 1:60*FS:42000;                                     % Create epochs of 1 minute for 6MWT

for i = 1:length(Controls)
    tic
    [~, EpochId] = histc(Turns.([char(Controls(i))])(:,1),Epochs);       % Classify the location of turns in the correct minute
    locsTurnsPerMinute = [Turns .([char(Controls(i))]) EpochId];
    LegLength = ControlsInformation.LegLength(i);                        % LegLength will be used for calculation of steplength 
    
    for IndexMinute = 1:Duration                            % Minute 1 till 6
        startPos = Epochs(IndexMinute)                      % Sample [1,6001,12001...]
        index = find(locsTurnsPerMinute(:,3)==IndexMinute)  % Find the turns within timeframe IndexMinute
        
        if isempty(index)                                   % If there is no turn detected, iterate only once (startPos - endPos)
            iteration = 1
        else                                                % Number of iterations depends on number of turns within timeframe IndexMinute
            iteration = size(index,1)
        end
        
        for k = 1:iteration
            
            if ~isempty(index)                              % There are turns detected in IndexMinute
                endPos = locsTurnsPerMinute(index(k),1)-1   % End position is index of start turn (column 1)
            else
                endPos = Epochs(IndexMinute+1)              % If no turns detected, inputdata should be the complete minute. 
                disp('There was no turn detected within this minute, please verify (i.e. visual check or change parameters turn detection!');
            end
            
            if endPos-startPos < 12*FS                      % Otherwise error in StrideFrequencyFrom3dAcc (AccXYZ must be at least ten seconds long), 
                disp('EndPos-StartPos < 12*FS = ERROR')
                continue                                    % Pass control to next iteration of for loop
            end
            
            inputData = ControlsData.([char(Controls(i))])(startPos:endPos,:);
            WindowLen = 10*FS;
            ApplyRealignment = false;                       % Created a new apply realignment statement for the controls data! 
            ApplyRemoveSteps = true;                        % See Comment prof. Bert Otten! Remove 1 step before and after the turn!
            [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps); % Store main outcome parameters in struct
        
            %if k ==1                                        % only the first time
            Parameters = fieldnames(ResultStruct);
            NrParameters = length(Parameters);
            %end
    
            for j = 1:NrParameters                          % only works if for every bin we get the same outcomes (which is the case in this script)
            DataStraight.([char(Parameters(j))])(k) = ResultStruct.([char(Parameters(j))]);
            end
       
        if ~isempty(index)                                  % Start position of the next iteration should be at the end of the turn (column 2 of locsTurnsPerMinute)      
            startPos = locsTurnsPerMinute(index(k),2)+1;
        else
            continue                                        % If there is no turn detected within a minute, pass control to next iteration of for loop 
        end
        
        end
        
        % Calculate mean for inside straight for every minute (Indexminute) and store for each participant(i) seperately in struct
        for j = 1:NrParameters
            OutcomesAcc.([char(Controls(i))]).([char(Parameters(j))])(IndexMinute) = nanmean(DataStraight.([char(Parameters(j))]));
        end
        clear ResultStruct;
    end
    clear DataStraight;
    
    % Calculate mean for inside straight (6 minutes) for each participant and store in OutcomesAcc struct at column 7 
    for j = 1:NrParameters
            OutcomesAcc.([char(Controls(i))]).([char(Parameters(j))])(7) = nanmean(OutcomesAcc.([char(Controls(i))]).([char(Parameters(j))]));
    end
    toc
end

%uisave % Save workspace variables to MAT-file 

%
% I prefer Python for data visulisation. Thus outcomes will be exported as .xls file and further analysis performed in Python.
%
%