Skip to content

Commit

Permalink
Final Version
Browse files Browse the repository at this point in the history
  • Loading branch information
@lepremiere
lepremiere committed Mar 1, 2021
1 parent ae132e6 commit d657d85
Show file tree
Hide file tree
Showing 34 changed files with 2,570 additions and 1,577 deletions.
58 changes: 37 additions & 21 deletions Check_CMJ.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@
% CMJdata = Check_CMJ()
%
% Selected files must be ".csv" of force plate data, with vertical ground
% reaction forces stored in variables "z1-z4" and time in last column.
% reaction forces stored in variables containing "z" and time in last column.
%
% OUTPUT
% CMJdata: Cell array (3xn) with filenames, -paths, and data for each selected file will be
% returned. Press button next to accept file or delete to drop it.
% returned. Press button next to accept file or discard to drop it.
%
%---------------------------------------------------------------------------------------------------
% Latest Edit: 31.August.2020
Expand All @@ -28,20 +28,21 @@
%% Looping through files
for ii = 1:length(files(:,1))
%% Determining total vertical ground reaction force (VGRF)
% Note: Force plate data must be table with VGRF stored in columns named 'z1-z4' and time stored in last column. Time must be consistent
% Note: Force plate data must be table with VGRF stored in columns with names containing "z" and time stored in last column. Time must be consistent
data = files{ii,3};
Fz = data.z1 + data.z2 + data.z3 + data.z4;
inds = contains(data.Properties.VariableNames, 'z'); % Finding columns with z-forces
Fz = sum(data{:,inds}, 2);

% Determines the the sampling frequency 'Fs'
Fs = 1/(data.(data.Properties.VariableNames{end})(11) - data.(data.Properties.VariableNames{end})(10));

%% Filtering the data
Fc = 50; % Cutt-off frequency
[b,a] = butter(3,(Fc/(Fs/2)/0.8022),'low'); % Calculating coefficients of second order lowpass butterworth filter. Adjusted by 0.8022
[b,a] = butter(3,(Fc/(Fs/2)/0.8022),'low'); % Calculating coefficients of third order lowpass butterworth filter. Adjusted by 0.8022
Fz_filtered = filtfilt(b,a,Fz); % Filtering VGRF with zero phase lowpass butterworth filter

%% Adjusting VGRF in case force plate was relocated
if (mean(Fz_filtered) < 50)
if (mean(Fz_filtered) < 250)
minimum = sortrows(Fz_filtered).*-1;
relocation = mean(minimum(1:500));
Fz_filtered = Fz_filtered + relocation;
Expand All @@ -65,19 +66,18 @@
% Finding window with smallest variability (std), assuming to be the best baseline and calculating bodyweight from its mean VGRF
% Also, sets variability as base line noise
[row ~] = find(q(:,1) == min(q(:,1)));
BW = mean(Fz_filtered(1:1500)./9.81); % m = F/a
noise = Fz - Fz_filtered;
threshold = range(noise(row:row + interval)); % Sets variability as base line noise
BW = mean(Fz_filtered(1:2000)./9.81); % m = F/a
threshold = range(Fz(row:row + interval)); % Sets variability as base line noise

