diff --git a/ControlsData.mat b/ControlsData.mat new file mode 100644 index 0000000..bbe2179 Binary files /dev/null and b/ControlsData.mat differ diff --git a/ControlsDetermineLocationTurnsFunc.m b/ControlsDetermineLocationTurnsFunc.m new file mode 100644 index 0000000..f2b38fa --- /dev/null +++ b/ControlsDetermineLocationTurnsFunc.m @@ -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 \ No newline at end of file diff --git a/FINAL_MAP_CLBP_Controls_23052021.xls b/FINAL_MAP_CLBP_Controls_23052021.xls new file mode 100644 index 0000000..47c1dd3 Binary files /dev/null and b/FINAL_MAP_CLBP_Controls_23052021.xls differ diff --git a/FINAL_NPRSBaseline_Minute_23052021.xls b/FINAL_NPRSBaseline_Minute_23052021.xls new file mode 100644 index 0000000..ad9c983 Binary files /dev/null and b/FINAL_NPRSBaseline_Minute_23052021.xls differ diff --git a/FINAL_NPRS_GaitCharacteristics_23052021.xls b/FINAL_NPRS_GaitCharacteristics_23052021.xls new file mode 100644 index 0000000..d78e92c Binary files /dev/null and b/FINAL_NPRS_GaitCharacteristics_23052021.xls differ diff --git a/FalseNearestNeighborsSR.m b/FalseNearestNeighborsSR.m new file mode 100644 index 0000000..7e2d4f8 --- /dev/null +++ b/FalseNearestNeighborsSR.m @@ -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. +%======================================================================== +% , 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 +% + +%======================================================================== +% 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 diff --git a/GaitOutcomesTrunkAccFuncLD.m b/GaitOutcomesTrunkAccFuncLD.m new file mode 100644 index 0000000..c1141cf --- /dev/null +++ b/GaitOutcomesTrunkAccFuncLD.m @@ -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 \ No newline at end of file diff --git a/GaitVarOutcomesCombined_230221.mat b/GaitVarOutcomesCombined_230221.mat new file mode 100644 index 0000000..92626ed Binary files /dev/null and b/GaitVarOutcomesCombined_230221.mat differ diff --git a/Hratio.m b/Hratio.m new file mode 100644 index 0000000..ba6c435 --- /dev/null +++ b/Hratio.m @@ -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)); diff --git a/IndexOfHarmonicity.m b/IndexOfHarmonicity.m new file mode 100644 index 0000000..b42933c --- /dev/null +++ b/IndexOfHarmonicity.m @@ -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); diff --git a/WindMainFreq.m b/WindMainFreq.m new file mode 100644 index 0000000..0aa398d --- /dev/null +++ b/WindMainFreq.m @@ -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 & f0 + 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)) 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 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 \ No newline at end of file diff --git a/struct2csv - Snelkoppeling.lnk b/struct2csv - Snelkoppeling.lnk new file mode 100644 index 0000000..9227850 Binary files /dev/null and b/struct2csv - Snelkoppeling.lnk differ