s3158985
/
MAP_Gait_Dynamics
Template
1
0
Fork 0
Code Issues Pull Requests 1 Releases Wiki Activity

add missing .m files for MAIN

This commit is contained in:
L. Dijk 2021-07-03 16:37:53 +02:00
parent 0be926128f
commit c64768e2bd
16 changed files with 1085 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
ControlsData.mat Normal file
View File

Binary file not shown.

91
ControlsDetermineLocationTurnsFunc.m Normal file
View File

@ -0,0 +1,91 @@
function [locsTurns,FilteredData,FootContacts] = ControlsDetermineLocationTurnsFunc(inputData,FS,ApplyRealignment,plotit,ResultStruct)
% Description: Determine the location of turns, plot for visual inspection
% Input: Acc Data (not yet realigned)
%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).
data = inputData(:,[3,2,1]); % NOT THE SAME AS IN THE CLBP DATA; 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
data = inputData; % Realignment performed in previous stage
dataAcc = inputData;
end
% Filter data
[B,A] = butter(2,20/(FS/2),'low'); % Filters data very strongly which is needed to determine turns correctly
dataStepDetection = filtfilt(B,A,dataAcc);
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; % (1:2:end,2); previous version
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*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.3*FS); % for step detection
numStepsOption2_filt = numel(pks); % counts number of locs for turn detection
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',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',0.2*FS); % values were initially 5
locsTurns = [locs(locs_diffLocs), locs(locs_diffLocs+1)];
%magNoGfilt_copy = magNoGfilt;
VTAcc_copy = data(:,1);
APAcc_copy = data(:,3);
for k = 1: size(locsTurns,1);
VTAcc_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
APAcc_copy(locsTurns(k,1):locsTurns(k,2)) = NaN;
end
% Visualising signal;
if plotit
figure()
subplot(2,1,1)
hold on;
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],[5,5,5,5,5,5;10,10,10,10,10,10],'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],[0,0,0,0,0,0;10,10,10,10,10,10],'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))])
%disp(['Total distance measured by researcher was = ' num2str(Distance)])
FilteredData = dataAcc;
end

BIN
FINAL_MAP_CLBP_Controls_23052021.xls Normal file
View File

Binary file not shown.

BIN
FINAL_NPRSBaseline_Minute_23052021.xls Normal file
View File

Binary file not shown.

BIN
FINAL_NPRS_GaitCharacteristics_23052021.xls Normal file
View File

Binary file not shown.

126
FalseNearestNeighborsSR.m Normal file
View File

