reads in dryft sample grf data, filters it, splits steps, normalizes to stance phase, then plots mean/sd waveform.

filter_split_normalize_mean_vGRF.m

% relies on data and functions from https://github.com/alcantarar/dryft.
% Please cite the paper if you use it.
%
% This script reads in vgrf data, filters it, identifies stance phase,
% separates each stance phase, normalizes to stance phase, then plots 
% the mean +/- SD vgrf waveform.
%
% Ryan Alcantara | ryansalcantara@gmail.com


%% Read in data from force plate
% ripped straight from `sample.m`...
clear
close
GRF = dlmread('custom_drift_S001runT25.csv');  %  example from https://github.com/alcantarar/dryft

%% Apply Butterworth filter
Fs = 300; % From Fukuchi et al. (2017) dataset
Fc = 50;
Fn = (Fs/2);
n_pass = 2;
order = 2;
C = (2^(1/n_pass)-1)^(1/(2*order)); % Correction factor per Research Methods in Biomechanics (2e) pg 288
Wn = (tan(pi*Fc/Fs))/C; % Apply correction factor to adjusted cutoff freq
Fc_corrected = atan(Wn)*Fs/pi; % Hz
[b, a] = butter(order, Fc_corrected/Fn);

GRF_filt = filtfilt(b, a, GRF);

%% Identify where stance phase occurs (foot on ground)
[stance_begin,stance_end, good_stances] = dryft.split_steps(GRF_filt,... %vertical GRF
    140,... %threshold
    Fs,... %Sampling Frequency
    0.2,... %min_tc
    0.4,... %max_tc
    0); %(d)isplay plots = True

%% keep the stance phases that meet min/max tc requirements
b = stance_begin(good_stances);  % b == begin of stance phase
e = stance_end(good_stances);  % e == end of stance phase

%% Store all the stance phases sections of the waveform together
all_waveforms = cell(length(b),1);

for i = 1:length(b)
    all_waveforms{i} = GRF_filt(b(i):e(i));
end

% we can see that each stance phase has slightly different lengths
disp('number of frames for each identified stance phase:')
disp(cellfun(@length, all_waveforms))

%% Loop through each stance phase and normalize to 101 data points (0-100% stance phase)

% prep matrix where we'll store normalized vGRF waveforms
norm_mat = NaN(size(all_waveforms,1), 101);

for i = 1:size(all_waveforms,1)
    temp = all_waveforms{i};
    % fill any missing values because interpolation will fail otherwise
    temp_filled = fillmissing(temp, 'nearest',1);
    % interpolate to 101 points using interp1
    norm_mat(i,:) = interp1(linspace(0,100,length(temp_filled)), temp_filled, 0:100);
end

%% plot each step
hold on
for i = 1:size(norm_mat,1)
    plot(norm_mat(i,:), 'color', [0.7, 0.7, 0.7])
end

% calculate mean of waveforms
mean_vgrf = mean(norm_mat);
sd_vgrf = std(norm_mat);

% plot mean
plot(mean_vgrf, 'LineWidth', 2, 'color', 'black')
%plot mean +/- 1 SD
plot(mean_vgrf+sd_vgrf, 'LineWidth', 0.5, 'color', 'black')
plot(mean_vgrf-sd_vgrf, 'LineWidth', 0.5, 'color', 'black')