-
Notifications
You must be signed in to change notification settings - Fork 0
/
mos.m
120 lines (108 loc) · 5.39 KB
/
mos.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
function result = mos(data, tmdata, detail)
% analysis of Margin of Stability
% data should be a data structure from a mocap file, from getdata.m
% tmdata should be the corresponding treadmill data, also from getdata.m
% use detail=1 to show and pause the results
%
% result is a table with one row and the following columns:
% MOS_AP_left, MOS_AP_right, MOS_ML_right, MOS_ML_left
% These are measured at each heelstrike, according to McAndrew Young et al 2012.
% if no input is specified, we use one particular file for testing
if nargin < 1
detail = 1;
[data,tmdata] = getdata('Par7_POST\Mocap0001.txt', 1); % we only need the mocap data
end
fprintf('MoS analysis for %s\n', data.name);
% find the heelstrikes and only keep those that occurred between 10 and 30 seconds
% (NOTE: we could add additional frames after the perturbation, if we know when the perturbation happened)
Lhs = data.Lhs;
Rhs = data.Rhs;
st = 10; % starting time
ed = 30; % ending time
time = data.data(:,1)-data.data(1,1); % time relative to start of file
Lhs = Lhs( (time(Lhs) >= st) & (time(Lhs) <= ed) );
Rhs = Rhs( (time(Rhs) >= st) & (time(Rhs) <= ed) );
% extract the required marker data
SACRx = getcolumn(data, 'SACR.PosX');
SACRy = getcolumn(data, 'SACR.PosY');
SACRz = getcolumn(data, 'SACR.PosZ');
LHEEy = getcolumn(data, 'LHEE.PosY');
LHEEz = getcolumn(data, 'LHEE.PosZ');
RHEEy = getcolumn(data, 'RHEE.PosY');
RHEEz = getcolumn(data, 'RHEE.PosZ');
LTOEz = getcolumn(data, 'LTOE.PosZ');
RTOEz = getcolumn(data, 'RTOE.PosZ');
LMT5x = getcolumn(data, 'LMT5.PosX');
RMT5x = getcolumn(data, 'RMT5.PosX');
% calculate the mean distance from heel marker to Sacrum (in the sagittal YZ plane)
frames = (time >= st) & (time <= ed); % use only these frames for length calculation
Rdistance = sqrt( (SACRz-RHEEz).^2 + (SACRy-RHEEy).^2 );
Ldistance = sqrt( (SACRz-LHEEz).^2 + (SACRy-LHEEy).^2 );
L = mean([Rdistance(frames) ; Ldistance(frames)], 'omitnan'); % average leg length, used for XCoM calculation
% find the treadmill belt speed between 10 and 30 seconds
tmtime = tmdata.data(:,1)-tmdata.data(1,1); % time relative to start of file
frames = (tmtime >= st) & (tmtime <= ed); % use only these frames for length calculation
beltspeed = mean(tmdata.data(frames,3)); % from column 3
% calculate the Sacrum velocity in X and Z direction, and the extrapolated center of mass
[SACRvx, SACRx] = velocity(time, SACRx);
[SACRvz, SACRz] = velocity(time, SACRz);
% the belt moves in the +Z direction, so we need to subtract belt speed to get SACR velocity relative to belt
SACRvz = SACRvz - beltspeed;
g = 9.81;
w0 = sqrt(g/L);
XCoMx = SACRx + SACRvx / w0;
XCoMz = SACRz + SACRvz / w0;
% find the MOS at the left and right heelstrikes
% sign convention: MOS is positive when XCoM is inside BOS
MOS_AP_left = ( XCoMz(Lhs) - LTOEz(Lhs) ); % because Z is posterior
MOS_AP_right = ( XCoMz(Rhs) - RTOEz(Rhs) ); % because Z is posterior
MOS_ML_left = ( XCoMx(Lhs) - LMT5x(Lhs) ); % because X is to the right
MOS_ML_right = -( XCoMx(Rhs) - RMT5x(Rhs) ); % because X is to the right
% calculate mean and SD for all variables, and store in result
MOS_AP_left_mean = mean(MOS_AP_left);
MOS_AP_left_SD = std(MOS_AP_left);
MOS_AP_right_mean = mean(MOS_AP_right);
MOS_AP_right_SD = std(MOS_AP_right);
MOS_ML_left_mean = mean(MOS_ML_left);
MOS_ML_left_SD = std(MOS_ML_left);
MOS_ML_right_mean = mean(MOS_ML_right);
MOS_ML_right_SD = std(MOS_ML_right);
result = table(MOS_AP_left_mean, MOS_AP_left_SD, MOS_AP_right_mean, MOS_AP_right_SD, ...
MOS_ML_left_mean, MOS_ML_left_SD, MOS_ML_right_mean, MOS_ML_right_SD);
if (detail)
% plot the mos variables for all heelstrikes
close all
figure(1)
subplot(2,1,1)
plot(MOS_AP_left);hold on; plot(MOS_AP_right);
ylabel('AP MOS (m)')
subplot(2,1,2)
plot(MOS_ML_left);hold on; plot(MOS_ML_right);
xlabel('step number')
ylabel('ML MOS (m)')
legend('Right','Left');
fprintf('Length L for XCoM calculation: %8.3f m.\n', L);
fprintf('Left AP MoS: %8.3f %s %8.3f s.\n',MOS_AP_left_mean, char(177),MOS_AP_left_SD);
fprintf('Right AP MoS: %8.3f %s %8.3f s.\n',MOS_AP_right_mean,char(177),MOS_AP_right_SD);
fprintf('Left ML MoS: %8.3f %s %8.3f s.\n',MOS_ML_left_mean, char(177),MOS_ML_left_SD);
fprintf('Right ML MoS: %8.3f %s %8.3f s.\n',MOS_ML_right_mean,char(177),MOS_ML_right_SD);
disp('Check Figure 1, and the results printed above, for problems');
disp('Hit ENTER to Continue')
pause
end
disp('MoS analysis (mos.m) completed.');
end
%==================================================
function [v,x] = velocity(t,x)
% find the velocity v(t) for the position data x(t)
% outputs the velocity v and the low-pass filtered position x
cutoff_freq = 6.0; % Hz, low pass filter to reduce noise in velocity
Fs = 1/mean(diff(t)); % sampling frequency
if std(diff(t)) > 0.001
error('velocity calculation: sampling rate is not constant');
end
[b, a] = butter(2, (cutoff_freq/(Fs/2)));
x = filtfilt(b, a, x);
v = diff(x)./diff(t);
v = [v(1) ; v]; % duplicate sample 1 so we have same number of samples in v and x
end