%% Determining start, takeoff, landing and airtime of CMJ
% Getting orientation in the jump
subzero = find(Fz_filtered<10); % Getting indices of no contact to force plate
takeoff = find(Fz_filtered(row:subzero(1)) < threshold,1,'first')+row; % Determines takeoff as first VGRF smaller than baseline noise. +row accounts for subzero offset
landing = find(Fz_filtered(takeoff:end) > threshold,1,'first'); % Determines landing as first VGRF greater than baseline noise after takeoff
airtime = landing*1/Fs; % Landing ontains number of sampling points for the jump
posStart = find(Fz_filtered(row:subzero(1)-1) > Fz_filtered(row) + threshold,1,'first') + row - 1000; % Determines start of movement by finding the first VGRF greater than baseline VGRF + noise. 1000 ms offset
negStart = find(Fz_filtered(row:subzero(1)-1) < Fz_filtered(row) - threshold,1,'first') + row - 1000; % Determines start of movement by finding the first VGRF greater than baseline VGRF - noise. 1000 ms offset
start = min([posStart,negStart]); % Finds true start
subzero = find(Fz_filtered < 0); % Getting indices of no contact to force plate
takeoff = find(Fz_filtered(row:subzero(1)) < threshold,1,'first')+row; % Determines takeoff as first VGRF smaller than baseline noise. +row accounts for subzero offset
landing = find(Fz_filtered(takeoff:end) > threshold,1,'first'); % Determines landing as first VGRF greater than baseline noise after takeoff
airtime = landing * 1/Fs; % Landing obtains number of sampling points for the jump. By multiplying with sampling frequency, airtime is calculated
posStart = find(Fz_filtered(row:subzero(1) - 1) > Fz_filtered(row) + threshold*2,1,'first') + row - 500; % Determines start of movement by finding the first VGRF greater than baseline VGRF + noise. 500 ms offset
negStart = find(Fz_filtered(row:subzero(1) - 1) < Fz_filtered(row) - threshold*2,1,'first') + row - 500; % Determines start of movement by finding the first VGRF greater than baseline VGRF - noise. 500 ms offset
start = min([posStart,negStart]); % Finds true start

%% Calculatng net force and net impulse
netForce = Fz_filtered - BW*9.81 ; % F(net) = VGRF - BW*g
Expand All @@ -100,7 +100,11 @@

% Calculates the RFDs during jump with respect to time
for i = 1:length(jump)-rfd_interval
rFD(i) = ((jump(i+rfd_interval)-jump(i))/rfd_interval)*Fs;
if(jump(i) > BW*9.81)
rFD(i) = ((jump(i+rfd_interval)-jump(i))/rfd_interval)*Fs;
else
rFD(i) = 0;
end
end

% Finds max RFD and time to max RFD
Expand All @@ -121,15 +125,27 @@

hold on
h1 = area([start:takeoff],netForce(start:takeoff)); % Plots impulse as filled area
h1.LineWidth = 2;

% Setting labels
titl = sprintf('File: %s \nJumpheight: %.2f / %.2f cm, Peak Power: %.2f W', string(files(ii,2)), jumpHeightImpulse, jumpHeightAirtime, peakPower);
title(titl, 'Interpreter', 'none');
ylabel(ax(1), 'Force (N)');
ylabel(ax(2), 'Velocity (m/s)');
xlabel('Time (ms)');
ylabel(ax(1), 'Force [N]');
ylabel(ax(2), 'Velocity [m/s]');
xlabel('Time [ms]');
legend([h(1), h(2), h1], {'VGRF', 'Velocity', 'Impulse'});

% Appereance behavior
h(1).LineWidth = 2;
ax(1).YTick = [round(min(ax(1).YLim),0):range([round(min(ax(1).YLim),0),round(max(ax(1).YLim),0)])/10:round(max(ax(1).YLim),0)];
ax(1).LineWidth = 2;
ax(1).FontWeight = 'bold';
ax(1).Box = 'off';
h(2).LineWidth = 2;
ax(2).YTick = [round(min(ax(2).YLim),0):range([round(min(ax(2).YLim),0),round(max(ax(2).YLim),0)])/10:round(max(ax(2).YLim),0)];
ax(2).LineWidth = 2;
ax(2).FontWeight = 'bold';

% Line markers for different events
xline(row,'','Baseline','LabelHorizontalAlignment','center','FontWeight','bold'); % Baseline start
xline(row + interval,'','LabelHorizontalAlignment','center','FontWeight','bold'); % Baseline end
Expand All @@ -154,5 +170,5 @@
files(f,:) = [];
clear('global')
CMJdata = files;
end
end

