LBP&Controls #1
@ -4,10 +4,10 @@ function [ResultStruct] = AutocorrStrides(data,FS, StrideTimeRange,ResultStruct
|
||||
% 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);
|
||||
[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]);
|
||||
@ -25,19 +25,13 @@ if isempty(IXRangeNew)
|
||||
|
||||
else
|
||||
StrideTimeSamples = Lags(IXRangeNew(AutocorrSum(IXRangeNew)==max(AutocorrSum(IXRangeNew))));
|
||||
StrideRegularity = AutocorrResult2(Lags==StrideTimeSamples,:)./AutocorrResult2(Lags==0,:); % Moe-Nilssen&Helbostatt,2004
|
||||
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_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;
|
||||
%ResultStruct.StrideTimeSeconds = StrideTimeSeconds;
|
||||
|
||||
end
|
||||
@ -42,10 +42,10 @@ function [L_Estimate,ExtraArgsOut] = CalcMaxLyapWolfFixedEvolv(ThisTimeSeries,FS
|
||||
% 23 October 2012, use fixed evolve time instead of adaptable
|
||||
|
||||
if nargin > 2
|
||||
if isfield(ExtraArgsIn,'J')
|
||||
if isfield(ExtraArgsIn,'J') % logical 0, based on MutualInformation.
|
||||
J=ExtraArgsIn.J;
|
||||
end
|
||||
if isfield(ExtraArgsIn,'m')
|
||||
if isfield(ExtraArgsIn,'m') % logical 1, embedding dimension of 7.
|
||||
m=ExtraArgsIn.m;
|
||||
end
|
||||
end
|
||||
@ -75,9 +75,9 @@ if ~(nanstd(ThisTimeSeries) > 0)
|
||||
return;
|
||||
end
|
||||
|
||||
%% Determine J
|
||||
%% Determine J (Embedding Delay as the first minimum of the average mutual information : Fraser, A.M., Swinney, H.L., 1986. Independent coordinates for strange attractors from mutual information. Phys. Rev. A 33, 1134–1140.)
|
||||
if ~exist('J','var') || isempty(J)
|
||||
% Calculate mutual information and take first local minimum Tau as J
|
||||
% Calculate mutual information and take first local minimum Tau as J; To obtain coordinates for time delayed phase-space embedding that are as independent as possible
|
||||
bV = min(40,floor(sqrt(size(ThisTimeSeries,1))));
|
||||
tauVmax = 70;
|
||||
[mutMPro,cummutMPro,minmuttauVPro] = MutualInformationHisPro(ThisTimeSeries,(0:tauVmax),bV,1); % (xV,tauV,bV,flag)
|
||||
@ -113,6 +113,9 @@ end
|
||||
ExtraArgsOut.m=m;
|
||||
|
||||
%% Create state space based upon J and m
|
||||
% the procedure of computing the LYE involves first the calculation of a phasespace reconstruction for each acceleration time-series,
|
||||
% by creating time-delayed copies of the time-series X(t) , X(t1 + tau),
|
||||
% x(t2 + 2tau), ...
|
||||
N_ss = size(ThisTimeSeries,1)-(m-1)*J;
|
||||
StateSpace=nan(N_ss,m*size(ThisTimeSeries,2));
|
||||
for dim=1:size(ThisTimeSeries,2),
|
||||
@ -121,13 +124,13 @@ for dim=1:size(ThisTimeSeries,2),
|
||||
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;
|
||||
%% Parameters for Lyapunov estimation % SEE: Rispens, S. M., Pijnappels, M. A. G. M., Van Dieën, J. H., Van Schooten, K. S., Beek, P. J., & Daffertshofer, A. (2014). A benchmark test of accuracy and precision in estimating dynamical systems characteristics from a time series. Journal of biomechanics, 47(2), 470-475.
|
||||
CriticalLen=J*m; % After a fixed period (here the embedding delay) the nearby point is replaced.
|
||||
max_dist = sqrt(sum(std(StateSpace).^2))/10; % The maximum distance to which the divergence may grow before replacing the neighbour.
|
||||
max_dist_mult = 5; % Whenever points cannot be found, the distance limit is increased stepwise to maximally five times the original limit.
|
||||
min_dist = max_dist/2; % The new nearby point must not be closer than 1/20 of the signal's SD to avoid noise-related artifacts.
|
||||
max_theta = 0.3; % The new nearby point has to point in the same direction in state space as the old nearby point, as seen from the current reference (at least within 0.3 rad)
|
||||
evolv = J; % Embedding delay (time step after which the neighbour is replaced) (approximately 10 % of gait cycle?)
|
||||
|
||||
%% Calculate Lambda
|
||||
[L_Estimate]=div_wolf_fixed_evolv(StateSpace, FS, min_dist, max_dist, max_dist_mult, max_theta, CriticalLen, evolv);
|
||||
|
||||
@ -1,64 +1,85 @@
|
||||
function [locsTurns,FilteredData] = DetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,Distance)
|
||||
function [locsTurns,FilteredData,FootContacts] = DetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,Distance,ResultStruct)
|
||||
% Description: Determine the location of turns, plot for visual inspection
|
||||
|
||||
% Input: Acc Data (not yet realigned)
|
||||
% Input: Acc Data (in previous stage realigned?)
|
||||
% Realigned sensor data to VT-ML-AP frame
|
||||
|
||||
%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).
|
||||
if ApplyRealignment % apply realignment 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;
|
||||
else % realignment/detrend performed in earlier stage!
|
||||
data = inputData;
|
||||
dataAcc = inputData;
|
||||
end
|
||||
% Filter data
|
||||
[B,A] = butter(2,3/(FS/2),'low'); % Filters data very strongly which is needed to determine turns correctly
|
||||
[B,A] = butter(2,20/(FS/2),'low'); % Filters data very strongly which is needed to determine turns correctly
|
||||
dataStepDetection = filtfilt(B,A,dataAcc);
|
||||
% Determine steps
|
||||
% 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,
|
||||
VTAcc = data(:,1);
|
||||
APAcc = data(:,3);
|
||||
|
||||
%% Determine Steps -
|
||||
% USE THE AP ACCELERATION: WIEBREN ZIJLSTRA: ASSESSMENT OF SPATIO-TEMPORAL PARAMTERS DURING UNCONSTRAINED WALKING (2004)
|
||||
|
||||
% In order to run the step detection script we first need to run an autocorrelation function;
|
||||
StrideTimeRange = [0.2 4.0]; % Range to search for stride time (seconds)
|
||||
|
||||
[ResultStruct] = AutocorrStrides(dataStepDetection,FS, StrideTimeRange,ResultStruct); % first guess of stridetimesamples
|
||||
|
||||
% StrideTimeSamples is needed as an input for the stepcountFunc;
|
||||
StrideTimeSamples = ResultStruct.StrideTimeSamples;
|
||||
|
||||
[PksAndLocsCorrected] = StepcountFunc(data,StrideTimeSamples,FS);
|
||||
FootContacts = PksAndLocsCorrected; % previous version LD: (1:2:end,2);
|
||||
numStepsOption2_filt = numel(FootContacts); % counts number of steps;
|
||||
|
||||
%% Determine Turns
|
||||
% 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 = 1.2*std(magNoGfilt); % based on visual inspection, parameter tuning was performed on standard deviation from MInPeak (used to be 1 SD)
|
||||
[pks, locs] = findpeaks(magNoGfilt, 'MINPEAKHEIGHT', minPeakHeight2); % for step detection
|
||||
numStepsOption2_filt = numel(pks); % counts number of steps;
|
||||
minPeakHeight2 = 1*std(magNoGfilt); % used to be 2 % based on visual inspection, parameter tuning was performed on *X*standard deviation from MInPeak (used to be 1 SD)
|
||||
[pks, locs] = findpeaks(magNoGfilt, 'MINPEAKHEIGHT', minPeakHeight2,'MINPEAKDISTANCE',0.4*FS); % for step detection
|
||||
numStepsOption2_filt = numel(pks); % counts number of locs for turn detection
|
||||
|
||||
diffLocs = diff(locs); % calculates difference in step location
|
||||
diffLocs = diff(locs); % calculates difference in location magnitude differences
|
||||
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',9); % these values have been chosen based on visual inspection of the signal
|
||||
findpeaks(diffLocs, 'MINPEAKHEIGHT', avg_diffLocs, 'MINPEAKDISTANCE',0.2*FS); % 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',10); % values were initially 5
|
||||
[pks_diffLocs, locs_diffLocs] = findpeaks(diffLocs, 'MINPEAKHEIGHT', avg_diffLocs,'MINPEAKDISTANCE',0.2*FS); % values were initially 5
|
||||
locsTurns = [locs(locs_diffLocs), locs(locs_diffLocs+1)];
|
||||
|
||||
magNoGfilt_copy = magNoGfilt;
|
||||
%magNoGfilt_copy = magNoGfilt;
|
||||
VTAcc_copy = data(:,1);
|
||||
APAcc_copy = data(:,3);
|
||||
for k = 1: size(locsTurns,1);
|
||||
magNoGfilt_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
|
||||
VTAcc_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
|
||||
APAcc_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
|
||||
end
|
||||
|
||||
% Visualising signal;
|
||||
if plotit
|
||||
figure;
|
||||
figure()
|
||||
subplot(2,1,1)
|
||||
hold on;
|
||||
plot(magNoGfilt,'b')
|
||||
plot(magNoGfilt_copy, 'r');
|
||||
title('Inside Straight: Filtered data with turns highlighted in blue')
|
||||
line([6000,12000,18000,24000,30000,36000;6000,12000,18000,24000,30000,36000],[-4,-4,-4,-4,-4,-4;4,4,4,4,4,4],'LineWidth',2,'Linestyle','--','color','k')
|
||||
% hold on;
|
||||
% for m = 1:size(locsTurns,1)
|
||||
% plot(locsTurns(m),DataStraight.([char(Participants(i))]).LyapunovPerStrideRosen_ML(:,m),'--gs',...
|
||||
% 'LineWidth',2,...
|
||||
% 'MarkerSize',10,...
|
||||
% 'MarkerEdgeColor','g',...
|
||||
% 'MarkerFaceColor',[0.5,0.5,0.5])
|
||||
% end
|
||||
plot(VTAcc,'b')
|
||||
plot(VTAcc_copy, 'r');
|
||||
title('Inside Straight: Filtered data VT with turns highlighted in blue')
|
||||
line([6000,12000,18000,24000,30000,36000;6000,12000,18000,24000,30000,36000],[-7,-7,-7,-7,-7,-7;7,7,7,7,7,7],'LineWidth',2,'Linestyle','--','color','k')
|
||||
hold on;
|
||||
subplot(2,1,2)
|
||||
plot(APAcc,'g')
|
||||
plot(APAcc_copy, 'r');
|
||||
title('Inside Straight: Filtered data AP with turns highlighted in blue')
|
||||
line([6000,12000,18000,24000,30000,36000;6000,12000,18000,24000,30000,36000],[-7,-7,-7,-7,-7,-7;7,7,7,7,7,7],'LineWidth',2,'Linestyle','--','color','k')
|
||||
hold off;
|
||||
end
|
||||
|
||||
% Check if number of turns * 20 m are making sense based on total
|
||||
% distance measured by researcher.
|
||||
disp(['Number of turns detected = ' num2str(size(locsTurns,1))])
|
||||
|
||||
@ -47,7 +47,6 @@ N_Harm = 12; % Number of harmonics used for harmonic ratio
|
||||
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];
|
||||
@ -59,20 +58,24 @@ ResultStruct = struct();
|
||||
% 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%
|
||||
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 is already detrended and in order 1 = AP, 2 = ML, 3 = VT in an earlier stage;
|
||||
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 = inputData(:, [3,2,1]);
|
||||
dataAcc_filt = filtfilt(B,A,dataAcc);
|
||||
dataAcc_filt = filtfilt(B,A,inputData);
|
||||
end
|
||||
|
||||
|
||||
%% Step dectection
|
||||
%% 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
|
||||
@ -84,15 +87,15 @@ if ApplyRemoveSteps
|
||||
StrideTimeSamples = ResultStruct.StrideTimeSamples;
|
||||
|
||||
% Calculate the number of steps;
|
||||
[PksAndLocsCorrected] = StepcountFunc(dataAcc_filt,StrideTimeSamples,FS);
|
||||
[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;
|
||||
LocsSteps = PksAndLocsCorrected(1:2:end,2);
|
||||
LocsStepsLD = PksAndLocsCorrected;
|
||||
|
||||
%% Cut data & remove currents results
|
||||
% Remove 1 step in the beginning and end of data
|
||||
dataAccCut = dataAcc(LocsSteps(1):LocsSteps(end-1),:);
|
||||
dataAccCut_filt = dataAcc_filt(LocsSteps(1):LocsSteps(end-1),:);
|
||||
dataAccCut = dataAcc(LocsStepsLD(1):LocsStepsLD(end-1),:);
|
||||
dataAccCut_filt = dataAcc_filt(LocsStepsLD(1):LocsStepsLD(end-1),:);
|
||||
|
||||
% Clear currently saved results from Autocorrelation Analysis
|
||||
|
||||
@ -100,9 +103,9 @@ if ApplyRemoveSteps
|
||||
clear PksAndLocsCorrected;
|
||||
clear LocsSteps;
|
||||
|
||||
% Change window length if ApplyRemoveSteps (16-2-2021 LD)
|
||||
% Change window length necessary if ApplyRemoveSteps? (16-2-2013 LD)
|
||||
|
||||
WindowLen = size(dataAccCut,1);
|
||||
WindowLen = 10*FS;
|
||||
|
||||
else;
|
||||
dataAccCut = dataAcc;
|
||||
@ -112,49 +115,10 @@ end
|
||||
%% Calculate stride parameters
|
||||
ResultStruct = struct; % create empty struct
|
||||
|
||||
% Run function AutoCorrStrides, Outcomeparameters: StrideRegularity,RelativeStrideVariability,StrideTimeSamples,StrideTime
|
||||
% 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 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
|
||||
@ -162,53 +126,45 @@ ResultStruct.GaitSymmIndex = nanmean(SI);
|
||||
|
||||
%% Measures derived from spectral analysis
|
||||
|
||||
AccVectorLen = sqrt(sum(dataAccCut_filt(:,1:3).^2,2));
|
||||
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);
|
||||
N_SkipBegin = ceil((size(dataAccCut,1)-N_Windows*WindowLen)/2);
|
||||
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);
|
||||
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
|
||||
%% 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.
|
||||
|
||||
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);
|
||||
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
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because it is too large
Load Diff
@ -1,103 +0,0 @@
|
||||
%% 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
|
||||
%
|
||||
%
|
||||
|
||||
|
||||
@ -1,6 +1,9 @@
|
||||
|
||||
function [ResultStruct] = HarmonicityFrequency(dataAccCut_filt,P,F, StrideFrequency,dF,LowFrequentPowerThresholds,N_Harm,FS,AccVectorLen,ResultStruct)
|
||||
|
||||
%LD 15-04-2021: Solve problem finding fundamental frequency in power
|
||||
%spectral density STRIDE - VS STEP frequency.
|
||||
|
||||
% Add sum of power spectra (as a rotation-invariant spectrum)
|
||||
P = [P,sum(P,2)];
|
||||
PS = sqrt(P);
|
||||
@ -8,6 +11,7 @@ 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;
|
||||
@ -21,12 +25,12 @@ for i=1:size(P,2);
|
||||
PHarm = zeros(N_Harm,1);
|
||||
PSHarm = zeros(N_Harm,1);
|
||||
for Harm = 1:N_Harm,
|
||||
FHarmRange = (Harm+[-0.1 0.1])*StrideFrequency;
|
||||
FHarmRange = (Harm+[-0.1 0.1])*(StrideFrequency); % In view of possible drift, which could lead to missing or widening peaks, the power spectral density of each peak was calculated within the frequency bands of +0.1 and −0.1 Hz of the peak frequency value.
|
||||
PHarm(Harm) = diff(interp1(FCumRel,PCumRel,FHarmRange));
|
||||
PSHarm(Harm) = diff(interp1(FCumRel,PSCumRel,FHarmRange));
|
||||
end
|
||||
|
||||
% Derive index of harmonicity
|
||||
% Derive index of harmonicity; which is the spectral power of the basic harmonic divided by the sum of the power of first six harmonics
|
||||
if i == 2 % for ML we expect odd instead of even harmonics
|
||||
IndexHarmonicity(i) = PHarm(1)/sum(PHarm(1:2:12));
|
||||
elseif i == 4
|
||||
@ -55,7 +59,8 @@ for i=1:size(P,2);
|
||||
FrequencyVariability(i) = nansum(StrideFreqFluctuation./(1:N_Harm)'.*PHarm)/nansum(PHarm);
|
||||
|
||||
if i<4,
|
||||
% Derive harmonic ratio (two variants)
|
||||
% Derive harmonic ratio (two variants) --> Ratio of the even and
|
||||
% odd harmonics!
|
||||
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
|
||||
@ -70,14 +75,6 @@ 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;
|
||||
Binary file not shown.
@ -22,7 +22,7 @@ function [mutM,cummutM,minmuttauV] = MutualInformationHisPro(xV,tauV,bV,flag)
|
||||
% 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.
|
||||
% - : 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/
|
||||
|
||||
BIN
Phyphoxdata.mat
BIN
Phyphoxdata.mat
Binary file not shown.
@ -1,17 +1,23 @@
|
||||
function [ResultStruct] = SpatioTemporalGaitParameters(dataAccCut_filt,StrideTimeSamples,ApplyRealignment,LegLength,FS,IgnoreMinMaxStrides,ResultStruct);
|
||||
%% Method Zijlstra & Hof
|
||||
% Mean step length and mean walking speed are estimated using the upward
|
||||
% and downward momvements of the trunk. Assuming a compass gait type, CoM
|
||||
% movements in the sagittal plane follow a circular trajectory during each
|
||||
% single support phase. In this inverted pendlum model, changes in height
|
||||
% of CoM depend on step length. Thus, when changes in height are known,
|
||||
% step length can be predicted from geometrical characteristics as in
|
||||
|
||||
Cutoff = 0.1;
|
||||
MinDist = floor(0.7*0.5*StrideTimeSamples); % Use StrideTimeSamples estimated above
|
||||
DatalFilt = dataAccCut_filt;
|
||||
|
||||
DatalFilt = dataAccCut_filt;
|
||||
% From acceleration to velocity
|
||||
Vel = cumsum(detrend(DatalFilt,'constant'))/FS;
|
||||
[B,A] = butter(2,Cutoff/(FS/2),'high');
|
||||
[B,A] = butter(2,Cutoff/(FS/2),'high'); % To avoid integration drift, position data were high-pass filtered
|
||||
Pos = cumsum(filtfilt(B,A,Vel))/FS;
|
||||
PosFilt = filtfilt(B,A,Pos);
|
||||
PosFiltVT = PosFilt(:,1);
|
||||
PosFiltVT = PosFilt(:,1); % Changes in vertical position were calculated by a double integration of Y.
|
||||
|
||||
% Find minima and maxima in vertical position
|
||||
|
||||
% if ~ApplyRealignment % Signals were not realigned, so it has to be done here
|
||||
% MeanAcc = mean(AccLoco);
|
||||
@ -28,7 +34,8 @@ if isempty(PosPks) && isempty(NegPks)
|
||||
else
|
||||
PksAndLocs = sortrows([PosPks,PosLocs,ones(size(PosPks)) ; NegPks,NegLocs,-ones(size(NegPks))], 2);
|
||||
end
|
||||
% Correct events for two consecutive maxima or two consecutive minima
|
||||
%% Correct events for two consecutive maxima or two consecutive minima
|
||||
|
||||
Events = PksAndLocs(:,2);
|
||||
NewEvents = PksAndLocs(:,2);
|
||||
Signs = PksAndLocs(:,3);
|
||||
@ -77,9 +84,9 @@ DH(DH>MaxDH) = MaxDH;
|
||||
% (Use delta h and delta t to calculate walking speed: use formula from
|
||||
% Z&H, but divide by 2 (skip factor 2)since we get the difference twice
|
||||
% each step, and multiply by 1.25 which is the factor suggested by Z&H)
|
||||
HalfStepLen = 1.25*sqrt(2*LegLength*DH-DH.^2);
|
||||
HalfStepLen = 1.25*sqrt(2*LegLength*DH-DH.^2); % Formulae is correct!
|
||||
Distance = sum(HalfStepLen);
|
||||
WalkingSpeedMean = Distance/(sum(DT)/FS);
|
||||
WalkingSpeedMean = Distance/(sum(DT)/FS); % Gait Speed (m/s) average walking speed
|
||||
% Estimate variabilities between strides
|
||||
StrideLengths = HalfStepLen(1:end-3) + HalfStepLen(2:end-2) + HalfStepLen(3:end-1) + HalfStepLen(4:end);
|
||||
StrideTimes = PksAndLocsCorrected(5:end,2)-PksAndLocsCorrected(1:end-4,2);
|
||||
@ -91,28 +98,25 @@ for i=1:4,
|
||||
WSS(i) = std(StrideSpeeds(i:4:end));
|
||||
end
|
||||
|
||||
StepLengthMean=mean(StrideLengths);
|
||||
StepLengthMean=mean(StrideLengths)/2; % 9-4-2021 LD, previous version with dividing by 2, THEN THIS SHOULD BE STRIDELENGTHMEAN!!
|
||||
StrideLengthMean = mean(StrideLengths);
|
||||
StrideTimeMean = mean(StrideTimes)/FS; % Samples to seconds
|
||||
StrideSpeedMean = mean(StrideSpeeds);
|
||||
|
||||
StrideTimeVariability = min(STS);
|
||||
StrideSpeedVariability = min(WSS);
|
||||
StrideLengthVariability = std(StrideLengths);
|
||||
% Estimate Stride time variability and stride speed variability by removing highest and lowest part
|
||||
if ~isinteger(IgnoreMinMaxStrides)
|
||||
IgnoreMinMaxStrides = ceil(IgnoreMinMaxStrides*size(StrideTimes,1));
|
||||
end
|
||||
StrideTimesSorted = sort(StrideTimes);
|
||||
StrideTimeVariabilityOmitOutlier = std(StrideTimesSorted(1+IgnoreMinMaxStrides:end-IgnoreMinMaxStrides));
|
||||
StrideSpeedSorted = sort(StrideSpeeds);
|
||||
StrideSpeedVariabilityOmitOutlier = std(StrideSpeedSorted(1+IgnoreMinMaxStrides:end-IgnoreMinMaxStrides));
|
||||
StrideLengthsSorted = sort(StrideLengths);
|
||||
StrideLengthVariabilityOmitOutlier = std(StrideLengthsSorted(1+IgnoreMinMaxStrides:end-IgnoreMinMaxStrides));
|
||||
|
||||
ResultStruct.StepLengthMean = StepLengthMean;
|
||||
ResultStruct.Distance = Distance;
|
||||
ResultStruct.WalkingSpeedMean = WalkingSpeedMean;
|
||||
ResultStruct.StrideTimeVariability = StrideTimeVariability;
|
||||
ResultStruct.StrideSpeedVariability = StrideSpeedVariability;
|
||||
ResultStruct.StrideLengthVariability = StrideLengthVariability;
|
||||
ResultStruct.StrideTimeVariabilityOmitOutlier = StrideTimeVariabilityOmitOutlier;
|
||||
ResultStruct.StrideSpeedVariabilityOmitOutlier = StrideSpeedVariabilityOmitOutlier;
|
||||
ResultStruct.StrideLengthVariabilityOmitOutlier = StrideLengthVariabilityOmitOutlier;
|
||||
% LD 16-03-2021 - Calculate CoV as ''coefficient of relative variability''
|
||||
CoVStrideTime = (nanstd(StrideTimes)/nanmean(StrideTimes))*100;
|
||||
CoVStrideSpeed = (nanstd(StrideSpeeds)/nanmean(StrideSpeeds))*100;
|
||||
CoVStrideLength = (nanstd(StrideLengths)/nanmean(StrideLengths))*100;
|
||||
|
||||
% ResultStruct
|
||||
ResultStruct.WalkingSpeedMean = WalkingSpeedMean;
|
||||
|
||||
ResultStruct.StrideTimeMean = StrideTimeMean;
|
||||
ResultStruct.CoVStrideTime = CoVStrideTime;
|
||||
|
||||
ResultStruct.StrideLengthMean = StrideLengthMean;
|
||||
ResultStruct.CoVStrideLength = CoVStrideLength;
|
||||
@ -2,14 +2,18 @@ function [ResultStruct] = SpectralAnalysisGaitfunc(dataAccCut_filt,WindowLen,FS,
|
||||
|
||||
P=zeros(0,size(dataAccCut_filt,2));
|
||||
|
||||
for i=1:size(dataAccCut_filt,2)
|
||||
for i=1:size(dataAccCut_filt,2) % pwelch = Welch’s power spectral density estimate
|
||||
[P1,~] = pwelch(dataAccCut_filt(:,i),hamming(WindowLen),[],WindowLen,FS);
|
||||
[P2,F] = pwelch(dataAccCut_filt(end:-1:1,i),hamming(WindowLen),[],WindowLen,FS);
|
||||
[P2,F] = pwelch(dataAccCut_filt(end:-1:1,i),hamming(WindowLen),[],WindowLen,FS); % data
|
||||
P(1:numel(P1),i) = (P1+P2)/2;
|
||||
end
|
||||
dF = F(2)-F(1);
|
||||
dF = F(2)-F(1); % frequencies
|
||||
|
||||
% Calculate stride frequency and peak widths
|
||||
[StrideFrequency, ~, PeakWidth, MeanNormalizedPeakWidth] = StrideFrequencyRispen(P,F);
|
||||
[ResultStruct] = HarmonicityFrequency(dataAccCut_filt, P,F, StrideFrequency,dF,LowFrequentPowerThresholds,N_Harm,FS,AccVectorLen,ResultStruct);
|
||||
|
||||
ResultStruct.IH_VT = IndexOfHarmonicity(dataAccCut_filt(:,1),(5*StrideFrequency),10); % expect fmax at step frequency (2* stridefrequency), what happens if you increase fmax (~ PSD)
|
||||
ResultStruct.IH_ML = IndexOfHarmonicity(dataAccCut_filt(:,2),(5*StrideFrequency),10); % expect fmax at stride frequency
|
||||
ResultStruct.IH_AP = IndexOfHarmonicity(dataAccCut_filt(:,3),(5*StrideFrequency),10); % expect fmax at step frequency (2* stridefrequency)
|
||||
end
|
||||
124
StepcountFunc.m
124
StepcountFunc.m
@ -1,65 +1,83 @@
|
||||
function [PksAndLocsCorrected] = StepcountFunc(dataAcc_filt,StrideTimeSamples,FS);
|
||||
|
||||
% Step count funciton extracted from other function called:
|
||||
function [PksAndLocsCorrected] = StepcountFunc(dataAcc,StrideTimeSamples,FS);
|
||||
|
||||
Cutoff = 0.1;
|
||||
MinDist = floor(0.7*0.5*StrideTimeSamples); % Use StrideTimeSamples estimated above
|
||||
DataLFilt = dataAcc_filt;
|
||||
|
||||
% From acceleration to
|
||||
Vel = cumsum(detrend(DataLFilt,'constant'))/FS;
|
||||
[B,A] = butter(2,Cutoff/(FS/2),'high');
|
||||
Pos = cumsum(filtfilt(B,A,Vel))/FS;
|
||||
PosFilt = filtfilt(B,A,Pos);
|
||||
PosFiltVT = PosFilt(:,1);
|
||||
%% USE THE AP ACCELERATION: WIEBREN ZIJLSTRA: ASSESSMENT OF SPATIO-TEMPORAL PARAMTERS DURING UNCONSTRAINED WALKING (2004)
|
||||
|
||||
% Find minima and maxima in vertical position
|
||||
[PosPks,PosLocs] = findpeaks(PosFiltVT(:,1),'minpeakdistance',MinDist);
|
||||
[NegPks,NegLocs] = findpeaks(-PosFiltVT(:,1),'minpeakdistance',MinDist);
|
||||
[B,A] = butter(4,5/(FS/2),'low'); % According to article Lamoth: Multiple gait parameters derived from iPod accelerometry predict age-related gait changes
|
||||
dataAcc = filtfilt(B,A,dataAcc);
|
||||
APAcceleration = dataAcc(:,3);
|
||||
|
||||
% Find minima and maxima in AP acceleration signal
|
||||
[PosPks,PosLocs] = findpeaks(APAcceleration(:,1),'minpeakdistance',MinDist);
|
||||
[NegPks,NegLocs] = findpeaks(-APAcceleration(:,1),'minpeakdistance',MinDist);
|
||||
NegPks = -NegPks;
|
||||
if isempty(PosPks) && isempty(NegPks)
|
||||
PksAndLocs = zeros(0,3);
|
||||
else
|
||||
PksAndLocs = sortrows([PosPks,PosLocs,ones(size(PosPks)) ; NegPks,NegLocs,-ones(size(NegPks))], 2);
|
||||
end
|
||||
% Correct events for two consecutive maxima or two consecutive minima
|
||||
Events = PksAndLocs(:,2);
|
||||
NewEvents = PksAndLocs(:,2);
|
||||
Signs = PksAndLocs(:,3);
|
||||
FalseEventsIX = find(diff(Signs)==0);
|
||||
PksAndLocsToAdd = zeros(0,3);
|
||||
PksAndLocsToAddNr = 0;
|
||||
for i=1:numel(FalseEventsIX),
|
||||
FIX = FalseEventsIX(i);
|
||||
if FIX <= 2
|
||||
% remove the event
|
||||
NewEvents(FIX) = nan;
|
||||
elseif FIX >= numel(Events)-2
|
||||
% remove the next event
|
||||
NewEvents(FIX+1) = nan;
|
||||
else
|
||||
StrideTimesWhenAdding = [Events(FIX+1)-Events(FIX-2),Events(FIX+3)-Events(FIX)];
|
||||
StrideTimesWhenRemoving = Events(FIX+3)-Events(FIX-2);
|
||||
if max(abs(StrideTimesWhenAdding-StrideTimeSamples)) < abs(StrideTimesWhenRemoving-StrideTimeSamples)
|
||||
% add an event
|
||||
[M,IX] = min(Signs(FIX)*PosFiltVT((Events(FIX)+1):(Events(FIX+1)-1)));
|
||||
PksAndLocsToAddNr = PksAndLocsToAddNr+1;
|
||||
PksAndLocsToAdd(PksAndLocsToAddNr,:) = [M,Events(FIX)+IX,-Signs(FIX)];
|
||||
else
|
||||
% remove an event
|
||||
if FIX >= 5 && FIX <= numel(Events)-5
|
||||
ExpectedEvent = (Events(FIX-4)+Events(FIX+5))/2;
|
||||
else
|
||||
ExpectedEvent = (Events(FIX-2)+Events(FIX+3))/2;
|
||||
end
|
||||
if abs(Events(FIX)-ExpectedEvent) > abs(Events(FIX+1)-ExpectedEvent)
|
||||
NewEvents(FIX) = nan;
|
||||
else
|
||||
NewEvents(FIX+1) = nan;
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
PksAndLocsCorrected = sortrows([PksAndLocs(~isnan(NewEvents),:);PksAndLocsToAdd],2);
|
||||
|
||||
end
|
||||
PksAndLocsCorrected = PosLocs;
|
||||
|
||||
end
|
||||
% DataLFilt = dataAcc_filt;
|
||||
|
||||
|
||||
% % From acceleration to
|
||||
% Vel = cumsum(detrend(DataLFilt,'constant'))/FS;
|
||||
% [B,A] = butter(2,Cutoff/(FS/2),'high');
|
||||
% Pos = cumsum(filtfilt(B,A,Vel))/FS;
|
||||
% PosFilt = filtfilt(B,A,Pos);
|
||||
% PosFiltVT = PosFilt(:,1);
|
||||
|
||||
% % Find minima and maxima in vertical position
|
||||
% [PosPks,PosLocs] = findpeaks(PosFiltVT(:,1),'minpeakdistance',MinDist);
|
||||
% [NegPks,NegLocs] = findpeaks(-PosFiltVT(:,1),'minpeakdistance',MinDist);
|
||||
% NegPks = -NegPks;
|
||||
% if isempty(PosPks) && isempty(NegPks)
|
||||
% PksAndLocs = zeros(0,3);
|
||||
% else
|
||||
% PksAndLocs = sortrows([PosPks,PosLocs,ones(size(PosPks)) ; NegPks,NegLocs,-ones(size(NegPks))], 2);
|
||||
% end
|
||||
|
||||
% % Correct events for two consecutive maxima or two consecutive minima
|
||||
% Events = PksAndLocs(:,2);
|
||||
% NewEvents = PksAndLocs(:,2);
|
||||
% Signs = PksAndLocs(:,3);
|
||||
% FalseEventsIX = find(diff(Signs)==0);
|
||||
% PksAndLocsToAdd = zeros(0,3);
|
||||
% PksAndLocsToAddNr = 0;
|
||||
% for i=1:numel(FalseEventsIX),
|
||||
% FIX = FalseEventsIX(i);
|
||||
% if FIX <= 2
|
||||
% % remove the event
|
||||
% NewEvents(FIX) = nan;
|
||||
% elseif FIX >= numel(Events)-2
|
||||
% % remove the next event
|
||||
% NewEvents(FIX+1) = nan;
|
||||
% else
|
||||
% StrideTimesWhenAdding = [Events(FIX+1)-Events(FIX-2),Events(FIX+3)-Events(FIX)];
|
||||
% StrideTimesWhenRemoving = Events(FIX+3)-Events(FIX-2);
|
||||
% if max(abs(StrideTimesWhenAdding-StrideTimeSamples)) < abs(StrideTimesWhenRemoving-StrideTimeSamples)
|
||||
% % add an event
|
||||
% [M,IX] = min(Signs(FIX)*APAcceleration((Events(FIX)+1):(Events(FIX+1)-1))); % APAcceleration used to be PosFiltVT
|
||||
% PksAndLocsToAddNr = PksAndLocsToAddNr+1;
|
||||
% PksAndLocsToAdd(PksAndLocsToAddNr,:) = [M,Events(FIX)+IX,-Signs(FIX)];
|
||||
% else
|
||||
% % remove an event
|
||||
% if FIX >= 5 && FIX <= numel(Events)-5
|
||||
% ExpectedEvent = (Events(FIX-4)+Events(FIX+5))/2;
|
||||
% else
|
||||
% ExpectedEvent = (Events(FIX-2)+Events(FIX+3))/2;
|
||||
% end
|
||||
% if abs(Events(FIX)-ExpectedEvent) > abs(Events(FIX+1)-ExpectedEvent)
|
||||
% NewEvents(FIX) = nan;
|
||||
% else
|
||||
% NewEvents(FIX+1) = nan;
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% end
|
||||
% PksAndLocsCorrected = sortrows([PksAndLocs(~isnan(NewEvents),:);PksAndLocsToAdd],2);
|
||||
|
||||
@ -5,7 +5,8 @@ function [StrideFrequency, QualityInd, PeakWidth, MeanNormalizedPeakWidth] = Str
|
||||
% pwelch spectral densities
|
||||
%
|
||||
% Input:
|
||||
% AccXYZ: a three-dimensional time series with trunk accelerations
|
||||
% AccXYZ: a three-dimensional time series with trunk accelerations (is
|
||||
% for error)
|
||||
% FS: the sample frequency of the time series
|
||||
% StrideFreqGuess: a first guess of the stride frequency
|
||||
%
|
||||
@ -35,6 +36,15 @@ function [StrideFrequency, QualityInd, PeakWidth, MeanNormalizedPeakWidth] = Str
|
||||
|
||||
%% History
|
||||
% February 2013, version 1.1, adapted from StrideFrequencyFrom3dAcc
|
||||
%% Additional information
|
||||
% Stride frequency was estimated from the median of the modal
|
||||
% frequencies (half the modal frequency for VT and AP) for all three
|
||||
% directions. Whenever the median of the three frequencies fell
|
||||
% outside the range 0.6–1.2 Hz, it was replaced by one of the other
|
||||
% two frequencies, if available, within this range. Finally, if all modal
|
||||
% frequencies were within 10% of an integer multiple of the resulting
|
||||
% frequency, it was replaced by the mean of the three associated base
|
||||
% frequencies. frequencies.
|
||||
|
||||
%% Check input
|
||||
if size(AccXYZ,2) ~= 3
|
||||
@ -44,7 +54,7 @@ elseif size(AccXYZ,1) < 10*F
|
||||
end
|
||||
|
||||
|
||||
%% Get PSD
|
||||
%% Get PSD (Power Spectral Density)
|
||||
if numel(F) == 1, % Calculate the PSD from time series AccXYZ, F is the sample frequency
|
||||
AccFilt = detrend(AccXYZ,'constant'); % Detrend data to get rid of DC component in most of the specific windows
|
||||
LenPSD = 10*F;
|
||||
@ -57,13 +67,13 @@ elseif numel(F)==size(AccXYZ,1), % F are the frequencies of the power spectrum A
|
||||
Fwf = F;
|
||||
Pwf = AccXYZ;
|
||||
end
|
||||
Pwf(:,4) = sum(Pwf,2);
|
||||
Pwf(:,4) = sum(Pwf,2); % if input = F
|
||||
|
||||
|
||||
%% Estimate stride frequency
|
||||
% set parameters
|
||||
HarmNr = [2 1 2];
|
||||
CommonRange = [0.6 1.2];
|
||||
HarmNr = [2 1 2]; % AP / V --> Biphasic
|
||||
CommonRange = [0.6 1.2]; % Whenever the median of the three frequencies fell outside the range 0.6–1.2 Hz, it was replaced by one of the other two frequencies,
|
||||
% Get modal frequencies and the 'mean freq. of the peak'
|
||||
for i=1:4,
|
||||
MF1I = find([zeros(5,1);Pwf(6:end,i)]==max([zeros(5,1);Pwf(6:end,i)]),1);
|
||||
|
||||
File diff suppressed because one or more lines are too long
@ -49,7 +49,7 @@ function [SE] = funcSampleEntropy(DataIn, m, r)
|
||||
if size(DataIn,1) ~= 1 && size(DataIn,2) ~= 1
|
||||
error('DataIn must be a vector');
|
||||
end
|
||||
DataIn = DataIn(:)/std(DataIn(:));
|
||||
DataIn = DataIn(:)/std(DataIn(:)); % The data is first normalized to unit variance, rendering the outcome scale-independent.
|
||||
N = size(DataIn,1);
|
||||
if N-m <= 0
|
||||
error('m must be smaller than the length of the time series DataIn');
|
||||
@ -58,22 +58,30 @@ end
|
||||
%% Create the vectors Xm to be compared
|
||||
Xm = zeros(N-m,m);
|
||||
for i = 1:m,
|
||||
Xm(:,i) = DataIn(i:end-1-m+i,1);
|
||||
Xm(:,i) = DataIn(i:end-1-m+i,1); % (1) length = DataIn - m, so in this case [X-5,m] double
|
||||
end
|
||||
|
||||
%% Count the numbers of matches for Xm and Xmplusone
|
||||
% 1) Divide up your data into vectors of length m
|
||||
% 2) Count your like matches (Vectors are considered a match if all numbers in the vector fall within +- r (0.3)
|
||||
% 3) Calculate your conditional probability
|
||||
% 4) Divide up your data into vectors of length m+1
|
||||
% 5) Count your like matches
|
||||
% 6) Calculate your conditional probability
|
||||
% 7) Entropy is a ratio of conditional probabilities
|
||||
|
||||
CountXm = 0;
|
||||
CountXmplusone = 0;
|
||||
XmDist = nan(size(Xm));
|
||||
for i = 1:N-m,
|
||||
for j=1:m,
|
||||
for j=1:m,
|
||||
XmDist(:,j)=abs(Xm(:,j)-Xm(i,j));
|
||||
end
|
||||
IdXmi = find(max(XmDist,[],2)<=r);
|
||||
CountXm = CountXm + length(IdXmi) - 1;
|
||||
CountXmplusone = CountXmplusone + sum(abs(DataIn(IdXmi+m)-DataIn(i+m))<=r) - 1;
|
||||
IdXmi = find(max(XmDist,[],2)<=r); % (2)
|
||||
CountXm = CountXm + length(IdXmi) - 1; % (3)
|
||||
CountXmplusone = CountXmplusone + sum(abs(DataIn(IdXmi+m)-DataIn(i+m))<=r) - 1; % (6)
|
||||
end
|
||||
|
||||
%% Return sample entropy
|
||||
SE = -log(CountXmplusone/CountXm);
|
||||
SE = -log(CountXmplusone/CountXm); % (7) A srictly periodic time-series is completely predictable and will have a SEn of zero. SEn is defined as the negative natural logarithm of an estimate of the condiotional probability of epochs of length m (m=5) that match point-wise within a tolerance r and repeates itself for m+1 points.
|
||||
|
||||
|
||||
104
matlab.sty
104
matlab.sty
@ -1,104 +0,0 @@
|
||||
\NeedsTeXFormat{LaTeX2e}
|
||||
\ProvidesPackage{matlab}
|
||||
|
||||
\RequirePackage{verbatim}
|
||||
\RequirePackage{fancyvrb}
|
||||
\RequirePackage{alltt}
|
||||
\RequirePackage{upquote}
|
||||
\RequirePackage[framemethod=tikz]{mdframed}
|
||||
\RequirePackage{hyperref}
|
||||
\RequirePackage{color}
|
||||
|
||||
|
||||
\newcommand{\maxwidth}[1]{\ifdim\linewidth>#1 #1\else\linewidth\fi}
|
||||
\newcommand{\mlcell}[1]{{\color{output}\verbatim@font#1}}
|
||||
|
||||
\definecolor{output}{gray}{0.4}
|
||||
|
||||
% Unicode character conversions
|
||||
\DeclareUnicodeCharacter{B0}{\ensuremath{^\circ}}
|
||||
\DeclareUnicodeCharacter{21B5}{\ensuremath{\hookleftarrow}}
|
||||
|
||||
% Paragraph indentation
|
||||
\setlength{\parindent}{0pt}
|
||||
|
||||
% Hyperlink style
|
||||
\hypersetup{
|
||||
colorlinks=true,
|
||||
linkcolor=blue,
|
||||
urlcolor=blue
|
||||
}
|
||||
|
||||
|
||||
% environment styles for MATLAB code and output
|
||||
\mdfdefinestyle{matlabcode}{%
|
||||
outerlinewidth=.5pt,
|
||||
linecolor=gray!20!white,
|
||||
roundcorner=2pt,
|
||||
innertopmargin=.5\baselineskip,
|
||||
innerbottommargin=.5\baselineskip,
|
||||
innerleftmargin=1em,
|
||||
backgroundcolor=gray!10!white
|
||||
}
|
||||
|
||||
\newenvironment{matlabcode}{\verbatim}{\endverbatim}
|
||||
\surroundwithmdframed[style=matlabcode]{matlabcode}
|
||||
|
||||
\newenvironment{matlaboutput}{%
|
||||
\Verbatim[xleftmargin=1.25em, formatcom=\color{output}]%
|
||||
}{\endVerbatim}
|
||||
|
||||
\newenvironment{matlabsymbolicoutput}{%
|
||||
\list{}{\leftmargin=1.25em\relax}%
|
||||
\item\relax%
|
||||
\color{output}\verbatim@font%
|
||||
}{\endlist}
|
||||
|
||||
\newenvironment{matlabtableoutput}[1]{%
|
||||
{\color{output}%
|
||||
\hspace*{1.25em}#1}%
|
||||
}{}
|
||||
|
||||
|
||||
% Table of Contents style
|
||||
\newcounter{multititle}
|
||||
\newcommand{\matlabmultipletitles}{\setcounter{multititle}{1}}
|
||||
|
||||
\newcounter{hastoc}
|
||||
\newcommand{\matlabhastoc}{\setcounter{hastoc}{1}}
|
||||
|
||||
\newcommand{\matlabtitle}[1]{
|
||||
\ifnum\value{multititle}>0
|
||||
\ifnum\value{hastoc}>0
|
||||
\addcontentsline{toc}{section}{#1}
|
||||
\fi
|
||||
\fi
|
||||
\section*{#1}
|
||||
}
|
||||
|
||||
\newcommand{\matlabheading}[1]{
|
||||
\ifnum\value{hastoc}>0
|
||||
\addcontentsline{toc}{subsection}{#1}
|
||||
\fi
|
||||
\subsection*{#1}
|
||||
}
|
||||
|
||||
\newcommand{\matlabheadingtwo}[1]{
|
||||
\ifnum\value{hastoc}>0
|
||||
\addcontentsline{toc}{subsubsection}{#1}
|
||||
\fi
|
||||
\subsubsection*{#1}
|
||||
}
|
||||
|
||||
\newcommand{\matlabheadingthree}[1]{
|
||||
\ifnum\value{hastoc}>0
|
||||
\addcontentsline{toc}{paragraph}{#1}
|
||||
\fi
|
||||
\paragraph*{#1}
|
||||
}
|
||||
|
||||
\newcommand{\matlabtableofcontents}[1]{
|
||||
\renewcommand{\contentsname}{#1}
|
||||
\tableofcontents
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user