170 lines
7.7 KiB
Matlab
170 lines
7.7 KiB
Matlab
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 LD feb 2021 (IH feb 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) 5360 [ 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)
|
||
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; % ALREADY REORDERD: 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 that data for CONTROLS; reorder data to 1 = V; 2 = ML; 3 = AP
|
||
%data = inputData(:,[3,2,1]);
|
||
%[RealignedAcc, ~] = RealignSensorSignalHRAmp(data, FS); might not be necessary
|
||
%dataAcc = RealignedAcc;
|
||
data = inputData;
|
||
dataAcc = inputData;
|
||
[B,A] = butter(2,20/(FS/2),'low');
|
||
dataAcc_filt = filtfilt(B,A,inputData);
|
||
end
|
||
|
||
|
||
%% Step detection
|
||
% Determines the number of steps in the signal so that the first 1 and last step 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(data,StrideTimeSamples,FS);
|
||
% This function selects steps based on negative and positive values.
|
||
% However to determine the steps correctly we only need one of these;
|
||
LocsStepsLD = PksAndLocsCorrected;
|
||
|
||
%% Cut data & remove currents results
|
||
% Remove 1 step in the beginning and end of data
|
||
dataAccCut = dataAcc(LocsStepsLD(1):LocsStepsLD(end-1),:);
|
||
dataAccCut_filt = dataAcc_filt(LocsStepsLD(1):LocsStepsLD(end-1),:);
|
||
|
||
% Clear currently saved results from Autocorrelation Analysis
|
||
|
||
clear ResultStruct;
|
||
clear PksAndLocsCorrected;
|
||
clear LocsSteps;
|
||
|
||
% Change window length necessary if ApplyRemoveSteps? (16-2-2013 LD)
|
||
|
||
WindowLen = 10*FS;
|
||
|
||
else;
|
||
dataAccCut = dataAcc;
|
||
dataAccCut_filt = dataAcc_filt;
|
||
end
|
||
|
||
%% Calculate stride parameters
|
||
ResultStruct = struct; % create empty struct
|
||
|
||
% Run function AutoCorrStrides, Outcomeparameters: StrideRegularity AP/VT ,StrideTimeSamples,StrideTime
|
||
[ResultStruct] = AutocorrStrides(dataAccCut_filt,FS, StrideTimeRange,ResultStruct);
|
||
StrideTimeSamples = ResultStruct.StrideTimeSamples; % needed as input for other functions
|
||
|
||
%% 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)); % WindowLen -> 10*Fs OR on InputSignal
|
||
[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); % Not sure if WindowLen should be different?
|
||
N_SkipBegin = ceil((size(dataAccCut,1)-N_Windows*WindowLen)/2);
|
||
LyapunovWolf = 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));
|
||
[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);
|
||
SampleEntropy = nanmean(SE,1);
|
||
|
||
ResultStruct.LyapunovWolf_V = LyapunovWolf(1);
|
||
ResultStruct.LyapunovWolf_ML = LyapunovWolf(2);
|
||
ResultStruct.LyapunovWolf_AP = LyapunovWolf(3);
|
||
ResultStruct.SampleEntropy_V = SampleEntropy(1);
|
||
ResultStruct.SampleEntropy_ML = SampleEntropy(2);
|
||
ResultStruct.SampleEntropy_AP = SampleEntropy(3);
|
||
|
||
%% Calculate RMS in each direction: added february 2021 by LD, CONSTRUCT: 'Pace'
|
||
% Sekine, M., Tamura, T., Yoshida, M., Suda, Y., Kimura, Y., Miyoshi, H., ... & Fujimoto, T. (2013). A gait abnormality measure based on root mean square of trunk acceleration. Journal of neuroengineering and rehabilitation, 10(1), 1-7.
|
||
|
||
Data_Centered = normalize(dataAcc_filt,'center','mean'); % The RMS coincides with the Sd since the Acc signals are transformed to give a mean equal to zero
|
||
RMS = rms(Data_Centered);
|
||
|
||
ResultStruct.RMS_V = RMS(1);
|
||
ResultStruct.RMS_ML = RMS(2);
|
||
ResultStruct.RMS_AP = RMS(3);
|
||
|
||
end |