69 changes: 46 additions & 23 deletions Check_ITT.m → Check_DYNO.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function ITTdata = Check_ITT(varargin)
function DYNOdata = Check_DYNO(varargin)
% This function calculates heart rate variability measures for a nx1 array of RR-intervals
% stored in a ".txt" file. Results will be visualized and files can be selected for
% further processing.
Expand All @@ -8,7 +8,7 @@
% ITTdata = Check_ITT('key', 'h')
%
% INPUT
% key: Character that defines which button press indicates a stimulation.
% key: Character that defines which button press indicates a stimulation. Default: "g".
%
% OUTPUT
% ITTdata: Cell array (3xn) with filenames, -paths, and data for each selected file will be
Expand Down Expand Up @@ -37,16 +37,17 @@
%% Looping through files
for ii = 1:length(files(:,1))
f = [];
data = files{ii,3}; % Getting raw data
Fs = 1/data.Torque.interval; % Getting sampling frequency
Fc = 20; % Cut-off frequency
[B,A] = butter(3,(Fc/(Fs/2)/0.8022),'low'); % Calculating coefficients of 3rd order lowpass butterworth filter
torque = filtfilt(B,A,data.Torque.values); % Filtering torque with zero phase butterworth filter
data = files{ii,3}; % Getting raw data
Fs = 1/data.Torque.interval; % Getting sampling frequency
Fc = 20; % Cut-off frequency
[B,A] = butter(3,(Fc/(Fs/2)/0.8022),'low'); % Calculating coefficients of 3rd order lowpass butterworth filter
torque = filtfilt(B,A,data.Torque.values); % Filtering torque with zero phase butterworth filter
key_presses = data.Keyboard.codes(:,1) == double(ip.key);
stim_time = round((fix(data.Keyboard.times(key_presses)*Fs)),0); % Getting the indices of stimulations indicated by keyboardpresses
time = [1:length(torque)]./Fs'; % Getting time in ms or n/Fs
max_stim_w = 250; % Setting the window following stimulation to analyse peak to peak torque in
min_stim_w = 50; % Setting the window following stimulation to analyse peak to peak torque in
stim_time = round((fix(data.Keyboard.times(key_presses)*Fs)),0); % Getting the indices of stimulations indicated by keyboardpresses
time = [1:length(torque)]./Fs'; % Getting time in ms or n/Fs
max_stim_w = 250; % Setting the window following stimulation to analyse peak to peak torque in
min_stim_w = 50; % Setting the window following stimulation to analyse peak to peak torque in
rfdWindow = 10;

