TEST remote repository

AutocorrRegSymmSteps.m

@@ -0,0 +1,7 @@
+% autocorrelation for step symmetrie
+function [c, lags] =  AutocorrRegSymmSteps(x) 
+[c,lags] = xcov(x,200,'unbiased');
+c = c/c(lags==0);
+c = c(lags>=0);
+lags = lags(lags>=0);
+

AutocorrStrides.m

@@ -0,0 +1,43 @@
+function  [ResultStruct] = AutocorrStrides(data,FS, StrideTimeRange,ResultStruct);
+
+%% Stride-times measures
+% Stride time and regularity from auto correlation (according to Moe-Nilssen and Helbostad, Estimation of gait cycle characteristics by trunk accelerometry. J Biomech, 2004. 37: 121-6.)
+RangeStart = round(FS*StrideTimeRange(1));
+RangeEnd = round(FS*StrideTimeRange(2));
+[AutocorrResult,Lags]=xcov(data,RangeEnd,'unbiased');
+AutocorrSum = sum(AutocorrResult(:,[1 5 9]),2); % This sum is independent of sensor re-orientation, as long as axes are kept orthogonal
+AutocorrResult2= [AutocorrResult(:,[1 5 9]),AutocorrSum];
+IXRange = (numel(Lags)-(RangeEnd-RangeStart)):numel(Lags);
+% check that autocorrelations are positive for any direction,
+
+AutocorrPlusTrans = AutocorrResult+AutocorrResult(:,[1 4 7 2 5 8 3 6 9]);
+IXRangeNew = IXRange( ...
+    AutocorrPlusTrans(IXRange,1) > 0 ...
+    & prod(AutocorrPlusTrans(IXRange,[1 5]),2) > prod(AutocorrPlusTrans(IXRange,[2 4]),2) ...
+    & prod(AutocorrPlusTrans(IXRange,[1 5 9]),2) + prod(AutocorrPlusTrans(IXRange,[2 6 7]),2) + prod(AutocorrPlusTrans(IXRange,[3 4 8]),2) ...
+    > prod(AutocorrPlusTrans(IXRange,[1 6 8]),2) + prod(AutocorrPlusTrans(IXRange,[2 4 9]),2) + prod(AutocorrPlusTrans(IXRange,[3 5 7]),2) ...
+    );
+if isempty(IXRangeNew)
+    StrideTimeSamples = Lags(IXRange(AutocorrSum(IXRange)==max(AutocorrSum(IXRange)))); % to be used in other estimations
+    StrideTimeSeconds = nan;
+    ResultStruct.StrideTimeSamples = StrideTimeSamples;
+    ResultStruct.StrideTimeSeconds = StrideTimeSeconds;
+    
+else
+    StrideTimeSamples = Lags(IXRangeNew(AutocorrSum(IXRangeNew)==max(AutocorrSum(IXRangeNew))));
+    StrideRegularity = AutocorrResult2(Lags==StrideTimeSamples,:)./AutocorrResult2(Lags==0,:); % Moe-Nilssen&Helbostatt,2004
+    RelativeStrideVariability = 1-StrideRegularity;
+    StrideTimeSeconds = StrideTimeSamples/FS;
+    
+    ResultStruct.StrideRegularity_V = StrideRegularity(1);
+    ResultStruct.StrideRegularity_ML = StrideRegularity(2);
+    ResultStruct.StrideRegularity_AP = StrideRegularity(3);
+    ResultStruct.StrideRegularity_All = StrideRegularity(4);
+    ResultStruct.RelativeStrideVariability_V = RelativeStrideVariability(1);
+    ResultStruct.RelativeStrideVariability_ML = RelativeStrideVariability(2);
+    ResultStruct.RelativeStrideVariability_AP = RelativeStrideVariability(3);
+    ResultStruct.RelativeStrideVariability_All = RelativeStrideVariability(4);
+    ResultStruct.StrideTimeSamples = StrideTimeSamples;
+    ResultStruct.StrideTimeSeconds = StrideTimeSeconds;
+    
+end

CalcMaxLyapConvGait.m

