What is wrong with my ankle joint moment estimation?

Hi all,

I collected human motion data (kinematics and ground reaction forces) during walking. I am trying to estimate ankle joint moments (3D) and later powers. I segmented (per stance phase on a force plate) 3D ankle joint angles, GRFs, ankle joint centes, foot and shank center of mass postions and accelerations, foot and shank angulat velocities and accelerations. I do not have center of pressure data, I approximate this. Below is my code that I use. However, I am not sure whether ankle moments are calculated correctly beacuse of the atypical spike in the sagital plane moment in the first 20 frames that you can see in tha image below. Can you help please?

% MATLAB Code for 3D Ankle Joint Torques Calculation with Segmented Data Structure

load Data_structure.mat 


% Constants 
g = 9.81; % m/s^2, gravitational acceleration

%Subject weight/height
mass = 84.6;
height = 1.55;

% segment mass (kg)
    massThigh = 0.142*mass; 
    massShank = 0.043*mass; 
    massFoot = 0.013*mass; 

% segment lenght (meters)
    L_thigh = 0.240*height; 
    L_shank = 0.262*height; 
    L_foot = 0.144*height; 

% Moment of inertia
I_foot  = diag(data.FootMoI);   % 3x3 matrix
I_shank = diag(data.ShankMoI);
I_thigh = diag(data.ThighMoI);

% Initialize moment arrays
nFrames = 101;
forcePlates = {'fp1', 'fp3', 'fp4'};
limbs = {'L', 'R'};
ankleMoments = struct();

% Loop through force plates
for fp = 1:length(forcePlates)
    fpName = forcePlates{fp};

    for limb = 1:length(limbs)
        limbName = limbs{limb};
        % Initialize moment data for each plate and limb
        ankleMoments.(fpName).(limbName) = zeros(3, nFrames);

        % Check if data exists for the current force plate and limb
        if isfield(data.Segmented_Angles, fpName) && isfield(data.Segmented_Angles.(fpName), limbName)
            if ~isempty(data.Segmented_Angles.(fpName).(limbName)) && isfield(data.Segmented_Angles.(fpName).(limbName){1, 1}, 'ankle')

                % Extract kinematic and GRF data
                ankleAngles = data.Segmented_Angles.(fpName).(limbName){1, 1}.ankle;
                ankleJC = [data.Segmented_Jointcenteres.(fpName).(limbName){1, 1}.ankle.ap;
                           data.Segmented_Jointcenteres.(fpName).(limbName){1, 1}.ankle.ml;
                           data.Segmented_Jointcenteres.(fpName).(limbName){1, 1}.ankle.ver];

                % COM positions and accelerations
                comFoot = data.Segmented_CoMpos.(fpName).(limbName){1, 1}.foot';
                foot_comAcc = data.Segmented_CoMacc.(fpName).(limbName){1, 1}.foot';

                % Angular velocities and accelerations
                foot_aVel = data.Segmented_aVel.(fpName).(limbName){1, 1}.foot';
                foot_aAcc = data.Segmented_aAcc.(fpName).(limbName){1, 1}.foot';

                % GRF data
                try
                    GRF_fp = [data.(['Kistler' num2str(fp) '_fx_stance']);
                              data.(['Kistler' num2str(fp) '_fy_stance']);
                              data.(['Kistler' num2str(fp) '_fz_stance'])];
                catch
                    GRF_fp = zeros(3, nFrames);
                end

                % Approximate COP position
                r_ankle_to_COP = zeros(3, nFrames);

                for i = 1:nFrames
                        % Get GRF vertical force to detect stance
                        grfZ = GRF_fp(3, :);
                        stanceFrames = find(grfZ > 10);
                    
                        % Skip if not in stance
                        if isempty(stanceFrames) || ~ismember(i, stanceFrames)
                            r_ankle_to_COP(:, i) = [0; 0; 0];
                            continue;
                        end
                    
                        % Use alpha based on stance progression
                        alpha = (i - stanceFrames(1)) / (stanceFrames(end) - stanceFrames(1));
                        alpha = min(max(alpha, 0), 1);  % clamp to [0, 1]
                    
                        % Assume foot runs anterior along Z (adjust if needed)
                        heel_to_toe_vector = [0; 0; 1] * L_foot;
                    
                        % CoP vector from ankle = moving from heel to toe
                        r_ankle_to_COP(:, i) = -0.5 * heel_to_toe_vector + alpha * heel_to_toe_vector;
                    end

                % Calculate Ankle Moments
                for i = 1:nFrames
                    % Inertial and gravitational forces
                    footWeight = [0; 0; -massFoot * g];
                    footInertialForce = -massFoot * foot_comAcc(:, i);

                    % Dynamic vector from ankle joint to CoM
                    r_ankle_to_CoM = comFoot(:, i) - ankleJC(:, i);

                    % Compute Ankle Moments
                    ankleMoments.(fpName).(limbName)(:, i) = cross(r_ankle_to_COP(:, i), GRF_fp(:, i)) + ...
                                                             cross(r_ankle_to_CoM, footWeight + footInertialForce) + ...
                                                             I_foot * foot_aAcc(:, i) + ...
                                                             cross(foot_aVel(:, i), I_foot * foot_aVel(:, i));
                end
            end
        end
    end