% Error catch function if keypresses do not contained desired key
if(length(stim_time) == 0)
Expand All @@ -64,33 +65,55 @@
max(torque(index:index + max_stim_w))]);
twitchTorqueRest = range([min(torque(index_2:index_2 + min_stim_w)),... % Peak to peak torque within resting stimulation + window
max(torque(index_2:index_2 + max_stim_w))]);
MVC = range([max(torque(index:index + max_stim_w)),... % Peak to peak torque within stimulation - window*20 (~5s)
MVC = range([max(torque(index - 5000:index)),... % Peak to peak torque within stimulation - window*20 (~5s)
min(torque(index_2 - min_stim_w:index_2))]);
VA = (1 - (twitchTorqueMVC*(torque(index)/MVC)/twitchTorqueRest))*100; % VA = (1 - (superimposed twitch*(Tb/MVC)*resting twitch^-1)) *100.

rest_contraction = torque(index_2:index_2 + max_stim_w*2);

% Rate of force development of the resting twitch
for j = 1:length(rest_contraction) - rfdWindow
rfd(j) = (rest_contraction(j+rfdWindow) - rest_contraction(j))/rfdWindow*Fs;
end
[peakRFD_rest, ind_rfd] = max(rfd);
ind_rfd = ind_rfd + index_2;

%% Plot
clf
hold on
index = stim_time(2*i-1) - screeningWindow;
index_2 = stim_time(2*i-1) + screeningWindow;
h = plot(time(index:min([index_2, length(time)])),... % Plots data of the specific stimulation
torque(index:index_2));

index = stim_time(2*i-1) - screeningWindow;
index_2 = stim_time(2*i) + screeningWindow/2;
h = plot(time(index:min([index_2, length(time)])),... % Plots data of the specific stimulation
torque(index:index_2), 'LineWidth', 2);

% Plotting markers that indicate the data points that are used for calculation
for j = [0, 1]
[~, ind] = min(torque(stim_time(2*i-j):stim_time(2*i-j) + min_stim_w));
plot(time(ind + stim_time(2*i-j)), torque(stim_time(2*i-j) + ind), 'xg', 'MarkerSize', 8, 'LineWidth', 2);
h(1) = plot(time(ind + stim_time(2*i-j)), torque(stim_time(2*i-j) + ind), 'xg', 'MarkerSize', 8, 'LineWidth', 2);
[~, ind] = max(torque(stim_time(2*i-j):stim_time(2*i-j) + max_stim_w));
plot(time(ind + stim_time(2*i-j)), torque(stim_time(2*i-j) + ind), 'xr', 'MarkerSize', 8, 'LineWidth', 2);
h(2) = plot(time(ind + stim_time(2*i-j)), torque(stim_time(2*i-j) + ind), 'xr', 'MarkerSize', 8, 'LineWidth', 2);
h(3) = plot(time(ind_rfd), torque(ind_rfd), 'xc', 'MarkerSize', 8, 'LineWidth', 2);
if(j == 1)
[~, ind] = max(torque(stim_time(2*i-j)-2000:stim_time(2*i-j)));
h(4) = plot(time(stim_time(2*i-j)-2000+ind), torque(stim_time(2*i-j)-2000+ind), 'xb', 'MarkerSize', 8, 'LineWidth', 2);
end
end

% Plot behavior
xlim([index, index_2]/Fs);
titl = sprintf('File: %s \n MVC: %.0f Nm, Superimposed: %.0f Nm, Resting: %0.f Nm, VA: %.0f %% ', [char(files{ii,2}), '_', num2str(i)], MVC, twitchTorqueMVC, twitchTorqueRest, VA );
ylim([min(torque(index:index_2))-50, max(torque(index:index_2))*1.1]);
titl = sprintf('File: %s \n MVT: %.0f Nm, Superimposed: %.0f Nm, Resting: %0.f Nm, Peak RTD: %.0f Nm/s VA: %.0f %% ',...
[char(files{ii,2}), '_', num2str(i)], MVC, twitchTorqueMVC, twitchTorqueRest, peakRFD_rest, VA );
title(titl, 'Interpreter', 'none');
xlabel('Time [s]');
ylabel('Torque [Nm]');
xline(time(stim_time(2*i-1)),'','Superimposed','LabelHorizontalAlignment','center','LabelVerticalAlignment','mid','FontWeight','bold'); % Sets vertical line to first (superimposed) stimulation
xline(time(stim_time(2*i)),'','Resting','LabelHorizontalAlignment','center','LabelVerticalAlignment','mid','FontWeight','bold'); % Sets vertical line to second (resting) stimulation
legend(h(1:4), {'Min', 'Max', 'PeakRFD', 'PeakTorque'})
ax = gca;
ax.FontWeight = 'bold';
ax.LineWidth = 2;
hold off

%% Control loop
c = uicontrol('Position',[10 80 75 50],'String', 'Accept', 'FontSize', 12, 'FontWeight', 'bold', 'BackgroundColor', 'g', 'Callback','uiresume(gcf);'); % Callback for button accept. uiresume stops script temporarily
d = uicontrol('Position',[10 5 75 50], 'String', 'Discard', 'FontSize', 12, 'FontWeight', 'bold', 'BackgroundColor', 'r', 'Callback','global f i; f = [f,i];uiresume()'); % Callback for button discard. Saves discarded file for later removement
Expand All @@ -111,5 +134,5 @@
%% Create outputs
files(F,:) = []; % Deletes all input files that do not conatain any selected stimulation
clear('global')
ITTdata = files; % Returns cleaned files
DYNOdata = files; % Returns cleaned files
end
Loading

0 comments on commit d657d85

Please sign in to comment.