@@ -0,0 +1,118 @@
+function [L_Estimate,ExtraArgsOut] = CalcMaxLyapConvGait(ThisTimeSeries,FS,ExtraArgsIn)
+if nargin > 2
+    if isfield(ExtraArgsIn,'J')
+        J=ExtraArgsIn.J;
+    end
+    if isfield(ExtraArgsIn,'m')
+        m=ExtraArgsIn.m;
+    end
+    if isfield(ExtraArgsIn,'FitWinLen')
+        FitWinLen=ExtraArgsIn.FitWinLen;
+    end
+end
+
+%% Initialize output args
+L_Estimate=nan;ExtraArgsOut.Divergence=nan;ExtraArgsOut.J=nan;ExtraArgsOut.m=nan;ExtraArgsOut.FitWinLen=nan;
+
+%% Some checks 
+% predefined J and m should not be NaN or Inf
+if (exist('J','var') && ~isempty(J) && ~isfinite(J)) || (exist('m','var') && ~isempty(m) && ~isfinite(m))
+    warning('Predefined J and m cannot be NaN or Inf');
+    return;
+end
+% multidimensional time series need predefined J and m
+if size(ThisTimeSeries,2) > 1 && (~exist('J','var') || ~exist('m','var') || isempty(J) || isempty(m)) 
+    warning('Multidimensional time series needs predefined J and m, can''t determine Lyapunov');
+    return;
+end
+%Check that there are no NaN or Inf values in the TimeSeries
+if any(~isfinite(ThisTimeSeries(:)))
+    warning('Time series contains NaN or Inf, can''t determine Lyapunov');
+    return;
+end
+%Check that there is variation in the TimeSeries
+if ~(nanstd(ThisTimeSeries) > 0)
+    warning('Time series is constant, can''t determine Lyapunov');
+    return;
+end
+
+%% Determine FitWinLen (=cycle time) of ThisTimeSeries
+if ~exist('FitWinLen','var') || isempty(FitWinLen)
+    if size(ThisTimeSeries,2)>1
+        for dim=1:size(ThisTimeSeries,2),
+            [Pd(:,dim),F] = pwelch(detrend(ThisTimeSeries(:,dim)),[],[],[],FS);
+        end
+        P = sum(Pd,2);
+    else
+        [P,F] = pwelch(detrend(ThisTimeSeries),[],[],[],FS);
+    end
+    MeanF = sum(P.*F)./sum(P);
+    CycleTime = 1/MeanF;
+    FitWinLen = round(CycleTime*FS);
+else
+    CycleTime = FitWinLen/FS;
+end
+ExtraArgsOut.FitWinLen=FitWinLen;
+
+%% Determine J
+if ~exist('J','var') || isempty(J)
+    % Calculate mutual information and take first local minimum Tau as J
+    bV = min(40,floor(sqrt(size(ThisTimeSeries,1))));
+    tauVmax = FitWinLen;
+    [mutMPro,cummutMPro,minmuttauVPro] = MutualInformationHisPro(ThisTimeSeries,(0:tauVmax),bV,1); % (xV,tauV,bV,flag)
+    if isnan(minmuttauVPro)
+        display(mutMPro);
+        warning('minmuttauVPro is NaN. Consider increasing tauVmax.');
+        return;
+    end
+    J=minmuttauVPro;
+end
+ExtraArgsOut.J=J;
+
+%% Determine m
+if ~exist('m','var') || isempty(m)
+    escape = 10;
+    max_m = 20;
+    max_fnnM = 0.02;
+    mV = 0;
+    fnnM = 1;
+    for mV = 2:max_m % for m=1, FalseNearestNeighbors is slow and lets matlab close if N>500000
+        fnnM = FalseNearestNeighborsSR(ThisTimeSeries,J,mV,escape,FS); % (xV,tauV,mV,escape,theiler)
+        if fnnM <= max_fnnM || isnan(fnnM)
+            break
+        end
+    end
+    if fnnM <= max_fnnM
+        m = mV;
+    else
+        warning('Too many false nearest neighbours');
+        return;
+    end
+end
+ExtraArgsOut.m=m;
+
+%% Create state space based upon J and m
+N_ss = size(ThisTimeSeries,1)-(m-1)*J;
+StateSpace=nan(N_ss,m*size(ThisTimeSeries,2)); 
+for dim=1:size(ThisTimeSeries,2),
+    for delay=1:m,
+        StateSpace(:,(dim-1)*m+delay)=ThisTimeSeries((1:N_ss)'+(delay-1)*J,dim);
+    end
+end
+
+%% Parameters for Lyapunov
+WindowLen = floor(min(N_ss/5,10*FitWinLen));
+if WindowLen < FitWinLen
+        warning('Not enough samples for Lyapunov estimation');
+        return;
+end    
+WindowLenSec=WindowLen/FS;
+
+%% Calculate divergence
+Divergence=div_calc(StateSpace,WindowLenSec,FS,CycleTime,0);
+ExtraArgsOut.Divergence=Divergence;
+
+%% Calculate slope of first FitWinLen samples of divergence curve
+p = polyfit((1:FitWinLen)/FS,Divergence(1:FitWinLen),1);
+L_Estimate = p(1);
+

CalcMaxLyapWolfFixedEvolv.m

@@ -0,0 +1,134 @@
+function [L_Estimate,ExtraArgsOut] = CalcMaxLyapWolfFixedEvolv(ThisTimeSeries,FS,ExtraArgsIn)
+
+%% Description
+% This function calculates the maximum Lyapunov exponent from a time 
+% series, based on the method described by Wolf et al. in 
+%    Wolf, A., et al., Determining Lyapunov exponents from a time series. 
+%    Physica D: 8 Nonlinear Phenomena, 1985. 16(3): p. 285-317.
+% 
+% Input:
+%   ThisTimeSeries: a vector or matrix with the time series
+%   FS: sample frequency of the ThisTimeSeries
+%   ExtraArgsIn: a struct containing optional input arguments 
+%       J (embedding delay)
+%       m (embedding dimension)
+% Output:
+%   L_Estimate: The Lyapunov estimate
+%   ExtraArgsOut: a struct containing the additional output arguments
+%       J (embedding delay)
+%       m (embedding dimension)
+
+%% Copyright
+%     COPYRIGHT (c) 2012 Sietse Rispens, VU University Amsterdam
+%
+%     This program is free software: you can redistribute it and/or modify
+%     it under the terms of the GNU General Public License as published by
+%     the Free Software Foundation, either version 3 of the License, or
+%     (at your option) any later version.
+% 
+%     This program is distributed in the hope that it will be useful,
+%     but WITHOUT ANY WARRANTY; without even the implied warranty of
+%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%     GNU General Public License for more details.
+% 
+%     You should have received a copy of the GNU General Public License
+%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+%% Author
+%  Sietse Rispens
+
+%% History
+%  April 2012, initial version of CalcMaxLyapWolf
+%  23 October 2012, use fixed evolve time instead of adaptable
+
+if nargin > 2
+    if isfield(ExtraArgsIn,'J')
+        J=ExtraArgsIn.J;
+    end
+    if isfield(ExtraArgsIn,'m')
+        m=ExtraArgsIn.m;
+    end
+end
+
+%% Initialize output args
+L_Estimate=nan;ExtraArgsOut.J=nan;ExtraArgsOut.m=nan;
+
+%% Some checks 
+% predefined J and m should not be NaN or Inf
+if (exist('J','var') && ~isempty(J) && ~isfinite(J)) || (exist('m','var') && ~isempty(m) && ~isfinite(m))
+    warning('Predefined J and m cannot be NaN or Inf');
+    return;
+end
+% multidimensional time series need predefined J and m
+if size(ThisTimeSeries,2) > 1 && (~exist('J','var') || ~exist('m','var') || isempty(J) || isempty(m)) 
+    warning('Multidimensional time series needs predefined J and m, can''t determine Lyapunov');
+    return;
+end
+%Check that there are no NaN or Inf values in the TimeSeries
+if any(~isfinite(ThisTimeSeries(:)))
+    warning('Time series contains NaN or Inf, can''t determine Lyapunov');
+    return;
+end
+%Check that there is variation in the TimeSeries
+if ~(nanstd(ThisTimeSeries) > 0)
+    warning('Time series is constant, can''t determine Lyapunov');
+    return;
+end
+
+%% Determine J
+if ~exist('J','var') || isempty(J)
+    % Calculate mutual information and take first local minimum Tau as J
+    bV = min(40,floor(sqrt(size(ThisTimeSeries,1))));
+    tauVmax = 70;
+    [mutMPro,cummutMPro,minmuttauVPro] = MutualInformationHisPro(ThisTimeSeries,(0:tauVmax),bV,1); % (xV,tauV,bV,flag)
+    if isnan(minmuttauVPro)
+        display(mutMPro);
+        warning('minmuttauVPro is NaN. Consider increasing tauVmax.');
+        return;
+    end
+    J=minmuttauVPro;
+end
+ExtraArgsOut.J=J;
+
+%% Determine m
+if ~exist('m','var') || isempty(m)
+    escape = 10;
+    max_m = 20;
+    max_fnnM = 0.02;
+    mV = 0;
+    fnnM = 1;
+    for mV = 2:max_m % for m=1, FalseNearestNeighbors is slow and lets matlab close if N>500000
+        fnnM = FalseNearestNeighborsSR(ThisTimeSeries,J,mV,escape,FS); % (xV,tauV,mV,escape,theiler)
+        if fnnM <= max_fnnM || isnan(fnnM)
+            break
+        end
+    end
+    if fnnM <= max_fnnM
+        m = mV;
+    else
+        warning('Too many false nearest neighbours');
+        return;
+    end
+end
+ExtraArgsOut.m=m;
+
+%% Create state space based upon J and m
+N_ss = size(ThisTimeSeries,1)-(m-1)*J;
+StateSpace=nan(N_ss,m*size(ThisTimeSeries,2)); 
+for dim=1:size(ThisTimeSeries,2),
+    for delay=1:m,
+        StateSpace(:,(dim-1)*m+delay)=ThisTimeSeries((1:N_ss)'+(delay-1)*J,dim);
+    end
+end
+
+%% Parameters for Lyapunov estimation
+CriticalLen=J*m;
+max_dist = sqrt(sum(std(StateSpace).^2))/10;
+max_dist_mult = 5;
+min_dist = max_dist/2;
+max_theta = 0.3;
+evolv = J;
+
+%% Calculate Lambda
+[L_Estimate]=div_wolf_fixed_evolv(StateSpace, FS, min_dist, max_dist, max_dist_mult, max_theta, CriticalLen, evolv);
+

CalculateNonLinearParametersFunc.m

@@ -0,0 +1,48 @@
+function [ResultStruct] = CalculateNonLinearParametersFunc(ResultStruct,dataAccCut,WindowLen,FS,Lyap_m,Lyap_FitWinLen,Sen_m,Sen_r)
+
+%% Calculation non-linear parameters;
+
+% cut into windows of size WindowLen
+N_Windows = floor(size(dataAccCut,1)/WindowLen);
+N_SkipBegin = ceil((size(dataAccCut,1)-N_Windows*WindowLen)/2);
+LyapunovWolf = nan(N_Windows,3);
+LyapunovRosen = nan(N_Windows,3);
+SE= nan(N_Windows,3);
+
+for WinNr = 1:N_Windows;
+    AccWin = dataAccCut(N_SkipBegin+(WinNr-1)*WindowLen+(1:WindowLen),:);
+    for j=1:3
+        [LyapunovWolf(WinNr,j),~] = CalcMaxLyapWolfFixedEvolv(AccWin(:,j),FS,struct('m',Lyap_m));
+        [LyapunovRosen(WinNr,j),outpo] = CalcMaxLyapConvGait(AccWin(:,j),FS,struct('m',Lyap_m,'FitWinLen',Lyap_FitWinLen));
+        [SE(WinNr,j)] = funcSampleEntropy(AccWin(:,j), Sen_m, Sen_r);
+        % no correction for FS; SE does increase with higher FS but effect is considered negligible as range is small (98-104HZ). Might consider updating r to account for larger ranges.
+    end
+end
+
+LyapunovWolf = nanmean(LyapunovWolf,1);
+LyapunovRosen = nanmean(LyapunovRosen,1);
+SampleEntropy = nanmean(SE,1);
+
+ResultStruct.LyapunovWolf_V = LyapunovWolf(1);
+ResultStruct.LyapunovWolf_ML = LyapunovWolf(2);
+ResultStruct.LyapunovWolf_AP = LyapunovWolf(3);
+ResultStruct.LyapunovRosen_V = LyapunovRosen(1);
+ResultStruct.LyapunovRosen_ML = LyapunovRosen(2);
+ResultStruct.LyapunovRosen_AP = LyapunovRosen(3);
+ResultStruct.SampleEntropy_V = SampleEntropy(1);
+ResultStruct.SampleEntropy_ML = SampleEntropy(2);
+ResultStruct.SampleEntropy_AP = SampleEntropy(3);
+
+if isfield(ResultStruct,'StrideFrequency')
+    LyapunovPerStrideWolf = LyapunovWolf/ResultStruct.StrideFrequency;
+    LyapunovPerStrideRosen = LyapunovRosen/ResultStruct.StrideFrequency;
+end
+
+ResultStruct.LyapunovPerStrideWolf_V = LyapunovPerStrideWolf(1);
+ResultStruct.LyapunovPerStrideWolf_ML = LyapunovPerStrideWolf(2);
+ResultStruct.LyapunovPerStrideWolf_AP = LyapunovPerStrideWolf(3);
+ResultStruct.LyapunovPerStrideRosen_V = LyapunovPerStrideRosen(1);
+ResultStruct.LyapunovPerStrideRosen_ML = LyapunovPerStrideRosen(2);
+ResultStruct.LyapunovPerStrideRosen_AP = LyapunovPerStrideRosen(3);
+
+end 

CalculateStrideParametersFunc.m

@@ -0,0 +1,48 @@
+function [ResultStruct] = CalculateStrideParametersFunc(dataAccCut_filt,FS,ApplyRemoveSteps,dataAcc,dataAcc_filt,StrideTimeRange)
+
+%% Calculate stride parameters
+ResultStruct = struct; % create empty struct
+
+% Run function AutoCorrStrides, Outcomeparameters: StrideRegularity,RelativeStrideVariability,StrideTimeSamples,StrideTime
+[ResultStruct] = AutocorrStrides(dataAccCut_filt,FS, StrideTimeRange,ResultStruct);
+StrideTimeSamples = ResultStruct.StrideTimeSamples; % needed as input for other functions
+
+% Calculate Step symmetry --> method 1
+ij = 1;
+dirSymm = [1,3]; % Gait Synmmetry is only informative in AP/V direction: See Tura A, Raggi M, Rocchi L, Cutti AG, Chiari L: Gait symmetry and regularity in transfemoral amputees assessed by trunk accelerations. J Neuroeng Rehabil 2010, 7:4.
+
+for jk=1:length(dirSymm)
+    [C, lags] = AutocorrRegSymmSteps(dataAccCut_filt(:,dirSymm(jk)));
+    [Ad,p] = findpeaks(C,'MinPeakProminence',0.2, 'MinPeakHeight', 0.2);
+    
+    if size(Ad,1) > 1
+        Ad1 = Ad(1);
+        Ad2 = Ad(2);
+        GaitSymm(:,ij) = abs((Ad1-Ad2)/mean([Ad1+Ad2]))*100;
+    else
+        GaitSymm(:,ij) = NaN;
+    end
+    ij = ij +1;
+end
+% Save outcome in struct;
+ResultStruct.GaitSymm_V = GaitSymm(1);
+ResultStruct.GaitSymm_AP = GaitSymm(2);
+
+% Calculate Step symmetry --> method 2
+[PksAndLocsCorrected] = StepcountFunc(dataAccCut_filt,StrideTimeSamples,FS);
+LocsSteps = PksAndLocsCorrected(2:2:end,2);
+
+if rem(size(LocsSteps,1),2) == 0; % is number of steps is even
+    LocsSteps2 = LocsSteps(1:2:end);
+else
+    LocsSteps2 = LocsSteps(3:2:end);
+end
+
+LocsSteps1 = LocsSteps(2:2:end);
+DiffLocs2 = diff(LocsSteps2);
+DiffLocs1 = diff(LocsSteps1);
+StepTime2 = DiffLocs2(1:end-1)/FS; % leave last one out because it is higher
+StepTime1 = DiffLocs1(1:end-1)/FS;
+SI = abs((2*(StepTime2-StepTime1))./(StepTime2+StepTime1))*100;
+ResultStruct.GaitSymmIndex = nanmean(SI);
+end 

FilterandRealignFunc.m

@@ -0,0 +1,18 @@
+function [dataAcc, dataAcc_filt] = FilterandRealignFunc(inputData,FS,ApplyRealignment)
+
+%% Filter and Realign Accdata
+
+% Apply Realignment & Filter data
+
+if ApplyRealignment % apply relignment as described in Rispens S, Pijnappels M, van Schooten K, Beek PJ, Daffertshofer A, van Die?n JH (2014).
+    data = inputData(:, [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);
+    dataAcc = RealignedAcc;
+    [B,A] = butter(2,20/(FS/2),'low');
+    dataAcc_filt = filtfilt(B,A,dataAcc);
+else % we asume tat data is already reorderd to 1 = V; 2= ML, 3 = AP in an earlier stage;
+    [B,A] = butter(2,20/(FS/2),'low');
+    dataAcc = inputData;
+    dataAcc_filt = filtfilt(B,A,dataAcc);
+end

GaitOutcomesTrunkAccFuncIH.m

@@ -0,0 +1,210 @@
+function [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps)
+
+% DESCRIPTON: Trunk analysis of Iphone data without the need for step detection
+% CL Nov 2019
+% Adapted IH feb-april 2020
+
+% koloms data of smartphone
+% 1st column is time data;
+% 2nd column is X, medio-lateral: + left, - right
+% 3rd column is Y, vertical: + downwards, - upwards
+% 4th column is Z, anterior- posterior : + forwards, - backwards
+
+%% Input Trunk accelerations during locomotion in VT, ML, AP direction
+% InputData: Acceleration signal with time and accelerations in VT,ML and
+% AP direction.
+% FS: sample frequency of the Accdata
+% LegLength: length of the leg of the participant in m;
+
+
+%% Output
+% ResultStruct: structure coninting all outcome measured calculated
+% Spectral parameters, spatiotemporal gait parameters, non-linear
+% parameters
+% fields and subfields: include the multiple measurements of a subject
+
+%% Literature
+% Richman & Moorman, 2000; [ sample entropy]
+% Bisi & Stagni Gait & Posture 2016,  47 (6) 37-42
+% Kavagnah et al., Eur J Appl Physiol 2005 94: 468?475; Human Movement Science 24(2005) 574?587 [ synchrony]
+% Moe-Nilsen J Biomech 2004 37, 121-126 [ autorcorrelation step regularity and symmetry
+% Kobsar et al. Gait & Posture 2014 39, 553?557 [ synchrony ]
+% Rispen et al; Gait & Posture 2014, 40, 187 - 192 [realignment axes]
+% Zijlstra & HofGait & Posture 2003  18,2, 1-10 [spatiotemporal gait variables]
+% Lamoth et al, 2002  [index of harmonicity]
+% Costa et al. 2003 Physica A 330 (2003) 53–60 [ multiscale entropy]
+% Cignetti F, Decker LM, Stergiou N. Ann Biomed Eng. 2012
+% May;40(5):1122-30. doi: 10.1007/s10439-011-0474-3. Epub 2011 Nov 25. [
+% Wofl vs. Rosenstein Lyapunov]
+
+
+%% Settings
+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();
+
+%% Filter and Realign Accdata
+
+% Apply Realignment & Filter data
+
+if ApplyRealignment % apply relignment as described in Rispens S, Pijnappels M, van Schooten K, Beek PJ, Daffertshofer A, van Die?n JH (2014).
+    data = inputData(:, [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);
+    dataAcc = RealignedAcc;
+    [B,A] = butter(2,20/(FS/2),'low');
+    dataAcc_filt = filtfilt(B,A,dataAcc);
+else % we asume tat data is already reorderd to 1 = V; 2= ML, 3 = AP in an earlier stage;
+    [B,A] = butter(2,20/(FS/2),'low');
+    dataAcc = inputData;
+    dataAcc_filt = filtfilt(B,A,dataAcc);
+end
+
+
+%% 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
+
+if ApplyRemoveSteps
+    
+    % In order to run the step detection script we first need to run an autocorrelation function;
+    [ResultStruct] = AutocorrStrides(dataAcc_filt,FS, StrideTimeRange,ResultStruct);
+    
+    % StrideTimeSamples is needed as an input for the stepcountFunc;
+    StrideTimeSamples = ResultStruct.StrideTimeSamples;
+    
+    % Calculate the number of steps;
+    [PksAndLocsCorrected] = StepcountFunc(dataAcc_filt,StrideTimeSamples,FS);
+    % This function selects steps based on negative and positive values.
+    % However to determine the steps correctly we only need one of these;
+    LocsSteps = PksAndLocsCorrected(1:2:end,2);
+    
+    %% Cut data & remove currents results
+    % Remove 20 steps in the beginning and end of data
+    dataAccCut = dataAcc(LocsSteps(31):LocsSteps(end-30),:);
+    dataAccCut_filt = dataAcc_filt(LocsSteps(31):LocsSteps(end-30),:);
+    
+    % Clear currently saved results from Autocorrelation Analysis
+    
+    clear ResultStruct;
+    clear PksAndLocsCorrected;
+    clear LocsSteps;
+    
+else;
+    dataAccCut = dataAcc;
+    dataAccCut_filt = dataAcc_filt;
+end
+
+%% Calculate stride parameters
+ResultStruct = struct; % create empty struct
+
+% Run function AutoCorrStrides, Outcomeparameters: StrideRegularity,RelativeStrideVariability,StrideTimeSamples,StrideTime
+[ResultStruct] = AutocorrStrides(dataAccCut_filt,FS, StrideTimeRange,ResultStruct);
+StrideTimeSamples = ResultStruct.StrideTimeSamples; % needed as input for other functions
+
+% Calculate Step symmetry --> method 1
+ij = 1;
+dirSymm = [1,3]; % Gait Synmmetry is only informative in AP/V direction: See Tura A, Raggi M, Rocchi L, Cutti AG, Chiari L: Gait symmetry and regularity in transfemoral amputees assessed by trunk accelerations. J Neuroeng Rehabil 2010, 7:4.
+
+for jk=1:length(dirSymm)
+    [C, lags] = AutocorrRegSymmSteps(dataAccCut_filt(:,dirSymm(jk)));
+    [Ad,p] = findpeaks(C,'MinPeakProminence',0.2, 'MinPeakHeight', 0.2);
+    
+    if size(Ad,1) > 1
+        Ad1 = Ad(1);
+        Ad2 = Ad(2);
+        GaitSymm(:,ij) = abs((Ad1-Ad2)/mean([Ad1+Ad2]))*100;
+    else
+        GaitSymm(:,ij) = NaN;
+    end
+    ij = ij +1;
+end
+% Save outcome in struct;
+ResultStruct.GaitSymm_V = GaitSymm(1);
+ResultStruct.GaitSymm_AP = GaitSymm(2);
+
+% Calculate Step symmetry --> method 2
+[PksAndLocsCorrected] = StepcountFunc(dataAccCut_filt,StrideTimeSamples,FS);
+LocsSteps = PksAndLocsCorrected(2:2:end,2);
+
+if rem(size(LocsSteps,1),2) == 0; % is number of steps is even
+    LocsSteps2 = LocsSteps(1:2:end);
+else
+    LocsSteps2 = LocsSteps(3:2:end);
+end
+
+LocsSteps1 = LocsSteps(2:2:end);
+DiffLocs2 = diff(LocsSteps2);
+DiffLocs1 = diff(LocsSteps1);
+StepTime2 = DiffLocs2(1:end-1)/FS; % leave last one out because it is higher
+StepTime1 = DiffLocs1(1:end-1)/FS;
+SI = abs((2*(StepTime2-StepTime1))./(StepTime2+StepTime1))*100;
+ResultStruct.GaitSymmIndex = nanmean(SI);
+
+%% Calculate spatiotemporal stride parameters
+
+% Measures from height variation by double integration of VT accelerations and high-pass filtering
+[ResultStruct] = SpatioTemporalGaitParameters(dataAccCut_filt,StrideTimeSamples,ApplyRealignment,LegLength,FS,IgnoreMinMaxStrides,ResultStruct);
+
+%% Measures derived from spectral analysis
+
+AccVectorLen = sqrt(sum(dataAccCut_filt(:,1:3).^2,2));
+[ResultStruct] = SpectralAnalysisGaitfunc(dataAccCut_filt,WindowLen,FS,N_Harm,LowFrequentPowerThresholds,AccVectorLen,ResultStruct);
+
+
+%% Calculation non-linear parameters;
+
+% cut into windows of size WindowLen
+N_Windows = floor(size(dataAccCut,1)/WindowLen);
+N_SkipBegin = ceil((size(dataAccCut,1)-N_Windows*WindowLen)/2);
+LyapunovWolf = nan(N_Windows,3);
+LyapunovRosen = nan(N_Windows,3);
+SE= nan(N_Windows,3);
+
+for WinNr = 1:N_Windows;
+    AccWin = dataAccCut(N_SkipBegin+(WinNr-1)*WindowLen+(1:WindowLen),:);
+    for j=1:3
+        [LyapunovWolf(WinNr,j),~] = CalcMaxLyapWolfFixedEvolv(AccWin(:,j),FS,struct('m',Lyap_m));
+        [LyapunovRosen(WinNr,j),outpo] = CalcMaxLyapConvGait(AccWin(:,j),FS,struct('m',Lyap_m,'FitWinLen',Lyap_FitWinLen));
+        [SE(WinNr,j)] = funcSampleEntropy(AccWin(:,j), Sen_m, Sen_r);
+        % no correction for FS; SE does increase with higher FS but effect is considered negligible as range is small (98-104HZ). Might consider updating r to account for larger ranges.
+    end
+end
+
+LyapunovWolf = nanmean(LyapunovWolf,1);
+LyapunovRosen = nanmean(LyapunovRosen,1);
+SampleEntropy = nanmean(SE,1);
+
+ResultStruct.LyapunovWolf_V = LyapunovWolf(1);
+ResultStruct.LyapunovWolf_ML = LyapunovWolf(2);
+ResultStruct.LyapunovWolf_AP = LyapunovWolf(3);
+ResultStruct.LyapunovRosen_V = LyapunovRosen(1);
+ResultStruct.LyapunovRosen_ML = LyapunovRosen(2);
+ResultStruct.LyapunovRosen_AP = LyapunovRosen(3);
+ResultStruct.SampleEntropy_V = SampleEntropy(1);
+ResultStruct.SampleEntropy_ML = SampleEntropy(2);
+ResultStruct.SampleEntropy_AP = SampleEntropy(3);
+
+if isfield(ResultStruct,'StrideFrequency')
+    LyapunovPerStrideWolf = LyapunovWolf/ResultStruct.StrideFrequency;
+    LyapunovPerStrideRosen = LyapunovRosen/ResultStruct.StrideFrequency;
+end
+
+ResultStruct.LyapunovPerStrideWolf_V = LyapunovPerStrideWolf(1);
+ResultStruct.LyapunovPerStrideWolf_ML = LyapunovPerStrideWolf(2);
+ResultStruct.LyapunovPerStrideWolf_AP = LyapunovPerStrideWolf(3);
+ResultStruct.LyapunovPerStrideRosen_V = LyapunovPerStrideRosen(1);
+ResultStruct.LyapunovPerStrideRosen_ML = LyapunovPerStrideRosen(2);
+ResultStruct.LyapunovPerStrideRosen_AP = LyapunovPerStrideRosen(3);
+
+end

GaitVariabilityAnalysisIH.m

@@ -0,0 +1,74 @@
+% Gait Variability Analysis 
+% Script created for BAP students 2020
+% Iris Hagoort
+% April 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
+load('Phyphoxdata.mat'); % loads accelerometer data, is stored in struct with name AccData
+load('ExcelInfo.mat');
+Participants = fields(AccData);
+
+%% Settings
+FS = 100; % sample frequency
+LegLengths = excel.data.GeneralInformation(:,5); % leglength info is in 5th column
+LegLengthsM = LegLengths/100; % convert to m
+
+%% Calculate parameters;
+for i = 1:  length(Participants);
+tic;   
+    LegLength = LegLengthsM(i);
+    WalkingConditions = fields(AccData.([char(Participants(i))]));
+    
+    for j = 1: length(WalkingConditions);
+        
+        if strcmp(char(WalkingConditions(j)),'Treadmill')
+            
+            SubConditions = fieldnames(AccData.([char(Participants(i))]).([char(WalkingConditions(j))]));
+            
+            for k = 1: length(SubConditions);
+               inputData = AccData.([char(Participants(i))]).([char(WalkingConditions(j))]).([char(SubConditions(k))]);
+               WindowLength = FS*10;
+               ApplyRealignment = true;
+               ApplyRemoveSteps = true;
+               [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps);
+                OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]).([char(SubConditions(k))]) = ResultStruct;
+            end
+            
+        elseif strcmp(char(WalkingConditions(j)),'Balance') || strcmp(char(WalkingConditions(j)),'TwoMWT') 
+              disp('Files are not used for current analysis');
+            
+        elseif strcmp(char(WalkingConditions(j)),'InsideStraight')
+            inputData = AccData.([char(Participants(i))]).([char(WalkingConditions(j))]);
+            ApplyRealignment = true;
+            ApplyRemoveSteps = false; % don't remove steps for the straight conditions
+            % function specific for the walking conditions containing a lot
+            % of turns
+            [ResultStruct] = GaitVariabilityAnalysisIH_WithoutTurns(inputData,FS,LegLength,ApplyRealignment,ApplyRemoveSteps);
+            OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]) = ResultStruct;
+            
+        else
+            inputData = AccData.([char(Participants(i))]).([char(WalkingConditions(j))]);
+            ApplyRealignment = true;
+            ApplyRemoveSteps = true;
+            WindowLen = FS*10;
+            [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps)
+            OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]) = ResultStruct;
+        end
+    end
+    
+ toc;  
+end
+
+% Save struct as .mat file
+save('GaitVarOutcomes30pril.mat', 'OutcomesAcc');
+
+
+