@ -0,0 +1,126 @@
function fnnM = FalseNearestNeighborsSR(xV,tauV,mV,escape,theiler)
% fnnM = FalseNearestNeighbors(xV,tauV,mV,escape,theiler)
% FALSENEARESTNEIGHBORS computes the percentage of false nearest neighbors
% for a range of delays in 'tauV' and embedding dimensions in 'mV'.
% INPUT
% xV : Vector of the scalar time series
% tauV : A vector of the delay times.
% mV : A vector of the embedding dimension.
% escape : A factor of escaping from the neighborhood. Default=10.
% theiler : the Theiler window to exclude time correlated points in the
% search for neighboring points. Default=0.
% OUTPUT:
% fnnM : A matrix of size 'ntau' x 'nm', where 'ntau' is the number of
% given delays and 'nm' is the number of given embedding
% dimensions, containing the percentage of false nearest
% neighbors.
%========================================================================
% <FalseNearestNeighbors.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/
%=========================================================================
% Updates:
% Sietse Rispens, January 2012: use a different kdtree algorithm,
% searching k nearest neighbours, to improve performance
%
%=========================================================================
fthres = 0.1; % A factor of the data SD to be used as the maximum radius
% for searching for valid nearest neighbor.
propthres = fthres; % Limit for the proportion of valid points, i.e. points
% for which a nearest neighbor was found. If the proporion
% of valid points is beyond this limit, do not compute
% FNN.
n = length(xV);
if isempty(escape), escape=10; end
if isempty(theiler), theiler=0; end
% Rescale to [0,1] and add infinitesimal noise to have distinct samples
xmin = min(xV);
xmax = max(xV);
xV = (xV - xmin) / (xmax-xmin);
xV = AddNoise(xV,10^(-10));
ntau = length(tauV);
nm = length(mV);
fnnM = NaN*ones(ntau,nm);
for itau = 1:ntau
tau = tauV(itau);
for im=1:nm
m = mV(im);
nvec = n-m*tau; % to be able to add the component x(nvec+tau) for m+1
if nvec-theiler < 2
break;
end
xM = NaN*ones(nvec,m);
for i=1:m
xM(:,m-i+1) = xV(1+(i-1)*tau:nvec+(i-1)*tau);
end
% k-d-tree data structure of the training set for the given m
TreeRoot=kdtree_build(xM);
% For each target point, find the nearest neighbor, and check whether
% the distance increase over the escape distance by adding the next
% component for m+1.
idxV = NaN*ones(nvec,1);
distV = NaN*ones(nvec,1);
k0 = 2; % The initial number of nearest neighbors to look for
kmax = min(2*theiler + 2,nvec); % The maximum number of nearest neighbors to look for
for i=1:nvec
tarV = xM(i,:);
iV = [];
k=k0;
kmaxreached = 0;
while isempty(iV) && ~kmaxreached
[neiindV, neidisV] = kdtree_k_nearest_neighbors(TreeRoot,tarV,k);
% [neiM,neidisV,neiindV]=kdrangequery(TreeRoot,tarV,rthres*sqrt(m));
[oneidisV,oneiindV]=sort(neidisV);
neiindV = neiindV(oneiindV);
neidisV = neidisV(oneiindV);
iV = find(abs(neiindV(1)-neiindV(2:end))>theiler);
if k >= kmax
kmaxreached = 1;
elseif isempty(iV)
k = min(kmax,k*2);
end
end
if ~isempty(iV)
idxV(i) = neiindV(iV(1)+1);
distV(i) = neidisV(iV(1)+1);
end
end % for i
iV = find(~isnan(idxV));
nproper = length(iV);
% Compute fnn only if there is a sufficient number of target points
% having nearest neighbor (in R^m) within the threshold distance
if nproper>propthres*nvec
nnfactorV = 1+(xV(iV+m*tau)-xV(idxV(iV)+m*tau)).^2./distV(iV).^2;
fnnM(itau,im) = length(find(nnfactorV > escape^2))/nproper;
end
kdtree_delete(TreeRoot); % Free the pointer to k-d-tree
end
end

210
GaitOutcomesTrunkAccFuncLD.m Normal file
View File

@ -0,0 +1,210 @@
function [ResultStruct] = GaitOutcomesTrunkAccFuncLD(inputData,FS,LegLength,WindowLen,ApplyRealignment,ApplyRemoveSteps)
% Description: 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 1 step in the beginning and end of data (comment Bert Otten)
dataAccCut = dataAcc(LocsSteps(1):LocsSteps(end-1),:);
dataAccCut_filt = dataAcc_filt(LocsSteps(1):LocsSteps(end-1),:);
% 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

BIN
GaitVarOutcomesCombined_230221.mat Normal file
View File

Binary file not shown.

55
Hratio.m Normal file
View File

