!pip freeze
Cython==0.19.1 -e git+git@github.com:moorepants/DynamicistToolKit.git@28505adf23b46bb5c3f3a0aa47342a6d8612dd65#egg=DynamicistToolKit-dev Jinja2==2.7.1 MarkupSafe==0.18 Pygments==1.6 Sphinx==1.2b1 argparse==1.2.1 coverage==3.6 docutils==0.11 ipython==1.0.0 matplotlib==1.3.0 nose==1.3.0 numpy==1.7.1 numpydoc==0.4 pandas==0.12.0 pyparsing==2.0.1 python-dateutil==2.1 pytz==2013b pyzmq==13.1.0 scipy==0.12.0 six==1.3.0 sympy==0.7.3 tornado==3.1 wsgiref==0.1.2
%pylab inline
Populating the interactive namespace from numpy and matplotlib
rcParams['figure.figsize'] = 16, 8
import pandas
from dtk import walk, process
Obinna provided me with some normal walking data. We exported the data from the D-Flow Motion Capture Module and included the HMB outputs. D-Flow sets missing values from HBM data to 0.000000
in the CSV file so I replace them with NA
when loading. I also set the resulting data frame's index to the TimeStamp
column.
obinna = pandas.read_csv('../data/obinna-walking.txt', delimiter='\t',
index_col="TimeStamp", na_values='0.000000')
The data frame has quite a few time series which are all floats except for the FrameNumber
column which is an integer.
obinna
<class 'pandas.core.frame.DataFrame'> Index: 5559 entries, 534.133513 to 589.712359 Columns: 617 entries, FrameNumber to HBM.COM.Z dtypes: float64(616), int64(1)
For this analysis we are only interested in the vertical ground reaction forces, the joint angles, the joint angular rates, and the joint torques. The D-Flow HBM computations give the angles and the torques. We will compute the rates ourselves.
[col for col in obinna.columns if col.endswith('.Ang') or col.endswith('ForY') or col.endswith('.Mom')]
['FP1.ForY', 'FP2.ForY', 'PelvisX.Ang', 'PelvisX.Mom', 'PelvisY.Ang', 'PelvisY.Mom', 'PelvisZ.Ang', 'PelvisZ.Mom', 'PelvisYaw.Ang', 'PelvisYaw.Mom', 'PelvisForwardPitch.Ang', 'PelvisForwardPitch.Mom', 'PelvisRightRoll.Ang', 'PelvisRightRoll.Mom', 'TrunkFlexion.Ang', 'TrunkFlexion.Mom', 'TrunkRightBend.Ang', 'TrunkRightBend.Mom', 'TrunkLeftTwist.Ang', 'TrunkLeftTwist.Mom', 'HeadFlexion.Ang', 'HeadFlexion.Mom', 'HeadRightBend.Ang', 'HeadRightBend.Mom', 'HeadLeftTwist.Ang', 'HeadLeftTwist.Mom', 'RShoulderUp.Ang', 'RShoulderUp.Mom', 'LShoulderUp.Ang', 'LShoulderUp.Mom', 'RShoulderForward.Ang', 'RShoulderForward.Mom', 'LShoulderForward.Ang', 'LShoulderForward.Mom', 'RShoulderInward.Ang', 'RShoulderInward.Mom', 'LShoulderInward.Ang', 'LShoulderInward.Mom', 'RShoulderFlexion.Ang', 'RShoulderFlexion.Mom', 'LShoulderFlexion.Ang', 'LShoulderFlexion.Mom', 'RShoulderAbduction.Ang', 'RShoulderAbduction.Mom', 'LShoulderAbduction.Ang', 'LShoulderAbduction.Mom', 'RShoulderInternalRotation.Ang', 'RShoulderInternalRotation.Mom', 'LShoulderInternalRotation.Ang', 'LShoulderInternalRotation.Mom', 'RElbowFlexion.Ang', 'RElbowFlexion.Mom', 'LElbowFlexion.Ang', 'LElbowFlexion.Mom', 'RForeArmPronation.Ang', 'RForeArmPronation.Mom', 'LForeArmPronation.Ang', 'LForeArmPronation.Mom', 'RWristFlexion.Ang', 'RWristFlexion.Mom', 'LWristFlexion.Ang', 'LWristFlexion.Mom', 'RHandAbduction.Ang', 'RHandAbduction.Mom', 'LHandAbduction.Ang', 'LHandAbduction.Mom', 'RHipFlexion.Ang', 'RHipFlexion.Mom', 'LHipFlexion.Ang', 'LHipFlexion.Mom', 'RHipAbduction.Ang', 'RHipAbduction.Mom', 'LHipAbduction.Ang', 'LHipAbduction.Mom', 'RHipInternalRotation.Ang', 'RHipInternalRotation.Mom', 'LHipInternalRotation.Ang', 'LHipInternalRotation.Mom', 'RKneeFlexion.Ang', 'RKneeFlexion.Mom', 'LKneeFlexion.Ang', 'LKneeFlexion.Mom', 'RAnklePlantarFlexion.Ang', 'RAnklePlantarFlexion.Mom', 'LAnklePlantarFlexion.Ang', 'LAnklePlantarFlexion.Mom', 'RFootPronation.Ang', 'RFootPronation.Mom', 'LFootPronation.Ang', 'LFootPronation.Mom', 'RToeFlexion.Ang', 'RToeFlexion.Mom', 'LToeFlexion.Ang', 'LToeFlexion.Mom']
The vertical ground reaction forces are stored as variables with the extension .ForY
and the right and left plates are designated by FP2
and FP1
respectively. The following plot shows the right foot vertical ground reaction force in Newtons.
obinna['FP2.ForY'].plot()
<matplotlib.axes.AxesSubplot at 0x5122bd0>
We will first use the ground reaction forces to find the heel strike and toe off times during gait. I'll select a portion of the data which has suitable GRF profiles for extracting the landmarks.
start = 500
stop = 3500
data = walk.WalkingData(obinna.iloc[start:stop].copy())
Now I compute all of the joint angular rates numerically with central differencing.
to_diff = [col for col in obinna.columns if col.endswith('.Ang')]
new_names = [name.split('.')[0] + '.Rate' for name in to_diff]
data.time_derivative(to_diff, new_names)
data.raw_data[new_names]
<class 'pandas.core.frame.DataFrame'> Index: 3000 entries, 539.133881 to 569.12273 Data columns (total 46 columns): PelvisX.Rate 2882 non-null values PelvisY.Rate 2882 non-null values PelvisZ.Rate 2882 non-null values PelvisYaw.Rate 2882 non-null values PelvisForwardPitch.Rate 2882 non-null values PelvisRightRoll.Rate 2882 non-null values TrunkFlexion.Rate 2882 non-null values TrunkRightBend.Rate 2882 non-null values TrunkLeftTwist.Rate 2882 non-null values HeadFlexion.Rate 2882 non-null values HeadRightBend.Rate 2882 non-null values HeadLeftTwist.Rate 2882 non-null values RShoulderUp.Rate 2882 non-null values LShoulderUp.Rate 2882 non-null values RShoulderForward.Rate 2882 non-null values LShoulderForward.Rate 2882 non-null values RShoulderInward.Rate 2882 non-null values LShoulderInward.Rate 2882 non-null values RShoulderFlexion.Rate 2882 non-null values LShoulderFlexion.Rate 2882 non-null values RShoulderAbduction.Rate 2882 non-null values LShoulderAbduction.Rate 2882 non-null values RShoulderInternalRotation.Rate 2882 non-null values LShoulderInternalRotation.Rate 2882 non-null values RElbowFlexion.Rate 2882 non-null values LElbowFlexion.Rate 2882 non-null values RForeArmPronation.Rate 2882 non-null values LForeArmPronation.Rate 2882 non-null values RWristFlexion.Rate 2882 non-null values LWristFlexion.Rate 2882 non-null values RHandAbduction.Rate 2882 non-null values LHandAbduction.Rate 2882 non-null values RHipFlexion.Rate 2882 non-null values LHipFlexion.Rate 2882 non-null values RHipAbduction.Rate 2882 non-null values LHipAbduction.Rate 2882 non-null values RHipInternalRotation.Rate 2882 non-null values LHipInternalRotation.Rate 2882 non-null values RKneeFlexion.Rate 2882 non-null values LKneeFlexion.Rate 2882 non-null values RAnklePlantarFlexion.Rate 2882 non-null values LAnklePlantarFlexion.Rate 2882 non-null values RFootPronation.Rate 2882 non-null values LFootPronation.Rate 2882 non-null values RToeFlexion.Rate 2882 non-null values LToeFlexion.Rate 2882 non-null values dtypes: float64(46)
The strikes and toe offs can then be computed for this data.
right_strikes, left_strikes, right_offs, left_offs = \
data.grf_landmarks('FP2.ForY', 'FP1.ForY', do_plot=True, threshold=28.0)
Now we can section the gait cycle based on the right foot and interpolate at 10 distinct equally spaced points in the gait cycle.
right_steps = data.split_at('right', num_samples=20)
data.plot_steps('FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom', marker='o')
array([<matplotlib.axes.AxesSubplot object at 0x6784b90>, <matplotlib.axes.AxesSubplot object at 0x60b0710>, <matplotlib.axes.AxesSubplot object at 0x5db1e50>, <matplotlib.axes.AxesSubplot object at 0x63273d0>], dtype=object)
We can see the mean gait cycle easily:
mean_right_steps = right_steps.mean(axis='items')
mean_right_steps[['FP2.ForY', 'RKneeFlexion.Ang', 'RKneeFlexion.Rate', 'RKneeFlexion.Mom']].plot(subplots=True)
array([<matplotlib.axes.AxesSubplot object at 0x65ff4d0>, <matplotlib.axes.AxesSubplot object at 0x512f9d0>, <matplotlib.axes.AxesSubplot object at 0x5c98510>, <matplotlib.axes.AxesSubplot object at 0x5cb0bd0>], dtype=object)
At this point we have the data in the form needed to identify a simple linear controller where we assume the sensors are fed back and subtract from a set point in the walking gait cycle, do give the error in the gait. This error vector is then multiplied by a gain matrix and added to the limit cycle control torques to power the plant in walking.
The following gives the equation that generates the controls given the error in the sensors at a single time step during a single foot step.
$$m_m(t) = m(t)_{lc} + K(t) s_e(t) = m(t)_{lc} + K(t) (s_{lc}(t) - s_m(t))$$Now rearrange the equations such that we have one linear in both the gains, $K(t)$, and in $m^*(t)$:
$$m_m(t) = m(t)_{lc} + K(t) s_{lc}(t) - K(t) s_m(t) = m^*(t) - K(t) s_m(t)$$Here I specify the hypothesized sensors and controls, and the data to initialized the solver:
sensors = ['RKneeFlexion.Ang',
'RKneeFlexion.Rate',
'RHipFlexion.Ang',
'RHipFlexion.Rate',
'LKneeFlexion.Ang',
'LKneeFlexion.Rate',
'LHipFlexion.Ang',
'LHipFlexion.Rate']
controls = ['RKneeFlexion.Mom',
'RHipFlexion.Mom',
'LKneeFlexion.Mom',
'LHipFlexion.Mom']
solver = walk.SimpleControlSolver(right_steps, sensors, controls)
gains, sensors = solver.solve()
gains.shape
(20, 4, 8)
solver.plot_gains(gains)
array([[<matplotlib.axes.AxesSubplot object at 0x6e93550>, <matplotlib.axes.AxesSubplot object at 0x6eb4bd0>, <matplotlib.axes.AxesSubplot object at 0x59d1f90>, <matplotlib.axes.AxesSubplot object at 0x59ebe10>, <matplotlib.axes.AxesSubplot object at 0x5b47bd0>, <matplotlib.axes.AxesSubplot object at 0x725de90>, <matplotlib.axes.AxesSubplot object at 0x727d2d0>, <matplotlib.axes.AxesSubplot object at 0x73e3c10>], [<matplotlib.axes.AxesSubplot object at 0x728dd90>, <matplotlib.axes.AxesSubplot object at 0x7559fd0>, <matplotlib.axes.AxesSubplot object at 0x75822d0>, <matplotlib.axes.AxesSubplot object at 0x76bccd0>, <matplotlib.axes.AxesSubplot object at 0x76e53d0>, <matplotlib.axes.AxesSubplot object at 0x7829710>, <matplotlib.axes.AxesSubplot object at 0x7996510>, <matplotlib.axes.AxesSubplot object at 0x79bc7d0>], [<matplotlib.axes.AxesSubplot object at 0x76ca210>, <matplotlib.axes.AxesSubplot object at 0x7b3ff50>, <matplotlib.axes.AxesSubplot object at 0x7ca6310>, <matplotlib.axes.AxesSubplot object at 0x7cc2d10>, <matplotlib.axes.AxesSubplot object at 0x7e06690>, <matplotlib.axes.AxesSubplot object at 0x7cd3850>, <matplotlib.axes.AxesSubplot object at 0x7f74990>, <matplotlib.axes.AxesSubplot object at 0x7f97c50>], [<matplotlib.axes.AxesSubplot object at 0x80ef750>, <matplotlib.axes.AxesSubplot object at 0x8258e10>, <matplotlib.axes.AxesSubplot object at 0x826e210>, <matplotlib.axes.AxesSubplot object at 0x83e3f50>, <matplotlib.axes.AxesSubplot object at 0x840d250>, <matplotlib.axes.AxesSubplot object at 0x72673d0>, <matplotlib.axes.AxesSubplot object at 0x857a050>, <matplotlib.axes.AxesSubplot object at 0x86d03d0>]], dtype=object)