%% Gait Variability Analysis CLBP
%% Script created for MAP 2020-2021 LD
%% Adapted from Claudine Lamoth and Iris Hagoort
%% Version 2 Februari 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
%
%% 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

plotit = false;                             % If true, plot turn detection

%
%% Determine epochs of 1 minute straight walking (for patients and controls seperately in function?)
% Determine the location of the turns

for i = 1 : length(Participants)
    inputData = AccData.([char(Participants(i))]);
    Distance = GeneralInformation.Distance(i);
    ApplyRealignment = true;
    [locsTurns,DataFilt] = DetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,Distance);
    Turns.([char(Participants(i))]) = locsTurns;
    FilteredData.([char(Participants(i))]) = DataFilt;      % This data is realigned and filtered with a low-pass butterworth filter
    clear locsTurns
    clear DataFilt
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:6                                   % 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 < 10*FS                      % Otherwise error in StrideFrequencyFrom3dAcc (AccXYZ must be at least ten seconds long), 
                continue                                    % Pass control to next iteration of for loop
            end
            
            inputData = AccData.([char(Participants(i))])(startPos:endPos,:);
            WindowLen = size(inputData,1);
            ApplyRealignment = true; 
            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('HealthyWalk3min.mat');                % Acc Data of healthy participants. Detrended and in order; AP ML V

% 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.P35 = cat(1,P35Walk1_IPod,P35Walk2_IPod);
ControlsData.P48 = cat(1,P48Walk1_IPod,P48Walk2_IPod);
ControlsData.P56 = cat(1,P56Walk1_IPod,P56Walk2_IPod);
ControlsData.P61 = cat(1,P61Walk1_IPod,P61Walk2_IPod);
ControlsData.P63 = cat(1,P63Walk1_IPod,P63Walk2_IPod);

Controls = fields(ControlsData);

for i = 1:size(Controls,1)
    figure(i);
    plot(ControlsData.([char(Controls(i))])(:,3)); % only plot the VT direction
    title(([char(Controls(i))]));
end

%
% Plot distributions of the data, create scatter plot with marginal histograms

Controls = fields(ControlsData);

% Example with P04

for i = 1 %:length(Controls)
    inputData = ControlsData.([char(Controls(i))]);
    
    % Realign sensor data to VT-ML-AP frame + Detrend
    data = inputData(:,[3,2,1]);
    [RealignedAcc, ~] = RealignSensorSignalHRAmp(data, FS);
    dataAcc = RealignedAcc;
    % Filter data
    [B,A] = butter(2,3/(FS/2),'low'); % Filters data very strongly which is needed to determine turns correctly
    dataAcc_filt = filtfilt(B,A,dataAcc);
    
    figure(10)
    plot(dataAcc_filt(:,1))
    title('Same preprocessing as in main script -- PP04.')
end
    
%     % Band pass frequency range is between 0.5 and 30 Hz. 
%     f1 = 0.5; f2=30;
%     % Normalized pass frequency ranges (Wn1, Wn2):
%     Wn1=f1/(0.5*FS);      
%     Wn2=f2/(0.5*FS);      
%     % Second order (N=2) filter passes all signal components within           
%     % a frequency range of: [Wn1 Wn2]
%     N=2;
%     % Compute second order filter coefficients:
%     [b, a]    = butter(N, [Wn1, Wn2], 'bandpass'); 
%     filter_F  = filter(b, a, dataStepDetection);     
% end

% Create scatter plot with marginal histograms

NRS = reshape(GeneralInformation.NRS,[],1);

for i = 1:length(Participants)
    plot = Out
end

%