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'); % The unbiased alternative is produced by dividing the raw autocorrelation coefficient by the number of samples representing the overlapping part of the time series and the time-lagged replication:
AutocorrSum = sum(AutocorrResult(:,[1 5 9]),2); % This sum is independent of sensor re-orientation, as long as axes are kept orthogonal. 1: VT-VT 5:ML-ML 9:AP-AP
AutocorrResult2= [AutocorrResult(:,[1 5 9]),AutocorrSum]; % column 4 represents the sum of autocorr [1 ... 9]
IXRange = (numel(Lags)-(RangeEnd-RangeStart)):numel(Lags); % Lags from -400 / 400, IXrange: 421 - 801
% 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 / stride regularity was shifted to the average stride time
    RelativeStrideVariability = 1-StrideRegularity;
    StrideTimeSeconds = StrideTimeSamples/FS;
    
    ResultStruct.StrideRegularity_V = StrideRegularity(1);
    ResultStruct.StrideRegularity_AP = StrideRegularity(3);
    ResultStruct.StrideTimeSamples = StrideTimeSamples;
    %ResultStruct.StrideTimeSeconds = StrideTimeSeconds;
    
end