GaitVariabilityAnalysisIH_WithoutTurns.m

@@ -0,0 +1,108 @@
+function [ResultStruct] = GaitVariabilityAnalysisIH_WithoutTurns(inputData,FS,LegLength,ApplyRealignment,ApplyRemoveSteps);
+
+
+% SCRIPT FOR ANAlysis straight parts
+% NOG GOEDE BESCHRIJVING TOEVOEGEN.
+
+%% Realign data
+data = inputData(:, [3,2,4]); % reorder data to  1 = V; 2= ML, 3 = AP
+
+%Realign sensor data to VT-ML-AP frame
+if ApplyRealignment % apply relignment as described in Rispens S, Pijnappels M, van Schooten K, Beek PJ, Daffertshofer A, van Die?n JH (2014).
+    % 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);
+    dataAcc = RealignedAcc;
+end
+
+%% Filter data strongly & Determine location of steps
+
+% Filter  data
+[B,A] = butter(2,3/(FS/2),'low'); % Filters data very strongly which is needed to determine turns correctly
+dataStepDetection = filtfilt(B,A,dataAcc);
+
+% Determine steps;
+
+%%%%%%% HIER MISSCHIEN ALTERNATIEF VOOR VAN RISPENS %%%%%%%%%%%%%
+
+% Explanation of method: https://nl.mathworks.com/help/supportpkg/beagleboneblue/ref/counting-steps-using-beagleboneblue-hardware-example.html
+% From website: To convert the XYZ acceleration vectors at each point in time into scalar values,
+% calculate the magnitude of each vector. This way, you can detect large changes in overall acceleration,
+% such as steps taken while walking, regardless of device orientation.
+
+magfilt = sqrt(sum((dataStepDetection(:,1).^2) + (dataStepDetection(:,2).^2) + (dataStepDetection(:,3).^2), 2));
+magNoGfilt = magfilt - mean(magfilt);
+minPeakHeight2 = std(magNoGfilt);
+[pks, locs] = findpeaks(magNoGfilt, 'MINPEAKHEIGHT', minPeakHeight2); % for step detection
+numStepsOption2_filt = numel(pks); % counts number of steps;
+
+%%%%%%%%%%%%%%%%%%%%%%%% TOT HIER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Determine locations of turns;
+
+diffLocs = diff(locs); % calculates difference in step location
+avg_diffLocs = mean(diffLocs); % average distance between steps
+std_diffLocs = std(diffLocs); % standard deviation of distance between steps
+
+figure;
+findpeaks(diffLocs, 'MINPEAKHEIGHT', avg_diffLocs, 'MINPEAKDISTANCE',5); % these values have been chosen based on visual inspection of the signal
+line([1 length(diffLocs)],[avg_diffLocs avg_diffLocs])
+[pks_diffLocs, locs_diffLocs] = findpeaks(diffLocs, 'MINPEAKHEIGHT', avg_diffLocs,'MINPEAKDISTANCE',5);
+locsTurns =  [locs(locs_diffLocs),  locs(locs_diffLocs+1)];
+
+%% Visualizing turns
+
+% Duplying signal + visualing
+% to make second signal with the locations of the turns filled with NaN, so
+% that both signals can be plotted above each other in a different colour
+
+magNoGfilt_copy = magNoGfilt;
+for k = 1: size(locsTurns,1);
+    magNoGfilt_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
+end
+
+
+% visualising signal;
+figure;
+subplot(2,1,1)
+hold on;
+plot(magNoGfilt,'b')
+plot(magNoGfilt_copy, 'r');
+title('Inside Straight: Filtered data with turns highlighted in blue')
+hold off;
+
+%% Calculation
+% VRAAG LAURENS zie blauwe blaadje
+
+startPos = 1;
+for i = 1: size(locsTurns,1);
+    endPos = locsTurns(i,1)-1;
+    
+    inputData = dataAcc(startPos:endPos,:);
+    WindowLen = size(inputData,1);
+    ApplyRealignment = false;
+    [ResultStruct] = GaitOutcomesTrunkAccFuncIH(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps); % Naam van deze moet nog aangepast.
+    
+    if i ==1 % only the firs 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))])(i) = ResultStruct.([char(Parameters(j))]);
+    end
+    startPos = locsTurns(i,2)+1;
+    
+end
+
+clear ResultStruct;
+
+% Calculate mean over the bins without turns to get 1 outcome value per parameter for inside
+% straight;
+
+for j = 1:NrParameters;
+    ResultStruct.([char(Parameters(j))]) = nanmean(DataStraight.([char(Parameters(j))]))
+end
+
+
+end
+