@ -0,0 +1,55 @@
function [Hratio,wb]=Hratio(x,fmax,numfreq)
% function [HIndex]=IndexOfHarmonicity(x,numfreq);
%
% INPUT, signal, fmax is max frequency that can be considered as main
% frequency = around 1 for stride and around 2 for steps.
% numfreq = number of superharmonics in analysis
% OUTPUT: Harmoncity Index and mainfrequency
% PLOT: with no output a plot is generated
% CALLS: WndMainFreq.m
% CL 2016
if nargin<3 || isempty(numfreq), numfreq = 10; end
fs = 100;% sample frequency
le=length(x);
temp=[];
n=10*length(x);
ws=fix(length(x)/2);
no=fix(0.99*ws);
[Px,fx]=spectrum2(x,n,no,boxcar(le),100,'mean');
Px = Px/sum(Px); % All power spectral densities were normalized by dividing the power by the sum of the total power spectrum, which equals the variance.
% select frequency
wb = WindMainFreq(x,100,fmax) ; % mainfrequency frequency
wbr= wb + 0.055; % % basis mainfrequency + tail
wbl = wb - 0.055;
wind=find(fx >=wbl & fx<= wbr); % index van basisfrequentie
PI_0=sum(Px(wind));
temp=PI_0;
for k=1:numfreq;
w=wb*(k+1);
wr=w+0.055;
wl=w-0.055;
wInd=find(fx >=wl & fx<=wr);
if nargout ==0; % controle
plot(fx,Px,'r'); hold on;
set(gca','XLim', [0 6],'FontSize', 8)
xlabel('Frequency (Hz)'); ylabel('Power');
plot(wb,max(Px),'r.','MarkerSize', 10),
plot(wbr,0,'c.'); plot(wbl,0,'c.');
plot(w,0,'r.','MarkerSize', 10)
plot(wl,0,'c.'); plot(wr,0,'c.');
end
PI=sum(Px(wInd));
temp=[temp; PI];
end
Hratio = sum(temp(1:2:end-1))/sum(temp(2:2:end));

54
IndexOfHarmonicity.m Normal file
View File

@ -0,0 +1,54 @@
function [HIndex,wb]=IndexOfHarmonicity(x,fmax,numfreq)
% function [HIndex]=IndexOfHarmonicity(x,numfreq);
%
% INPUT, signal, fmax is max frequency that can be considered as main
% frequency = around 1 for stride and around 2 for steps.
% numfreq = number of superharmonics in analysis
% OUTPUT: Harmoncity Index and mainfrequency
% PLOT: with no output a plot is generated
% CALLS: WndMainFreq.m
% CL 2016
if nargin<3 || isempty(numfreq), numfreq = 10; end
fs = 100;% sample frequency
le=length(x);
temp=[];
n=10*length(x);
ws=fix(length(x)/2);
no=fix(0.99*ws);
[Px,fx]=spectrum2(x,n,no,boxcar(le),100,'mean');
Px = Px/sum(Px); % All power spectral densities were normalized by dividing the power by the sum of the total power spectrum, which equals the variance.
% select frequency
wb = WindMainFreq(x,100,fmax) ; % mainfrequency frequency
wbr= wb + 0.055; % % basis mainfrequency + tail
wbl = wb - 0.055;
wind=find(fx >=wbl & fx<= wbr); % index van basisfrequentie
PI_0=sum(Px(wind));
temp=PI_0;
for k=1:numfreq;
w=wb*(k+1);
wr=w+0.055;
wl=w-0.055;
wInd=find(fx >=wl & fx<=wr);
if nargout ==0; % controle
plot(fx,Px,'r'); hold on;
set(gca','XLim', [0 6],'FontSize', 8)
xlabel('Frequency (Hz)'); ylabel('Power');
plot(wb,max(Px),'r.','MarkerSize', 10),
plot(wbr,0,'c.'); plot(wbl,0,'c.');
plot(w,0,'r.','MarkerSize', 10)
plot(wl,0,'c.'); plot(wr,0,'c.');
end
PI=sum(Px(wInd));
temp=[temp; PI];
end
HIndex=PI_0/sum(temp);

26
WindMainFreq.m Normal file
View File

@ -0,0 +1,26 @@
function Omega=WindMainFreq(x,fs,fmax);
% function Omega=WindMainFreq(x,fs,fmax);
%INPUT: x = signal; fs = sample frequency; fmax = max frequency for
%mainfrequency
% OUTPUT = mainfrequency Omega
% CALSS: spectrum2.m
% PLOT : with no outtput [] generates plot
if nargin<3 || isempty(fmax), fmax = 4; end
n=10*length(x);
ws=fix(length(x)/2);
no=fix(0.99*ws);
[p,f]=spectrum2(x,n,no,hamming(ws),fs,'mean');
p=p/sum(p);
aa=find(f>0.5 & f<fmax);
[a,b]=max(p(aa));
Omega=f(b+aa(1)-1);
if nargout ==0;
ii=find(f(:,1)<=fmax);
f=f(1:length(ii),:);
p=p(1:length(ii),:);
plot(f,p,'b');hold on;
plot(Omega,max(p(aa)),'m.','MarkerSize', 15); %check w0
end

68
cross_sampen.m Normal file
View File

@ -0,0 +1,68 @@
function [e,A,B]=cross_sampen(x,y,M,r,sflag);
%function [e,A,B]=cross_sampen(x,y,M,r);
%
% The Cross-sample Entropy (Cross-SampEn) quantified the degree of synchronization between AP and ML, AP and V, and ML and V accelerations.
% Cross-SampEn is the negative natural logarithm of the conditional probability that epochs with length m that match point-wise
% in the two related signals, repeat itself for m+1 points, within a tolerance of r (in the present study m = 5 and r = 0.3)
%Input
%
%x,y input data
%M maximum template length
%r matching tolerance
%sflag flag to standardize signals(default yes/sflag=1)
%
%Output
%
%e sample entropy estimates for m=0,1,...,M-1
%A number of matches for m=1,...,M
%B number of matches for m=0,...,M-1 excluding last point
if ~exist('sflag')|isempty(sflag),sflag=1;end
y=y(:);
x=x(:);
ny=length(y);
nx=length(x);
N=nx;
if sflag>0
y=y-mean(y);
sy=sqrt(mean(y.^2));
y=y/sy;
x=x-mean(x);
sx=sqrt(mean(x.^2));
x=x/sx;
end
lastrun=zeros(nx,1);
run=zeros(nx,1);
A=zeros(M,1);
B=zeros(M,1);
p=zeros(M,1);
e=zeros(M,1);
for i=1:ny
for j=1:nx
if abs(x(j)-y(i))<r
run(j)=lastrun(j)+1;
M1=min(M,run(j));
for m=1:M1
A(m)=A(m)+1;
if (i<ny)&(j<nx)
B(m)=B(m)+1;
end
end
else
run(j)=0;
end
end
for j=1:nx
lastrun(j)=run(j);
end
end
N=ny*nx;
B=[N;B(1:(M-1))];
p=A./B;
e=-log(p);

BIN
delta_NRPS_Gaitoutcomes.xlsx Normal file
View File

Binary file not shown.

28
msentropy.m Normal file
View File

@ -0,0 +1,28 @@
% Function for calculating multiscale entropy
% input: signal
% m: match point(s)
% r: matching tolerance
% factor: number of scale factor
% sampenc is available at http://people.ece.cornell.edu/land/PROJECTS/Complexity/sampenc.m
%
% Multi-scale sample Entropy (Mscale-En) is an indicator of gait predictability. Multi-scale entropy takes the
% complexity of a system into account by calculating the
% predictability of a signal over time scales with increasing
% length. A ‘coarse-graining’ process is applied to the acceleration signals; non-overlapping windows of data
% points with an increasing length τ are constructed, with
% τ representing the time scale with a tolerance of r (in the
% present study τ = 7 and r = 0.2). A complete predictable
% signal will adopt a Mscale-En value of 0 [29].
function e = msentropy(input,m,r,factor)
y=input;
y=y-mean(y);
y=y/std(y);
for i=1:factor
s=coarsegraining(y,i); % dont have this function yet.
sampe=sampenc(s,m+1,r);
e(i)=sampe(m+1);
end
e=e';

427
spectrum2.m Normal file
View File

@ -0,0 +1,427 @@
function [Spec,f,SpecConf] = spectrum2(varargin)
%SPECTRUM Power spectrum estimate of one or two data sequences.
% P=SPECTRUM2(X,NFFT,NOVERLAP,WIND) estimates the Power Spectral Density of
% signal vector(s) X using Welch's averaged periodogram method. X must be a
% set of column vectors. They are divided into overlapping sections, each of
% which is detrended and windowed by the WINDOW parameter, then zero padded
% to length NFFT. The magnitude squared of the length NFFT DFTs of the sec-
% tions are averaged to form Pxx. P is a two times SIZE(X,2) column matrix
% P = [Pxx Pxxc]; the second half of columns Pxxc are the 95% confidence
% intervals. The number of rows of P is NFFT/2+1 for NFFT even, (NFFT+1)/2
% for NFFT odd, or NFFT if the signal X is complex. If you specify a scalar
% for WINDOW, a Hanning window of that length is used.
%
% [P,F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs) given a sampling frequency
% Fs returns a vector of frequencies the same length as Pxx at which the
% PSD is estimated. PLOT(F,P(:,1:end/2)) plots the power spectrum estimate
% versus true frequency.
%
% [P, F] = SPECTRUM2(X,NFFT,NOVERLAP,WINDOW,Fs,Pr) where Pr is a scalar
% between 0 and 1, overrides the default 95% confidence interval and
% returns the Pr*100% confidence interval for Pxx instead.
%
% SPECTRUM2(X) with no output arguments plots the PSD in the current
% figure window, with confidence intervals.
%
% The default values for the parameters are NFFT = 256 (or SIZE(X,1),
% whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2,
% and Pr = .95. You can obtain a default parameter by leaving it out
% or inserting an empty matrix [], e.g. SPECTRUM(X,[],128).
%
% P = SPECTRUM2(X,Y) performs spectral analysis of the two times SIZE(X,2)
% sequences X and Y using the Welch method. SPECTRUM returns the 8 times
% SIZE(X,2)column array
% P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]
% where
% Pxx = X-vector power spectral density
% Pyy = Y-vector power spectral density
% Pxy = Cross spectral density
% Txy = Complex transfer function from X to Y = Pxy./Pxx
% Cxy = Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
% Pxxc,Pyyc,Pxyc = Confidence range.
% All input and output options are otherwise exactly the same as for the
% single input case.
%
% SPECTRUM2(X,Y) with no output arguments will plot Pxx, Pyy, abs(Txy),
% angle(Txy) and Cxy in sequence, pausing between plots.
%
% SPECTRUM2(X,...,DFLAG), where DFLAG can be 'linear', 'mean' or 'none',
% specifies a detrending mode for the prewindowed sections of X (and Y).
% DFLAG can take the place of any parameter in the parameter list
% (besides X) as long as it is last, e.g. SPECTRUM(X,'none');
%
% See also SPECTRUM, PSD, CSD, TFE, COHERE, SPECGRAM, SPECPLOT, DETREND,
% PMTM, PMUSIC.
% ETFE, SPA, and ARX in the Identification Toolbox.
% The units on the power spectra Pxx and Pyy are such that, using
% Parseval's theorem:
%
% SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/size(x,1) = COV(X)
%
% The RMS value of the signal is the square root of this.
% If the input signal is in Volts as a function of time, then
% the units on Pxx are Volts^2*seconds = Volt^2/Hz.
%
% Here are the covariance, RMS, and spectral amplitude values of
% some common functions:
% Function Cov=SUM(Pxx)/LENGTH(Pxx) RMS Pxx
% a*sin(w*t) a^2/2 a/sqrt(2) a^2*LENGTH(Pxx)/4
%Normal: a*rand(t) a^2 a a^2
%Uniform: a*rand(t) a^2/12 a/sqrt(12) a^2/12
%
% For example, a pure sine wave with amplitude A has an RMS value
% of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)).
%
% See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal
% Processing, Prentice-Hall, 1975.
error(nargchk(1,8,nargin))
[msg,x,y,nfft,noverlap,window,Fs,p,dflag]=specchk2(varargin);
error(msg)
if isempty(p),
p = .95; % default confidence interval even if not asked for
end
n = size(x,1); % Number of data points
ns = size(x,2); % Number of signals
nwind = length(window);
if n < nwind % zero-pad x (and y) if length less than the window length
x(nwind,:)=0; n=nwind;
if ~isempty(y), y(nwind)=0; end
end
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
index = 1:nwind;
KMU = k*norm(window)^2; % Normalizing scale factor ==> asymptotically unbiased
% KMU = k*sum(window)^2;% alt. Nrmlzng scale factor ==> peaks are about right
if (isempty(y)) % Single sequence case.
Pxx = zeros(nfft,ns); Pxx2 = zeros(nfft,ns);
for i=1:k
for l=1:ns
if strcmp(dflag,'linear')
xw = window.*detrend(x(index,l));
elseif strcmp(dflag,'none')
xw = window.*(x(index,l));
else
xw = window.*detrend(x(index,l),0);
end
Xx = abs(fft(xw,nfft)).^2;
Pxx(:,l) = Pxx(:,l) + Xx;
Pxx2(:,l) = Pxx2(:,l) + abs(Xx).^2;
end
index = index + (nwind - noverlap);
end
% Select first half
if ~any(any(imag(x)~=0)), % if x and y are not complex
if rem(nfft,2), % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
else
select = 1:nfft;
end
Pxx = Pxx(select,:);
Pxx2 = Pxx2(select,:);
cPxx = zeros(size(Pxx));
if k > 1
c = (k.*Pxx2-abs(Pxx).^2)./(k-1);
c = max(c,zeros(size(Pxx)));
cPxx = sqrt(c);
end
ff = sqrt(2)*erfinv(p); % Equal-tails.
Pxxc = ff.*cPxx/KMU;
P = Pxx/KMU;
Pc = Pxxc;
else
Pxx = zeros(nfft,ns); % Dual sequence case.
Pxy = Pxx; Pxx2 = Pxx; Pxy2 = Pxx;
Pyy = zeros(nfft,1); Pyy2 = Pyy;
for i=1:k
if strcmp(dflag,'linear')
yw = window.*detrend(y(index));
elseif strcmp(dflag,'none')
yw = window.*(y(index));
else
yw = window.*detrend(y(index),0);
end
Yy = fft(yw,nfft);
Yy2 = abs(Yy).^2;
Pyy = Pyy + Yy2;
Pyy2 = Pyy2 + abs(Yy2).^2;
for l=1:ns
if strcmp(dflag,'linear')
xw = window.*detrend(x(index,l));
elseif strcmp(dflag,'none')
xw = window.*(x(index,l));
else
xw = window.*detrend(x(index,l),0);
end
Xx = fft(xw,nfft);
Xx2 = abs(Xx).^2;
Pxx(:,l) = Pxx(:,l) + Xx2;
Pxx2(:,l) = Pxx2(:,l) + abs(Xx2).^2;
Xy = Yy .* conj(Xx);
Pxy(:,l) = Pxy(:,l) + Xy;
Pxy2(:,l) = Pxy2(:,l) + Xy .* conj(Xy);
end
index = index + (nwind - noverlap);
end
% Select first half
if ~any(any(imag([x y])~=0)), % if x and y are not complex
if rem(nfft,2), % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
else
select = 1:nfft;
end
Pxx = Pxx(select,:);
Pxy = Pxy(select,:);
Pxx2 = Pxx2(select,:);
Pxy2 = Pxy2(select,:);
Pyy = Pyy(select);
Pyy2 = Pyy2(select);
cPxx = zeros(size(Pxx));
cPyy = zeros(size(Pyy));
cPxy = cPxx;
if k > 1
c = max((k.*Pxx2-abs(Pxx).^2)./(k-1),zeros(size(Pxx)));
cPxx = sqrt(c);
c = max((k.*Pyy2-abs(Pyy).^2)./(k-1),zeros(size(Pyy)));
cPyy = sqrt(c);
c = max((k.*Pxy2-abs(Pxy).^2)./(k-1),zeros(size(Pxx)));
cPxy = sqrt(c);
end
Txy = Pxy./Pxx;
Cxy = (abs(Pxy).^2)./(Pxx.*repmat(Pyy,1,size(Pxx,2)));
ff = sqrt(2)*erfinv(p); % Equal-tails.
Pxx = Pxx/KMU;
Pyy = Pyy/KMU;
Pxy = Pxy/KMU;
Pxxc = ff.*cPxx/KMU;
Pxyc = ff.*cPxy/KMU;
Pyyc = ff.*cPyy/KMU;
P = [Pxx Pyy Pxy Txy Cxy];
Pc = [Pxxc Pyyc Pxyc];
end
freq_vector = (select - 1)'*Fs/nfft;
if nargout == 0, % do plots
newplot;
if Fs==2, xl='Frequency'; else, xl = 'f [Hz]'; end
nplot=1+(1-isempty(y))*4;
subplot(nplot,1,1);
c = [max(Pxx-Pxxc,0) Pxx+Pxxc];
c = c.*(c>0);
h=semilogy(freq_vector,Pxx,...
freq_vector,c(:,1:size(c,2)/2),'--',...
freq_vector,c(:,size(c,2)/2+1:end),'--');
title('\bf X Power Spectral Density')
ylabel('P_x')
xlabel(xl)
if length(h)>3
s={};
for k=1:length(h)/3
c=get(h(k),'Color');
set(h(k+length(h)/3),'Color',c);
set(h(k+2*length(h)/3),'Color',c);
s{k}=['x_' num2str(k)];
end
legend(h(1:length(h)/3),s);
end
if (isempty(y)), % single sequence case
return
end
subplot(nplot,1,2);
c = [max(Pyy-Pyyc,0) Pyy+Pyyc];
c = c.*(c>0);
h=semilogy(freq_vector,Pyy,...
freq_vector,c(:,1),'--',...
freq_vector,c(:,2),'--');
if size(Pxx,2)>1
for k=1:length(h)/3
c=get(h(k),'Color');
set(h(k+length(h)/3),'Color',c);
set(h(k+2*length(h)/3),'Color',c);
end
end
title('\bf Y Power Spectral Density')
ylabel('P_y')
xlabel(xl)
subplot(nplot,1,3);
semilogy(freq_vector,abs(Txy));
title('\bf Transfer function magnitude')
ylabel('T_{xy}^{(m)}')
xlabel(xl)
subplot(nplot,1,4);
plot(freq_vector,(angle(Txy))), ...
title('\bf Transfer function phase')
ylabel('T_{xy}^{(p)}')
xlabel(xl)
subplot(nplot,1,5);
plot(freq_vector,Cxy);
title('\bf Coherence')
ylabel('C_{xy}')
xlabel(xl)
if exist ('niceaxes') ~= 0
niceaxes(findall(gcf,'Type','axes'));
end
elseif nargout ==1,
Spec = P;
elseif nargout ==2,
Spec = P;
f = freq_vector;
elseif nargout ==3,
Spec = P;
f = freq_vector;
SpecConf = Pc;
end
function [msg,x,y,nfft,noverlap,window,Fs,p,dflag] = specchk2(P)
%SPECCHK Helper function for SPECTRUM
% SPECCHK(P) takes the cell array P and uses each cell as
% an input argument. Assumes P has between 1 and 7 elements.
msg = [];
if size(P{1},1)<=1
if max(size(P{1}))==1
msg = 'Input data must be a vector, not a scalar.';
else
msg = 'Requires column vector input.';
end
x=[];
y=[];
elseif (length(P)>1),
if (all(size(P{1},1)==size(P{2})) & (size(P{1},1)>1) ) | ...
size(P{2},1)>1, % 0ne signal or 2 present?
% two signals, x and y, present
x = P{1}; y = P{2};
% shift parameters one left
P(1) = [];
else
% only one signal, x, present
x = P{1}; y = [];
end
else % length(P) == 1
% only one signal, x, present
x = P{1}; y = [];
end
% now x and y are defined; let's get the rest
if length(P) == 1
nfft = min(size(x,1),256);
window = hanning(nfft);
noverlap = 0;
Fs = 2;
p = [];
dflag = 'linear';
elseif length(P) == 2
if isempty(P{2}), dflag = 'linear'; nfft = min(size(x,1),256);
elseif isstr(P{2}), dflag = P{2}; nfft = min(size(x,1),256);
else dflag = 'linear'; nfft = P{2}; end
window = hanning(nfft);
noverlap = 0;
Fs = 2;
p = [];
elseif length(P) == 3
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
if isempty(P{3}), dflag = 'linear'; noverlap = 0;
elseif isstr(P{3}), dflag = P{3}; noverlap = 0;
else dflag = 'linear'; noverlap = P{3}; end
window = hanning(nfft);
Fs = 2;
p = [];
elseif length(P) == 4
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
if isstr(P{4})
dflag = P{4};
window = hanning(nfft);
else
dflag = 'linear';
window = P{4}; window = window(:); % force window to be a column
if length(window) == 1, window = hanning(window); end
if isempty(window), window = hanning(nfft); end
end
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
Fs = 2;
p = [];
elseif length(P) == 5
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
window = P{4}; window = window(:); % force window to be a column
if length(window) == 1, window = hanning(window); end
if isempty(window), window = hanning(nfft); end
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
if isstr(P{5})
dflag = P{5};
Fs = 2;
else
dflag = 'linear';
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
end
p = [];
elseif length(P) == 6
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
window = P{4}; window = window(:); % force window to be a column
if length(window) == 1, window = hanning(window); end
if isempty(window), window = hanning(nfft); end
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
if isstr(P{6})
dflag = P{6};
p = [];
else
dflag = 'linear';
if isempty(P{6}), p = .95; else p = P{6}; end
end
elseif length(P) == 7
if isempty(P{2}), nfft = min(size(x,1),256); else nfft=P{2}; end
window = P{4}; window = window(:); % force window to be a column
if length(window) == 1, window = hanning(window); end
if isempty(window), window = hanning(nfft); end
if isempty(P{3}), noverlap = 0; else noverlap=P{3}; end
if isempty(P{5}), Fs = 2; else Fs = P{5}; end
if isempty(P{6}), p = .95; else p = P{6}; end
if isstr(P{7})
dflag = P{7};
else
msg = 'DFLAG parameter must be a string.'; return
end
end
% NOW do error checking
if isempty(msg)
if (nfft<length(window)),
msg = 'Requires window''s length to be no greater than the FFT length.';
end
if (noverlap >= length(window)),
msg = 'Requires NOVERLAP to be strictly less than the window length.';
end
if (nfft ~= abs(round(nfft)))|(noverlap ~= abs(round(noverlap))),
msg = 'Requires positive integer values for NFFT and NOVERLAP.';
end
if ~isempty(p),
if (prod(size(p))>1)|(p(1,1)>1)|(p(1,1)<0),
msg = 'Requires confidence parameter to be a scalar between 0 and 1.';
end
end
if (min(size(y))~=1)&(~isempty(y)),
msg = 'Requires column vector input as second signal.';
end
if (size(x,1)~=length(y))&(~isempty(y)),
msg = 'Requires X and Y to have the same number of rows.';
end
end

BIN
struct2csv - Snelkoppeling.lnk Normal file
View File

Binary file not shown.