end

% Plotting Ankle Moments in the Sagittal Plane from one force plate

figure 
plot(ankleMoments.fp1.R(1,:))
xlabel('Frame');
ylabel('Ankle sagital Moment (Nm)');

[Matt J]

Your problem isn't a problem of Matlab knowledge. It's a problem of knowledge about ankle moments. You're likely the only one on here who knows anything about that, I'm afraid.

[Sam Chak]

@Tomaszzz, Your code is incomplete with missing values or parameters as it doesn't run when is clicked.

[Tomaszzz]

Thanks guys; @Sam Chak Sorry Sam, I have updated the code and here is the link to the input paramteres because the server does now allow me to attach it: Data structure

[Sam Chak]

I see, @Tomaszzz. It is currently unable to upload the .mat file to the server. I'm not an expert in kinesiology. However, Prof. @William Rose may be able to advise on one or two things. Do you apply Euler's equations to compute the ankle moments?

[William Rose]

@Sam Chak, yes, you do. But the moment M in your equation is the total moment acting on the (foot) segment. In this case we want to find the net moment acting at the ankle joint, which is considered the proximal end of the foot segment. The total moment acting on the foot segment arises from 3 sources: 1. ground reaction force, represented by force acting at a point usually near the toes (although the location varies during the stance phase, and is zero during the swing phase, when the foot is off the ground), 2. force pulling on the segment at the ankle joint, which is due to a combination of muscles and the weight of th body, and 3. torque at the ankle joint, due to muscles in the lower leg that pull on the foot. We measure (1) and we solve for (2) and (3), using the measured linear and rotational acceleraiton of the foot segment. Tomaszz is calculating and plotting (3).

In normal subjects, most linear and angular movement is in the sagittal plane, so we can get fairly accurate results with a 2D analysis (Eng JJ & Winter DA, 1995; McClay I & Manal K, 1999).

[Sam Chak]

Hi @William Rose

Thank you for your response. The information you provided in the Answer and subsequent comments is not only impressive but also helpful to the OP in understanding the methodologies used in modeling human ankle kinematics. 👍

[William Rose]

THank you @Sam Chak. Means a lot coming from you!

[William Rose]

@Sam Chak,

Recently you answered a question about design of a feedback block in a control system, root locus analysis, etc. It was a really good answer and it made me wish I remembered more from my one controls course. Your answer prompted me to go back and spend an hour or so with my 1980s copy of Design of Feedback Control Systems by Hostetter.

[William Rose]

@Sam Chak,

And you did again, tonight!

Answer 1

[William Rose]

@Tomaszzz,

I agree with you that the early peak in your calculated ankle moment is anomalous and probably indicates a problem in the data or the calculations. My initial reaction is that this is probably a consequence of the fact that you do not measure the center of pressure (COP) location. Errors in the estimated COP location could easily cause significant errors in ankle joint moment.

Is plantarflexion positive in your plot? The broad late peak makes me think it is.

Here are published examples of sagittal plane ankle moment.

In the plot above, plantar flexion is positive, moment is normalized by body weight, and a full stride is shown.

In the plot above, dorsiflexion is positive, moment is normalized by body weight, and the stance phase is shown. Focus on the normal traces.

In the plot above, dorsiflexion is positive, moment is normalized by body weight, and the stance phase is shown. The trace for walking is thicker and has a negative peak that occurs later.

In the plot above, dorsiflexion is positive, moment is not normalized by body weight, and the stance phase is shown.

Your plot says the units are N-m, but I wonder if you have normalized the moment by body weight. If N-m is correct (i.e. you have not normalized), then the broad peak in plantar flexion in the second half of stance, in your plot, is much smaller than what others have reported.

Hansen et al. (2004) J Biomech use a simpler method to estimate ankle moment - see their Fig 2. They say that RP Wells (1981) found that the simpler method gave very similar results to a full inverse dynamics approach.

I have not checked your inverse dynamics calculations.

[William Rose]

@Tomaszzz,

Please deifine the positive directions for the x,y,z axes, for the lab and for the foot coordinate system, and explain which vectors are expressed in which frames of reference. Your script includes