GaitVariabilityAnalysisLD_v1.m

@@ -0,0 +1,103 @@
+%% 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
+% 
+% 
+
+

GaitVariabilityStructToExcelIH.m

@@ -0,0 +1,95 @@
+%% Description file
+
+%% Clear and close
+clear;
+close all;
+
+%% Load data;
+load('GaitVarOutcomes30april.mat')
+Participants = fields(OutcomesAcc);
+
+%% Settings;
+counter = 1;
+run = 1;
+ParticipantNr = 0;
+ConditionNames = {'TwoMWT','InsidePath','InsideStraight','Outside','Treadmill_Comfortable','Treadmill_Condition1','Treadmill_Condition2','Treadmill_Condition3','Treadmill_Condition4','Treadmill_Condition5','Treadmill_Condition6','Treadmill_Condition7','Treadmill_Condition8'};
+
+%% Reorder data;
+for i = 1: length(Participants);
+    % Participant Nr;
+    ParticipantNr = ParticipantNr + 1;
+    
+    % Group Nr;
+    if contains(Participants(i),'Y');
+        Group = 1;
+    else
+        Group = 2;
+    end
+    
+    WalkingConditions = fields(OutcomesAcc.([char(Participants(i))]));
+    
+    for j = 1: length(WalkingConditions);
+        
+        if strcmp(char(WalkingConditions(j)),'Treadmill')
+            
+            SubConditions = fieldnames(OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]));
+            
+            for k = 1: length(SubConditions);
+                
+                NameTreadmill = [char(WalkingConditions(j)),'_',char(SubConditions(k))];
+                ConditionNr = find(strcmp(ConditionNames, NameTreadmill));
+                Parameters = fieldnames(OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]).([char(SubConditions(k))]));
+                
+                for l = 1: length(Parameters);
+                    
+                    Data(counter+1,1) = {([char(Participants(i))])};
+                    Data(counter+1,2) = {ParticipantNr};
+                    Data(counter+1,3) = {Group};
+                    Data(counter+1,4) = {([char(NameTreadmill)])};
+                    Data(counter+1,5) = {ConditionNr};
+                    Data(counter+1,l+5) = {OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]).([char(SubConditions(k))]).([char(Parameters(l))])(1)};
+                    
+                end
+                counter = counter+1;
+                
+            end
+            
+        else
+            
+            ConditionNr = find(strcmp(ConditionNames, WalkingConditions(j)));
+            Parameters = fieldnames(OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]));
+            
+            for l = 1: length(Parameters);
+                
+                Data(counter+1,1) = {([char(Participants(i))])};
+                Data(counter+1,2) = {ParticipantNr};
+                Data(counter+1,3) = {Group};
+                Data(counter+1,4) = {([char(WalkingConditions(j))])};
+                Data(counter+1,5) = {ConditionNr};
+                Data(counter+1,l+5) = {OutcomesAcc.([char(Participants(i))]).([char(WalkingConditions(j))]).([char(Parameters(l))])(1)};
+                
+            end
+            
+            counter = counter+1;
+            
+        end
+        if run == 1;
+            Data(1,1) = {'ParticipantCode'};
+            Data(1,2) = {'ParticipantNr'};
+            Data(1,3) = {'Group'};
+            Data(1,4) = {'ConditionName'};
+            Data(1,5) = {'ConditionNr'};
+            Data(1,6:(size(Parameters,1)+5)) = Parameters'; % First row with variable names;
+            run = run+1;
+        end
+        
+    end
+end
+
+if ispc;
+    xlswrite('GaitVariabilityOutcomes.xls', 'Data');
+elseif ismac;
+    filename = 'GaitVariabilityOutcomes.xlsx';
+    writecell(Data,filename);
+else
+end

