Calculate joint stiffness using Hamill et al. 2014 method

calculate_joint_stiffness.m

function k_joint = calculate_joint_stiffness(ankle_angle, joint_angle, joint_moment, joint_name)

%Calculates joint stiffness per the Hamill et al 2014 method
% John J Davis IV
%IU Biomechanics Lab
%   @JDruns
%30 June 2020

%                     ---  INPUT: --- 
% ankle_angle - (vector) ankle angle data in degrees, respecting right hand rule

% joint_angle - (vector) angle data for the joint in question. Might be a duplicate of ankle_angle
%               if you are calculating ankle stiffness. In degrees; respects right hand rule

% joint_moment - (vector) joint moment data for the joint in question. Respects right hand rule.

% joint_name - (string) one of 'ankle', 'knee', 'hip'

%                     ---  OUTPUT:  ---

% k_joint - (scalar) stiffness for the joint in question. Units will be:
%                         (whatever units moment is in) / degree

% Please note that I (JJD) do not believe that hip stiffness has a meaningful interpretation
% for data collected in runners.


%Flip signs as needed so stiffness at ankle and knee follow literature traditions of being positive
switch lower(joint_name)
    case 'ankle'
        %something
        joint_moment = joint_moment * -1;
        %Flip moment to positive
    case 'knee'
        %Something else
        joint_angle = joint_angle * -1;
        %Flip so knee extension is positive angle
    case 'hip'
        %flip so positive = hip EXTENSION moment
        joint_moment = joint_moment * -1;
    otherwise
        error('Joint not supported! Check name');
end

% --- Get index of max ankle dorsiflexion
[~, index_of_max_ankle_flexion] = max(ankle_angle);

%Catch cases where max dorsiflex is too early (15% of stance)
if index_of_max_ankle_flexion < 0.15*length(ankle_angle)
    n_rest_of_stance = round(0.15*length(ankle_angle));
    [~, index_of_max_ankle_flexion] = max(ankle_angle(n_rest_of_stance:end));
    index_of_max_ankle_flexion = index_of_max_ankle_flexion+n_rest_of_stance-1;
    warning('High initial contact dorsiflexion detected, searching for max in last 85% of stance...');
end

joint_angle_during_absorption = joint_angle(1:index_of_max_ankle_flexion);
joint_moment_during_absorption = joint_moment(1:index_of_max_ankle_flexion);


%Fit linear regression to moment as a function of angle
joint_lm = fitlm(joint_angle_during_absorption, joint_moment_during_absorption);

%The slope is the estimate for stiffness
k_joint = joint_lm.Coefficients.Estimate(2);

end