Marcos Duarte
Laboratory of Biomechanics and Motor Control
Federal University of ABC, Brazil
import sys, os
sys.path.insert(1, r'./../functions')
import numpy as np
import pandas as pd
pd.set_option('precision', 4)
# matplotlib configuration:
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib inline
import seaborn as sns
sns.set_context("notebook", font_scale=1.4,
rc={'font.size': 16, 'lines.linewidth': 2,\
'lines.markersize': 10, 'axes.titlesize': 'x-large'})
matplotlib.rc('legend', numpoints=1, fontsize=16)
# IPython:
from IPython.display import display, Image, IFrame
import ipywidgets
from ipywidgets import FloatProgress, interactive
The control system for equilibrium of the human body includes two subsystems:
Gurfinkel et al. (1995) Kinesthetic reference for human orthograde posture. Neuroscience, 68, 229.
The supporting surface where the subjects stood was tilted slowly (0.04 $^o/s$).
During the tilting, small high frequency oscillations of the body were superimposed on a large slow body movements.
The usual process of stabilization of the body continued, but the instant equilibrium was maintained relative to a slowly changing position, rather than around a fixed set point.
See Zatsiorsky & Duarte (2000), for a review about findings supporting a two-subsystems control.
Image(filename='../images/forcevectors.png', width=500)
display(Image(filename='../images/prolongedstanding.png', width=600))
display(Image(filename='../images/fractal.png', width=500))
print('Duarte & Zatsiorsky (2000)')
Duarte & Zatsiorsky (2000)
The reference point migration is called rambling and the COP migration around the reference was coined trembling.
IFrame('http://pesquisa.ufabc.edu.br/bmclab/pubs/mc99b.pdf', width='100%', height='500px')
IFrame('http://pesquisa.ufabc.edu.br/bmclab/pubs/mc00.pdf', width='100%', height='500px')
Image(filename='../images/rambtremb.png', width=600)
The rambling and trembling components of the COP trajectory were computed in the following way:
def detect_zerocross(y, method='no_zero'):
"""
Zero-crossing detection, index before a zero cross.
"""
import numpy as np
y = np.asarray(y)
if method == 'no_zero':
# doesn't detect if a value is exactly zero:
inds0 = np.where(y[:-1] * y[1:] < 0)[0]
else:
# detects if a value is exactly zero:
inds0 = np.where(np.diff(np.signbit(y)))[0]
return inds0
def iep_decomp(time, force, cop):
"""
Center of pressure decomposition based on the IEP hypothesis.
"""
import numpy as np
from scipy.interpolate import UnivariateSpline
force = force - np.mean(force)
# detect zeros
inds0 = detect_zerocross(force)
# select data between first and last zeros
time = time[inds0[0]:inds0[-1]+1]
force = force[inds0[0]:inds0[-1]+1]
cop = cop[inds0[0]:inds0[-1]+1]
# IEP0 & IEP:
iep0 = inds0 - inds0[0]
spl = UnivariateSpline(time[iep0], cop[iep0], k=3, s=0)
rambling = spl(time)
trembling = cop - rambling
return rambling, trembling, iep0, inds0
# data available at https://doi.org/10.6084/m9.figshare.3394432.v2
# load some data
# locally:
path2 = r'./../../../X/BDB/'
filename = os.path.join(path2, 'BDS00038.txt')
grf = pd.read_csv(filename, delimiter='\t', skiprows=1, header=None,
names=['Time','Fx','Fy','Fz','Mx','My','Mz','COPx','COPy'],
engine='c')
t, Fx, Fy, Fz, Mx, My, Mz, COPx, COPy = [_ for _ in grf.values.T]
def plot_grf(grf):
t, Fx, Fy, Fz, Mx, My, Mz, COPx, COPy = [_ for _ in grf.values.T]
Funits, Munits, COPunits = 'N', 'Nm', 'cm'
plt.figure(figsize=(12, 7))
gs1 = gridspec.GridSpec(3, 1)
gs1.update(bottom=0.52, top=0.96, hspace=0.12, wspace=.2)
ax1, ax2, ax3 = plt.subplot(gs1[0]), plt.subplot(gs1[1]), plt.subplot(gs1[2])
gs2 = gridspec.GridSpec(3, 3)
gs2.update(bottom=0.08, top=0.42, wspace=0.3)
ax4, ax5 = plt.subplot(gs2[:, :-1]), plt.subplot(gs2[:, 2])
ax1.set_ylabel('Fx (%s)' %Funits)
ax1.set_xticklabels([]), ax1.locator_params(axis='y', nbins=4)
ax1.yaxis.set_label_coords(-.07, 0.5)
ax2.set_ylabel('Fy (%s)' %Funits)
ax2.set_xticklabels([]), ax2.locator_params(axis='y', nbins=4)
ax2.yaxis.set_label_coords(-.07, 0.5)
ax3.set_ylabel('Fz (%s)' %Funits)
ax3.locator_params(axis='y', nbins=4)
ax3.yaxis.set_label_coords(-.07, 0.5)
ax3.set_xlabel('Time (s)')
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('COP (%s)' %COPunits)
ax5.set_xlabel('COPml (%s)' %COPunits)
ax5.set_ylabel('COPap (%s)' %COPunits)
ax1.plot(t, Fx), ax2.plot(t, Fy), ax3.plot(t, Fz)
ax4.plot(t, COPx, 'b', label='COP ap')
ax4.plot(t, COPy, 'r', label='COP ml')
ax4.yaxis.set_label_coords(-.1, 0.5)
ax4.legend(fontsize=12, loc='best', framealpha=.5)
ax5.plot(COPy, COPx)
ax5.locator_params(axis='both', nbins=5)
plt.suptitle('Ground reaction force data during quiet standing', fontsize=20, y=1)
plt.show()
plot_grf(grf)
rambling, trembling, iep0, inds0 = iep_decomp(t, Fx, COPx)
t2 = t[inds0[0]:inds0[-1]+1]
Fx2 = Fx[inds0[0]:inds0[-1]+1] - np.mean(Fx)
COPx2 = COPx[inds0[0]:inds0[-1]+1]
def plot_rambtremb(t, Fx, COPx, rambling, trembling, iep0, inds0):
Funits, Munits, COPunits = 'N', 'Nm', 'cm'
plt.figure(figsize=(12, 6))
gs1 = gridspec.GridSpec(4, 1)
gs1.update(bottom=0.01, top=0.96, hspace=0.12, wspace=.15)
ax1 = plt.subplot(gs1[0])
ax2 = plt.subplot(gs1[1])
ax3 = plt.subplot(gs1[2])
ax4 = plt.subplot(gs1[3])
ax1.set_ylabel('F (%s)' %Funits)
ax1.set_xticklabels([]), ax1.locator_params(axis='y', nbins=4)
ax1.yaxis.set_label_coords(-.05, 0.5)
ax2.set_ylabel('COP (%s)' %COPunits)
ax2.set_xticklabels([]), ax2.locator_params(axis='y', nbins=4)
ax2.yaxis.set_label_coords(-.05, 0.5)
ax3.set_ylabel('Rambling'), ax3.locator_params(axis='y', nbins=4)
ax3.set_xticklabels([]), ax3.locator_params(axis='y', nbins=4)
ax3.yaxis.set_label_coords(-.05, 0.5)
ax4.set_ylabel('Trembling'), ax4.locator_params(axis='y', nbins=4)
ax4.yaxis.set_label_coords(-.05, 0.5)
ax4.set_xlabel('Time (s)')
ax1.plot(t, Fx)
ax1.plot(t[iep0], Fx[iep0], 'ro', markersize=5)
ax2.plot(t, COPx, linewidth=3)
ax3.plot(t, rambling, linewidth=3)
ax3.plot(t, COPx, 'k', linewidth=.5)
ax4.plot(t, trembling, linewidth=3)
plt.suptitle('Rambling & Trembling decomposition during quiet standing', fontsize=18, y=1)
plt.show()
plot_rambtremb(t2, Fx2, COPx2, rambling, trembling, iep0, inds0)
Santos DA, Duarte M. (2016) A public data set of human balance evaluations. PeerJ4:e2648 https://doi.org/10.7717/peerj.2648.
The data set comprises signals from the force platform (raw data for the force, moments of forces, and centers of pressure) of 163 subjects plus one file with information about the subjects and balance conditions and the results of the other evaluations.
Subject’s balance was evaluated by posturography during standing still for 60 s in four different conditions where vision and the standing surface were manipulated.
path2 = r'./../../../X/BDB/'
fname = os.path.join(path2, 'BDSinfo.txt')
BDSinfo = pd.read_csv(fname, sep='\t', header=0, index_col=None, engine='c', encoding='utf-8')
print("Information of %s subjects loaded (%s rows, %s columns)."
%(len(pd.unique(BDSinfo.Subject)), BDSinfo.shape[0], BDSinfo.shape[1]))
Information of 163 subjects loaded (1930 rows, 64 columns).
fp = FloatProgress(min=0, max=len(BDSinfo.Trial)-1)
display(fp)
freq = 100
for i, fname in enumerate(BDSinfo.Trial):
filename = os.path.join(path2, fname + '.txt')
fp.description = 'Reading data from file %s (%s/%s)/n' %(os.path.basename(filename),
i+1, len(BDSinfo.Trial))
fp.value = i
grf = pd.read_csv(filename, delimiter='\t', skiprows=1, header=None,
names=['Time','Fx','Fy','Fz','Mx','My','Mz','COPx','COPy'],
engine='c')
t, Fx, Fy, Fz, Mx, My, Mz, COPx, COPy = [_ for _ in grf.values.T]
# rambling-trembling decomposition:
rambling, trembling, iep0, inds0 = iep_decomp(t, Fx, COPx)
t2 = t[inds0[0]:inds0[-1]+1]
Fx2 = Fx[inds0[0]:inds0[-1]+1] - np.mean(Fx)
COPx2 = COPx[inds0[0]:inds0[-1]+1]
# calculate some summary statistics of the trajectories:
COPsd, Rsd, Tsd = np.std(COPx2), np.std(rambling), np.std(trembling)
BDSinfo.loc[i, 'COP_sd'] = COPsd
BDSinfo.loc[i, 'Rambling_sd'] = Rsd
BDSinfo.loc[i, 'Trembling_sd'] = Tsd
BDSinfo.to_csv(os.path.join(path2, 'BDSinfoIEP.txt'), sep='\t',
encoding='utf-8',index=False)
print('Data from %d files were processed.' %len(BDSinfo.Trial))
FloatProgress(value=0.0, max=1929.0)
Data from 1930 files were processed.
fname = os.path.join(path2, 'BDSinfoIEP.txt')
BDSinfo = pd.read_csv(fname, sep='\t', header=0, index_col=None, engine='c', encoding='utf-8')
print("Information from %s files successfully loaded (total of %s subjects)."
%(len(BDSinfo), len(pd.unique(BDSinfo.Subject))))
Information from 1930 files successfully loaded (total of 163 subjects).
Because we have 3 trials per condition for each subject, let's take the median across trials for each subject as representative of the subject:
BDSinfo = BDSinfo.groupby(['Subject','Vision','Surface','Illness','Disability',\
'AgeGroup'], as_index=False).median()
print('%s subjects.' %len(pd.unique(BDSinfo.Subject)))
163 subjects.
We may choose to not consider the subjects with disabilities because their numbers are unbalanced in the age groups:
print('Before selection: %s subjects.' %len(pd.unique(BDSinfo.Subject)))
#BDSinfo = BDSinfo.loc[BDSinfo.Disability=='No'] # uncomment this line
display(BDSinfo.drop_duplicates(subset='Subject')\
[['AgeGroup', 'Subject']].groupby(['AgeGroup']).count())
print('After selection: %s subjects.' %len(pd.unique(BDSinfo.Subject)))
Before selection: 163 subjects.
Subject | |
---|---|
AgeGroup | |
Old | 76 |
Young | 87 |
After selection: 163 subjects.
pd.set_option('precision', 2)
BDSinfo.groupby(['AgeGroup', 'Vision', 'Surface'])\
[['COP_sd','Rambling_sd','Trembling_sd']].agg([np.mean, np.std])
COP_sd | Rambling_sd | Trembling_sd | ||||||
---|---|---|---|---|---|---|---|---|
mean | std | mean | std | mean | std | |||
AgeGroup | Vision | Surface | ||||||
Old | Closed | Firm | 0.48 | 0.21 | 0.42 | 0.20 | 0.18 | 0.10 |
Foam | 1.28 | 0.27 | 0.98 | 0.24 | 0.57 | 0.14 | ||
Open | Firm | 0.48 | 0.22 | 0.44 | 0.22 | 0.16 | 0.08 | |
Foam | 1.14 | 0.20 | 0.91 | 0.21 | 0.49 | 0.11 | ||
Young | Closed | Firm | 0.51 | 0.23 | 0.49 | 0.23 | 0.12 | 0.06 |
Foam | 1.01 | 0.21 | 0.83 | 0.21 | 0.40 | 0.14 | ||
Open | Firm | 0.45 | 0.19 | 0.43 | 0.20 | 0.11 | 0.05 | |
Foam | 0.96 | 0.25 | 0.86 | 0.26 | 0.31 | 0.12 |
def plot_summary(data):
g0 = sns.catplot(x='Vision', y='COP_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g0.set_xticklabels(''), g0.set_axis_labels('', 'COP (cm)')
g1 = sns.catplot(x='Vision', y='Rambling_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g1.set_axis_labels('', 'Rambling (cm)')
g1.set_titles('',''), g1.set_xticklabels('')
g2 = sns.catplot(x='Vision', y='Trembling_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g2.set_axis_labels('Vision', 'Trembling (cm)'), g2.set_titles('','')
plt.show()
plot_summary(BDSinfo)
It's possible to estimate the COG vertical projection (COGv or GL) from the COP displacement based on the inverted-pendulum model of the body on quiet standing.
See the notebook The inverted pendulum model of the human standing posture
Let's do some analysis with that.
from cogve import cogve
fig, ax = plt.subplots(1, 1, figsize=(12, 6))
cogv = cogve(COPx2, freq=100, mass=np.mean(Fz)/10, height=170, ax=ax, show=True)
cop_cogv = COPx2 - cogv
fname = os.path.join(path2, 'BDSinfo.txt')
BDSinfo = pd.read_csv(fname, sep='\t', header=0, index_col=None, engine='c',
encoding='utf-8')
print("Information of %s subjects loaded (%s rows, %s columns)."
%(len(pd.unique(BDSinfo.Subject)), BDSinfo.shape[0], BDSinfo.shape[1]))
display(BDSinfo.drop_duplicates(subset='Subject')\
[['AgeGroup', 'Subject']].groupby(['AgeGroup']).count())
#BDSinfo = BDSinfo.loc[BDSinfo.Disability=='No'] # uncomment this line
print('After selection: %s subjects.' %len(pd.unique(BDSinfo.Subject)))
Information of 163 subjects loaded (1930 rows, 64 columns).
Subject | |
---|---|
AgeGroup | |
Old | 76 |
Young | 87 |
After selection: 163 subjects.
fp = FloatProgress(min=0, max=len(BDSinfo.Trial)-1)
display(fp)
freq = 100
for i, fname in enumerate(BDSinfo.Trial):
filename = os.path.join(path2, fname + '.txt')
fp.description = 'Reading data from file %s (%s/%s)/n' %(os.path.basename(filename),
i+1, len(BDSinfo.Trial))
fp.value = i
grf = pd.read_csv(filename, delimiter='\t', skiprows=1, header=None,
names=['Time','Fx','Fy','Fz','Mx','My','Mz','COPx','COPy'],
engine='c')
t, Fx, Fy, Fz, Mx, My, Mz, COPx, COPy = [_ for _ in grf.values.T]
# COP-COGv decomposition:
t2 = t[inds0[0]:inds0[-1]+1]
Fx2 = Fx[inds0[0]:inds0[-1]+1] - np.mean(Fx)
COPx2 = COPx[inds0[0]:inds0[-1]+1]
# calculate some summary statistics of the trajectories:
cogv = cogve(COPx2, freq=100, mass=np.mean(Fz)/10, height=170, ax=ax, show=False)
copcogv = COPx2 - cogv
COPsd, COGVsd, COPCOGVsd = np.std(COPx2), np.std(cogv), np.std(copcogv)
BDSinfo.loc[i, 'COP2_sd'] = COPsd
BDSinfo.loc[i, 'COGv_sd'] = COGVsd
BDSinfo.loc[i, 'COPCOGv_sd'] = COPCOGVsd
BDSinfo.to_csv(os.path.join(path2, 'BDSinfoCOPCOGv.txt'), sep='\t',
encoding='utf-8', index=False)
print('Data from %d files were processed.' %len(BDSinfo.Trial))
FloatProgress(value=0.0, max=1929.0)
Data from 1930 files were processed.
def plot_summary2(data):
g0 = sns.catplot(x='Vision', y='COP2_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g0.set_xticklabels(''), g0.set_axis_labels('', 'COP (cm)')
g1 = sns.catplot(x='Vision', y='COGv_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g1.set_axis_labels('', 'COGv (cm)')
g1.set_titles('',''), g1.set_xticklabels('')
g2 = sns.catplot(x='Vision', y='COPCOGv_sd', hue='AgeGroup', order=['Open', 'Closed'],
data=data, estimator=np.mean, ci=95, col='Surface',
kind='point', dodge=True, sharey=False, height=3, aspect=1.6)
g2.set_axis_labels('Vision', 'COP-COGv (cm)'), g2.set_titles('','')
plt.show()
plot_summary2(BDSinfo)
from scipy import signal
freq = 100
t = np.arange(0, 1, .01)
w = 2*np.pi*1 # 1 Hz
rambling = np.sin(w*t)
trembling = 0.1*np.sin(10*w*t)
cop = rambling + trembling
rms = lambda x : np.sqrt(np.mean(x*x))
RMScop = rms(cop)
RMSrmb = rms(rambling)
RMStmb = rms(trembling)
# plot
fig, ax1 = plt.subplots(1, 1, figsize=(12, 6))
ax1.plot(t, cop, 'b', linewidth=3, label = 'COP: RMS=%.3f'%RMScop)
ax1.plot(t, rambling, 'g', linewidth=3, label = 'Rambling: RMS=%.3f'%RMSrmb)
ax1.plot(t, trembling, 'r', linewidth=3, label = 'Trembling: RMS=%.3f'%RMStmb)
ax1.legend(frameon=False, fontsize=18)
ax1.set_title("Fictitious COP, Rambling and Trembling and their 'amplitudes'",
fontsize=18)
ax1.set_xlabel("Time [s]")
ax1.set_ylabel("Amplitude")
plt.show()