HarmonicityFrequency.m

@@ -0,0 +1,83 @@
+
+function [ResultStruct] = HarmonicityFrequency(dataAccCut_filt,P,F, StrideFrequency,dF,LowFrequentPowerThresholds,N_Harm,FS,AccVectorLen,ResultStruct)
+
+% Add sum of power spectra (as a rotation-invariant spectrum)
+P = [P,sum(P,2)];
+PS = sqrt(P);
+
+% Calculate the measures for the power per separate dimension
+for i=1:size(P,2);
+    % Relative cumulative power and frequencies that correspond to these cumulative powers
+    PCumRel = cumsum(P(:,i))/sum(P(:,i));
+    PSCumRel = cumsum(PS(:,i))/sum(PS(:,i));
+    FCumRel = F+0.5*dF;
+    
+    % Derive relative cumulative power for threshold frequencies
+    Nfreqs = size(LowFrequentPowerThresholds,2);
+    LowFrequentPercentage(i,1:Nfreqs) = interp1(FCumRel,PCumRel,LowFrequentPowerThresholds)*100;
+    
+    % Calculate relative power of first 10 harmonics, taking the power of each harmonic with a band of + and - 10% of the first
+    % harmonic around it
+    PHarm = zeros(N_Harm,1);
+    PSHarm = zeros(N_Harm,1);
+    for Harm = 1:N_Harm,
+        FHarmRange = (Harm+[-0.1 0.1])*StrideFrequency;
+        PHarm(Harm) = diff(interp1(FCumRel,PCumRel,FHarmRange));
+        PSHarm(Harm) = diff(interp1(FCumRel,PSCumRel,FHarmRange));
+    end
+    
+    % Derive index of harmonicity
+    if i == 2 % for ML we expect odd instead of even harmonics
+        IndexHarmonicity(i) = PHarm(1)/sum(PHarm(1:2:12));
+    elseif i == 4
+        IndexHarmonicity(i) = sum(PHarm(1:2))/sum(PHarm(1:12));
+    else
+        IndexHarmonicity(i) = PHarm(2)/sum(PHarm(2:2:12));
+    end
+    
+    % Calculate the phase speed fluctuations
+    StrideFreqFluctuation = nan(N_Harm,1);
+    StrSamples = round(FS/StrideFrequency);
+    for h=1:N_Harm,
+        CutOffs = [StrideFrequency*(h-(1/3)) , StrideFrequency*(h+(1/3))]/(FS/2);
+        if all(CutOffs<1) % for Stride frequencies above FS/20/2, the highest harmonics are not represented in the power spectrum
+            [b,a] = butter(2,CutOffs);
+            if i==4 % Take the vector length as a rotation-invariant signal
+                AccFilt = filtfilt(b,a,AccVectorLen);
+            else
+                AccFilt = filtfilt(b,a,dataAccCut_filt(:,i));
+            end
+            Phase = unwrap(angle(hilbert(AccFilt)));
+            SmoothPhase = sgolayfilt(Phase,1,2*(floor(FS/StrideFrequency/2))+1); % This is in fact identical to a boxcar filter with linear extrapolation at the edges
+            StrideFreqFluctuation(h) = std(Phase(1+StrSamples:end,1)-Phase(1:end-StrSamples,1));
+        end
+    end
+    FrequencyVariability(i) = nansum(StrideFreqFluctuation./(1:N_Harm)'.*PHarm)/nansum(PHarm);
+    
+    if i<4,
+        % Derive harmonic ratio (two variants)
+        if i == 2 % for ML we expect odd instead of even harmonics
+            HarmonicRatio(i) = sum(PSHarm(1:2:end-1))/sum(PSHarm(2:2:end)); % relative to summed 3d spectrum
+            HarmonicRatioP(i) = sum(PHarm(1:2:end-1))/sum(PHarm(2:2:end)); % relative to own spectrum
+        else
+            HarmonicRatio(i) = sum(PSHarm(2:2:end))/sum(PSHarm(1:2:end-1));
+            HarmonicRatioP(i) = sum(PHarm(2:2:end))/sum(PHarm(1:2:end-1));
+        end
+    end
+    
+end
+
+    ResultStruct.IndexHarmonicity_V = IndexHarmonicity(1); % higher smoother more regular patter
+    ResultStruct.IndexHarmonicity_ML = IndexHarmonicity(2); % higher smoother more regular patter
+    ResultStruct.IndexHarmonicity_AP = IndexHarmonicity(3); % higher smoother more regular patter
+    ResultStruct.IndexHarmonicity_All = IndexHarmonicity(4);
+    ResultStruct.HarmonicRatio_V = HarmonicRatio(1);
+    ResultStruct.HarmonicRatio_ML = HarmonicRatio(2);
+    ResultStruct.HarmonicRatio_AP = HarmonicRatio(3);
+    ResultStruct.HarmonicRatioP_V = HarmonicRatioP(1);
+    ResultStruct.HarmonicRatioP_ML = HarmonicRatioP(2);
+    ResultStruct.HarmonicRatioP_AP = HarmonicRatioP(3);
+    ResultStruct.FrequencyVariability_V = FrequencyVariability(1); 
+    ResultStruct.FrequencyVariability_ML = FrequencyVariability(2);
+    ResultStruct.FrequencyVariability_AP = FrequencyVariability(3); 
+    ResultStruct.StrideFrequency = StrideFrequency;

MutualInformationHisPro.m

@@ -0,0 +1,146 @@
+function [mutM,cummutM,minmuttauV] = MutualInformationHisPro(xV,tauV,bV,flag)
+% [mutM,cummutM,minmuttauV] =  MutualInformationHisPro(xV,tauV,bV,flag)
+% MUTUALINFORMATIONHISPRO computes the mutual information on the time
+% series 'xV' for given delays in 'tauV'. The estimation of mutual
+% information is based on 'b' partitions of equal probability at each dimension. 
+% A number of different 'b' can be given in the input vector 'bV'.
+% According to a given flag, it can also compute the cumulative mutual 
+% information for each given lag, as well as the time of the first minimum 
+% of the mutual information.
+% INPUT
+% - xV      : a vector for the time series
+% - tauV    : a vector of the delays to be evaluated for
+% - bV      : a vector of the number of partitions of the histogram-based
+%             estimate. 
+% - flag    : if 0-> compute only mutual information,
+%           : if 1-> compute the mutual information, the first minimum of
+%             mutual information and the cumulative mutual information. 
+%             if 2-> compute (also) the cumulative mutual information
+%             if 3-> compute (also) the first minimum of mutual information
+% OUTPUT
+% - mutM    : the vector of the mutual information values s for the given
+%             delays. 
+% - cummutM : the vector of the cumulative mutual information values for
+%             the given delays 
+% - minmuttauV : the time of the first minimum of the mutual information.
+%========================================================================
+%     <MutualInformationHisPro.m>, v 1.0 2010/02/11 22:09:14  Kugiumtzis & Tsimpiris
+%     This is part of the MATS-Toolkit http://eeganalysis.web.auth.gr/
+
+%========================================================================
+% Copyright (C) 2010 by Dimitris Kugiumtzis and Alkiviadis Tsimpiris 
+%                       <dkugiu@gen.auth.gr>
+
+%========================================================================
+% Version: 1.0
+
+% LICENSE:
+%     This program is free software; you can redistribute it and/or modify
+%     it under the terms of the GNU General Public License as published by
+%     the Free Software Foundation; either version 3 of the License, or
+%     any later version.
+%
+%     This program is distributed in the hope that it will be useful,
+%     but WITHOUT ANY WARRANTY; without even the implied warranty of
+%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%     GNU General Public License for more details.
+%
+%     You should have received a copy of the GNU General Public License
+%     along with this program. If not, see http://www.gnu.org/licenses/>.
+
+%=========================================================================
+% Reference : D. Kugiumtzis and A. Tsimpiris, "Measures of Analysis of Time Series (MATS): 
+% 	          A Matlab  Toolkit for Computation of Multiple Measures on Time Series Data Bases",
+%             Journal of Statistical Software, in press, 2010
+
+% Link      : http://eeganalysis.web.auth.gr/
+%========================================================================= 
+nsam = 1;
+n = length(xV);
+if nargin==3
+    flag = 1;
+elseif nargin==2 
+    flag = 1;
+    bV = round(sqrt(n/5));
+end
+if isempty(bV)
+    bV = round(sqrt(n/5));
+end
+bV(bV==0)=round(sqrt(n/5));
+tauV = sort(tauV);
+ntau = length(tauV);
+taumax = tauV(end);
+nb = length(bV);
+[oxV,ixV]=sort(xV);
+[tmpV,ioxV]=sort(ixV);
+switch flag
+    case 0
+        % Compute only the mutual information for the given lags
+        mutM = NaN*ones(ntau,nb);
+        for ib=1:nb
+            b = bV(ib);
+            if n<2*b
+                break;
+            end
+            mutM(:,ib)=mutinfHisPro(xV,tauV,b,ioxV,ixV);
+        end % for ib
+        cummutM=[];
+        minmuttauV=[];
+    case 1
+        % Compute the mutual information for all lags up to the
+        % largest given lag, then compute the lag of the first minimum of
+        % mutual information and the cumulative mutual information for the
+        % given lags.
+        mutM = NaN*ones(ntau,nb);
+        cummutM = NaN*ones(ntau,nb);
+        minmuttauV = NaN*ones(nb,1);
+        miM = NaN*ones(taumax+1,nb);
+        for ib=1:nb
+            b = bV(ib);
+            if n<2*b
+                break;
+            end
+            miM(:,ib)=mutinfHisPro(xV,[0:taumax]',b,ioxV,ixV);
+            mutM(:,ib) = miM(tauV+1,ib);
+            minmuttauV(ib) = findminMutInf(miM(:,ib),nsam);
+            % Compute the cumulative mutual information for the given delays
+            for i=1:ntau
+                cummutM(i,ib) = sum(miM(1:tauV(i)+1,ib));
+            end
+        end % for ib
+    case 2
+        % Compute the mutual information for all lags up to the largest 
+        % given lag and then sum up to get the cumulative mutual information 
+        % for the given lags.
+        cummutM = NaN*ones(ntau,nb);
+        miM = NaN*ones(taumax+1,nb);
+        for ib=1:nb
+            b = bV(ib);
+            if n<2*b
+                break;
+            end
+            miM(:,ib)=mutinfHisPro(xV,[0:taumax]',b,ioxV,ixV);
+            % Compute the cumulative mutual information for the given delays
+            for i=1:ntau
+                cummutM(i,ib) = sum(miM(1:tauV(i)+1,ib));
+            end
+        end % for ib
+        mutM = [];
+        minmuttauV=[];
+    case 3
+        % Compute the mutual information for all lags up to the largest 
+        % given lag and then compute the lag of the first minimum of the
+        % mutual information.
+        minmuttauV = NaN*ones(nb,1);
+        miM = NaN*ones(taumax+1,nb);
+        for ib=1:nb
+            b = bV(ib);
+            if n<2*b
+                break;
+            end
+            miM(:,ib)=mutinfHisPro(xV,[0:taumax]',b,ioxV,ixV);
+            minmuttauV(ib) = findminMutInf(miM(:,ib),nsam);
+        end % for ib
+        mutM = [];
+        cummutM=[];
+end

RealignSensorSignalHRAmp.m

@@ -0,0 +1,275 @@
+function [RealignedAcc,RotationMatrixT,Flags] = RealignSensorSignalHRAmp(Acc, FS)
+
+% History
+% 2013/05/31 SR solve problem of ML-AP rotation of appr. 90 degrees ("Maximum number of function evaluations has been exceeded") 
+% 2013/08/15 SR use literature HR (i.e. amplitude not power, and harmonic frequencies not sin^4 windows) 
+% 2013/08/16 SR correct: use abs(X) instead of sqrt(X.*X)
+% 2013/09/11 SR use only relevant F (containing odd or even harmonics) and speed up determination of harmonic frequencies 
+
+%% Define VT direction (RVT) as direction of mean acceleration
+MeanAcc = mean(Acc);
+RVT = MeanAcc'/norm(MeanAcc);
+RMLGuess = cross([0;0;1],RVT); RMLGuess = RMLGuess/norm(RMLGuess);
+RAPGuess = cross(RVT,RMLGuess); RAPGuess = RAPGuess/norm(RAPGuess);
+
+%% Estimate stride frequency
+StrideFrequency = StrideFrequencyFrom3dAcc(Acc, FS);
+AccDetrend = detrend(Acc,'constant');
+
+N = size(AccDetrend,1);
+
+%% calculate complex DFT
+F = FS*(0:(N-1))'/N;
+
+dF = FS/N;
+FHarmRange = [(1:20)'-0.1, (1:20)'+0.1]*StrideFrequency;
+
+EvenRanges = FHarmRange(2:2:end,:);
+OddRanges = FHarmRange(1:2:end,:);
+
+EvenHarmonics = zeros(size(F));
+for j=1:size(EvenRanges,1)
+    IXRange = [EvenRanges(j,1)/dF, EvenRanges(j,2)/dF]+1;
+    IXRangeRound = min(N,round(IXRange));
+    EvenHarmonics(IXRangeRound(1):IXRangeRound(2)) = EvenHarmonics(IXRangeRound(1):IXRangeRound(2)) + 1;
+    EvenHarmonics(IXRangeRound(1)) = EvenHarmonics(IXRangeRound(1)) - (IXRange(1)-IXRangeRound(1)+0.5);
+    EvenHarmonics(IXRangeRound(2)) = EvenHarmonics(IXRangeRound(2)) - (IXRangeRound(2)-IXRange(2)+0.5);
+end
+
+OddHarmonics = zeros(size(F));
+for j=1:size(OddRanges,1)
+    IXRange = [OddRanges(j,1)/dF, OddRanges(j,2)/dF]+1;
+    IXRangeRound = min(N,round(IXRange));
+    OddHarmonics(IXRangeRound(1):IXRangeRound(2)) = OddHarmonics(IXRangeRound(1):IXRangeRound(2)) + 1;
+    OddHarmonics(IXRangeRound(1)) = OddHarmonics(IXRangeRound(1)) - (IXRange(1)-IXRangeRound(1)+0.5);
+    OddHarmonics(IXRangeRound(2)) = OddHarmonics(IXRangeRound(2)) - (IXRangeRound(2)-IXRange(2)+0.5);
+end
+
+%% select only frequencies that contain odd or even harmonics
+RelevantF = OddHarmonics~=0 | EvenHarmonics~=0;
+OddHarmonics = OddHarmonics(RelevantF);
+EvenHarmonics = EvenHarmonics(RelevantF);
+DFT = fft(AccDetrend,N,1);
+DFT = DFT(RelevantF,:);
+DFTMLG = DFT*RMLGuess;
+DFTAPG = DFT*RAPGuess;
+
+%% Initial estimation of ML-AP realignment (educated guess)
+% We start of by finding the directions ML = MLguess + x*APguess for which 
+% the 'harmonic ratio' for the ML directions is maximal or minimal, and we 
+% do this by varying x in still the power variant of HR instead of amplitude:
+ HRML = @(x) real(((OddHarmonics'.*(DFTMLG'+x*DFTAPG'))*(DFTMLG+x*DFTAPG))...
+      /((EvenHarmonics'.*(DFTMLG'+x*DFTAPG'))*(DFTMLG+x*DFTAPG)));
+% or by substititutions
+a = real((OddHarmonics.*DFTAPG)'*DFTAPG);
+b = real((OddHarmonics.*DFTMLG)'*DFTAPG + (OddHarmonics.*DFTAPG)'*DFTMLG);
+c = real((OddHarmonics.*DFTMLG)'*DFTMLG);
+d = real((EvenHarmonics.*DFTAPG)'*DFTAPG);
+e = real((EvenHarmonics.*DFTMLG)'*DFTAPG + (EvenHarmonics.*DFTAPG)'*DFTMLG);
+f = real((EvenHarmonics.*DFTMLG)'*DFTMLG);
+% and by setting d/dx(HRML(x))==0 we get the quadratic formula
+aq = a*e-b*d;
+bq = 2*a*f-2*c*d;
+cq = b*f-c*e;
+% DT = bq^2-4*aq*cq = 4*(a*f-*c*d)*(a*f-c*d) - 4*(a*e-b*d)*(b*f-c*e)
+%                   = 4aaff +4ccdd -8acdf -4abef -4bcde +4acee + 4bbdf 
+%                   = 4aaff +4ccdd -4abef -4bcde +4ac(ee-df) +4(bb-ac)df 
+x = (-bq +[1 -1]*sqrt(bq.^2-4*aq*cq))/(2*aq);
+% with hopefully two solutions for x. 
+
+% from here on use amplitude HR:
+HRML = @(x) real((OddHarmonics'*abs(DFTMLG+x*DFTAPG))...
+     /(EvenHarmonics'*abs(DFTMLG+x*DFTAPG)));
+% now we take the x for which the HR is maximal, and the x rotated 90
+% degrees fo which it is minimal, and take the mean angle between them
+HR1 = HRML(x(1));
+HR2 = HRML(x(2));
+if HR1 > HR2
+    thetas = atan([x(1),-1/x(2)]);
+else
+    thetas = atan([-1/x(1),x(2)]);
+end
+theta = mean(thetas);
+xeg = tan(theta);
+if abs(diff(thetas)) > pi/2 % thetas are more than 90 degrees apart, and the 'mean of the two axes' should be rotated 90 degrees 
+    xeg = -1/xeg;
+end
+
+%% Numerical search for best harmonic ratios product
+% Now we want to improve the solution, since the two above directions are 
+% not necessarily orthogonal. Thus we want to find the orthogonal 
+% directions ML = MLguess + x*APguess, and AP = APguess - x*MLguess, 
+% for which the product of 'harmonic ratios' for the ML and AP directions 
+% is maximal, and we do this by varying x, starting with the educated guess
+% above using power instead of amplitude:
+%   HRML = ((OddHarmonics'.*(DFTMLG'+x*DFTAPG'))*(DFTMLG+x*DFTAPG))...
+%     /((EvenHarmonics'.*(DFTMLG'+x*DFTAPG'))*(DFTMLG+x*DFTAPG));
+%   HRAP = ((EvenHarmonics'.*(-x*DFTMLG'+DFTAPG'))*(-x*DFTMLG+DFTAPG))...
+%     /((OddHarmonics'.*(-x*DFTMLG'+DFTAPG'))*(-x*DFTMLG+DFTAPG));
+% and
+%   HR_Prod = HRML*HRAP = ...
+
+% Now use amplitude
+Minus_HR_Prod = @(x) -real(...
+                (OddHarmonics'*abs(DFTMLG+x*DFTAPG))...
+               /(EvenHarmonics'*abs(DFTMLG+x*DFTAPG))...
+               *(EvenHarmonics'*abs(-x*DFTMLG+DFTAPG))...
+               /(OddHarmonics'*abs(-x*DFTMLG+DFTAPG)));
+[xopt1,Minus_HR_Prod_opt1,exitflag1]=fminsearch(Minus_HR_Prod,xeg);
+if exitflag1 == 0 && abs(xopt1) > 100 % optimum for theta beyond +- pi/2
+    [xopt1a,Minus_HR_Prod_opt1a,exitflag1a]=fminsearch(Minus_HR_Prod,-xopt1);
+    if ~(exitflag1a == 0 && abs(xopt1a) > 100)
+        xopt1org = xopt1;
+        xopt1 = xopt1a;
+        Minus_HR_Prod_opt1 = Minus_HR_Prod_opt1a;
+    end
+end
+options = optimset('LargeScale','off','GradObj','off','Display','notify');
+[xopt2,Minus_HR_Prod_opt2,exitflag2]=fminunc(Minus_HR_Prod,xopt1,options);
+if exitflag2 == 0 && abs(xopt2) > 100 % optimum for theta beyond +- pi/2
+    [xopt2a,Minus_HR_Prod_opt2a,exitflag2a]=fminsearch(Minus_HR_Prod,-xopt2);
+    if ~(exitflag2a == 0 && abs(xopt2a) > 100)
+        xopt2 = xopt2a;
+        Minus_HR_Prod_opt2 = Minus_HR_Prod_opt2a;
+    end
+end
+
+RML = (RMLGuess+xopt2*RAPGuess)/sqrt(1+xopt2^2);
+RAP = cross(RVT,RML); RAP = RAP/norm(RAP);
+
+RT = [RVT,RML,RAP];
+if nargout > 1
+    RotationMatrixT = RT;
+end
+RealignedAcc = Acc*RT;
+
+if nargout > 2
+    Flags = [exitflag1,exitflag2];
+end
+
+if nargin > 2 && plotje
+    Thetas = pi*(-0.49:0.01:0.49)';
+    Xs = tan(Thetas);
+    HR_Prod = nan(size(Xs));
+    HRMLs = nan(size(Xs));
+    for i=1:numel(Xs),
+        HR_Prod(i,1) = -Minus_HR_Prod(Xs(i));
+        HRMLs(i,1) = HRML(Xs(i));
+    end
+    figure(); 
+    semilogy(Thetas,[HR_Prod,HRMLs.^2]);
+    hold on;
+    semilogy(atan(x),[HR1,HR2].^2,'r.');
+    semilogy(atan([xopt1;xopt2]),-[Minus_HR_Prod_opt1;Minus_HR_Prod_opt2],'g.');
+    semilogy(atan(xeg),-Minus_HR_Prod(xeg),'kx');
+end
+
+%% attempt to solve it analytically:
+% % We can rewrite HR_Prod as:
+% % HR_Prod = 
+% %    ( (x.^2)*((OddHarmonics.*DFTAPG)'*DFTAPG) ...
+% %     +  x*((OddHarmonics.*DFTMLG)'*DFTAPG + (OddHarmonics.*DFTAPG)'*DFTMLG)...
+% %     +  ((OddHarmonics.*DFTMLG)'*DFTMLG) )...
+% %  /
+% %     ( (x.^2)*((EvenHarmonics.*DFTAPG)'*DFTAPG) ...
+% %     +  x*((EvenHarmonics.*DFTMLG)'*DFTAPG + (EvenHarmonics.*DFTAPG)'*DFTMLG)...
+% %     +  ((EvenHarmonics.*DFTMLG)'*DFTMLG) )
+% %  *
+% %    ( (x.^2)*((EvenHarmonics.*DFTMLG)'*DFTMLG) ...
+% %     -  x*((EvenHarmonics.*DFTMLG)'*DFTAPG + (EvenHarmonics.*DFTAPG)'*DFTMLG)...
+% %     +  ((EvenHarmonics.*DFTAPG)'*DFTAPG) )...
+% %  /
+% %     ( (x.^2)*((OddHarmonics.*DFTMLG)'*DFTMLG) ...
+% %     -  x*((OddHarmonics.*DFTMLG)'*DFTAPG + (OddHarmonics.*DFTAPG)'*DFTMLG)...
+% %     +  ((OddHarmonics.*DFTAPG)'*DFTAPG) )
+% % 
+% % or by making these substitutions:
+% a2 = (OddHarmonics.*DFTAPG)'*DFTAPG;
+% a1 = (OddHarmonics.*DFTMLG)'*DFTAPG + (OddHarmonics.*DFTAPG)'*DFTMLG;
+% a0 = (OddHarmonics.*DFTMLG)'*DFTMLG;
+% c2 = (EvenHarmonics.*DFTAPG)'*DFTAPG;
+% c1 = (EvenHarmonics.*DFTMLG)'*DFTAPG + (EvenHarmonics.*DFTAPG)'*DFTMLG;
+% c0 = (EvenHarmonics.*DFTMLG)'*DFTMLG;
+% b2 = c0;
+% b1 = c1;
+% b0 = c2;
+% d2 = a0;
+% d1 = a1;
+% d0 = a2;
+% % we get:
+% %   HR_Prod = (x.^2*a2 + x*a1 + a0) * (x.^2*b2 + x*b1 + b0) / (x.^2*c2 + x*c1 + c0) * (x.^2*d2 + x*d1 + d0)
+% % or:
+% %   HR_Prod =  N/D 
+% % with
+% %   N =      ( x.^4*(a2*b2)...
+% %             +x.^3*(a2*b1+a1*b2)...
+% %             +x.^2*(a2*b0+a1*b1+a0*b2)...
+% %             +x   *(a1*b0+a0*b1)...
+% %             +      a0*b0              )...
+% % and with 
+% %   D =      ( x.^4*(c2*d2)...
+% %             +x.^3*(c2*d1+c1*d2)...
+% %             +x.^2*(c2*d0+c1*d1+c0*d2)...
+% %             +x   *(c1*d0+c0*d1)...
+% %             +      c0*d0              )
+% %   or D =   ( x.^4*(a0*b0)...
+% %             +x.^3*(a1*b0+a0*b1)...
+% %             +x.^2*(a2*b0+a1*b1+a0*b2)...
+% %             +x   *(a2*b1+a1*b2)...
+% %             +      a2*b2              )
+% % or by substition of
+% n4 = a2*b2;
+% n3 = a2*b1+a1*b2;
+% n2 = a2*b0+a1*b1+a0*b2;
+% n1 = a1*b0+a0*b1;
+% n0 = a0*b0;
+% % d4 = n0;
+% % d3 = n1;
+% % d2 = n2;
+% % d1 = n3;
+% % d0 = n4;
+% % we get 
+% % N = sum(ni*x^i) and D=sum(di*x^i)
+% % 
+% % Then HR_Prod reaches a maximum or minimum when d/dx (N/D) = 0 
+% % or (dN/dx)/D - N/(D^2)*(dD/dx) = 0
+% % or (dN/dx)*D - (dD/dx)*N = 0
+% % or
+% %    x^7*(4*n4*d4 - 4*n4*d4)
+% %  + x^6*(4*n4*d3 - 4*n3*d4 + 3*n3*d4 - 3*n4*d3)
+% %  + x^5*(4*n4*d2 - 4*n2*d4 + 3*n3*d3 - 3*n3*d3 + 2*n2*d4 - 2*n4*d2)
+% %  + x^4*(4*n4*d1 - 4*n1*d4 + 3*n3*d2 - 3*n2*d3 + 2*n2*d3 - 2*n3*d2 + 1*n1*d4 - 1*n4*d1)
+% %  + x^3*(4*n4*d0 - 4*n0*d4 + 3*n3*d1 - 3*n1*d3 + 2*n2*d2 - 2*n2*d2 + 1*n1*d3 - 1*n3*d1)
+% %  + x^2*(                    3*n3*d0 - 3*n0*d3 + 2*n2*d1 - 2*n1*d2 + 1*n1*d2 - 1*n2*d1)
+% %  + x^1*(                                        2*n2*d0 - 2*n0*d2 + 1*n1*d1 - 1*n1*d1)
+% %  + x^0*(                                                            1*n1*d0 - 1*n0*d1)
+% %  = 0 
+% % or
+% %    x^6*(1*n4*d3 - 1*n3*d4)
+% %  + x^5*(2*n4*d2 - 2*n2*d4)