% Assume foot runs anterior along Z (adjust if needed)
heel_to_toe_vector = [0; 0; 1] * L_foot;

which suggests that +z points along the long axis of the foot, which corresponds to forward in the lab. But the script also has

% Inertial and gravitational forces
footWeight = [0; 0; -massFoot * g];

which suggests that +z points upward, since the weight vector is negative and along z.

disp(data.FootMoI')
0.0037    0.0034    0.0009

This suggests that the +Z axis is the long axis of the foot, because the minimum moment of inertia for the foot will be the moment about the foot's long axis.

You plot the sagittal plane ankle moment with

plot(ankleMoments.fp1.R(1,:))

which indicates that the x axis is normal to the sagittal plane. Does +x point to the subject's right, and +y forward, and +z upward?

What are the units for FootMoI? Is the FootMoI about the c.o.m. of the foot or about the proximal joint (ankle joint)?

You have foot wt=0.013*M=1.01 kg (M=body mass=84.6 kg), and foot length=0.144*H=0.223 m (H=height=1.55 m). (These values are similar to Winter 4th ed., Table 4.1 and Fig. 4.1, which have foot wt=0.0145*M and foot length=0.152*H.) Winter Table 4.1 says foot radius of gyration (in sagittal plane, about the c.o.m.)=0.475*footlength=0.106 m. Then the MoI=mass(foot)* = 0.0113 kg-m^2. This is 3 times larger than data.FootMoI(1). How do you explain the factor-of-three difference between the MoI derived from Winter's anthropometric table and the value of data.FootMoI?

Please include comments thorughout your code indicating the units. When data are read in from an external file, as in this case, include explanatory comments about th units of the data in the file. The more the merrier. It is hard to underrstand your own code or find errors if the inits are not clear - normalized or not, etc.

[William Rose]

@Tomaszzz,

You said that you approximate the COP position relaitve to the ankle joint. The script has

% Use alpha based on stance progression
alpha = (i - stanceFrames(1)) / (stanceFrames(end) - stanceFrames(1));
alpha = min(max(alpha, 0), 1);  % clamp to [0, 1]

% Assume foot runs anterior along Z (adjust if needed)
heel_to_toe_vector = [0; 0; 1] * L_foot;

% CoP vector from ankle = moving from heel to toe
r_ankle_to_COP(:, i) = -0.5 * heel_to_toe_vector + alpha * heel_to_toe_vector;

In other words, alpha progresses linearly from 0 at beginning of stance to 1 at the end of stance, and the COP position progreeses linearly from half-a-foot length behind the ankle at start of stance to half-a-foot-length in front of the ankle at the end of stance.

That seems wrong to me, and I wonder if it acocunts for the abnormal shape of the ankle moment versus time plot. I suspect the COP is less far behind the ankle at the start and may be farther ahead of the ankle at the end of stance. How did you validate or check your assumption about the COP position relative to the ankle joint during stance?

.

[William Rose]

@Tomaszzz,

You can compare the COP position relative to ankle, which you estimate with your alpha-model, to the ankle to the data published in Winter, Biomechanics and Motor Control of Human Movement, 4th ed., 2009. Appendix A of the book has the data, and is available as a PDF here. Table A.2(c), pp.311-315, includes right ankle X and Y (sagittal plane) coordinates . Table A.5(a), pp.346-349, has ground reaction forces RX, RY, and the x-coordinate of the COP. (The Y coordinate is zero, of course.) The data in the tables shows that at the time of right heel contact (frame 28, "HCR"), ankleX=1.254 and COPX=1.227, i.e. the COP is 0.027 m behind the ankle joint, i.e. about 0.1*footlength behind the ankle. Your alpha-model for COP has the COP 0.5*foot length behind the ankle, at the start of stance. At the end of stance (frame 69 in the appendix data, last frame before "TOR"=toe-off right), ankle X=1.462 and COPX=1.486, so the COP is 0.024 m in front of the ankle, about 0.1 foot length. Of course the ankle has risen by this point, so the foot is not flat on the ground. You can compare these data to your alpha-model of how the COP moves during stance. I am becoming more confident that the unusual shape of your calculated ankle monet versus time plot is due to the way the COP is estim,ated. Assuming that the claculations are otherwise correct, which I have not checked.

[Tomaszzz]

@William Rose; Dear William, many thanks for spending time familirazing yourself with the problem and writing back. Your answers have helped me better understand what I am doing. I am yet to cross check all the matters you have raised but I have just quickly modified the CoP estimations based on your input and indeed it seems the problem was here. Many thanks!

[William Rose]

@Tomaszzz,

That shape looks a lot better! And now that you have changed the units to N-m/kg, the magnitude also looks reasonable. Good luck with your research.