The choice of the sampling interval and the selection of the order of the model in time series analysis are considered. Band limited (up to 15 Hz) random torque perturbations are applied to the human ankle joint. The applied torque input, the angular rotation output, and the electromyographic activity using surface electrodes from the extensor and flexor muscles of the ankle joint are recorded. Autoregressive moving average models are developed. A parameter constraining technique is applied to develop more reliable models. The asymptotic behavior of the system must be taken into account during parameter optimization to develop predictive models

Agarwal, G. C.

Goodarzi, S. M.

Gottlieb, G. L.

Oneill, W. D.

English

NASA Technical Reports Server

, N82-13698 -)23
COMPUTATIONAL PROBLEMS IN AUTOREGRESSIVE MOVING AVERAGE (AR_) MODELS
!
Gyan C. Agarwal, Shirin M. Goodarzi, William D. O'Neill and Gerald L. Gottlieb
College of Engineering I:i
|
University of Illinois at Chicago Circle
Chicago, Iii. 60680
and
Department of Physiology
Rush Hedical Center
Chicago, Iii. 60612 /
Abstract !i
In developing mathematical models of systems from a given input-output i
I.
data sequence, the choice of the sampling interval and the selection of the I_
order of the model in time-series analysis pose difficult problems. Band- [_
limited (up to i5 Hz) random torque perturbations were applied to the human If,
ankle joint. The applied torque input, the angular rotation output, and the
electromyographic activity using surface electrodes from the extensor and flexor
muscles of the ankle joint were recorded. Autoregressive moving average models
were developed. A parameter constraining technique is applied to develop more
reliable models. It is sho_ that the asymptotic behavior of the system must
be taken into account during parameter optimization to develop predictive models.
INTRODUCTION
In a series of previous papers (Agarwal and Gottlieb, 1977 a,b; Gottlieb
and Agarwal, 1978; Gottlieb, Agarwal, and Penn, 1978) we have attemped to
describe quantitatively the neuromuscluar system dynamics to applied sinusoidal
and band-limited gaussian torque perturbations. In these studies, the compliance
of the joint was calculated using Fourier series analysis for sinusoidal and
power spectral density methods for random perturbations. Although linear
. analysis methods were used, the system is known to be nonlinear and the parameter
= values such as the joint viscous and stiffness coefficients are functions of
the level of neuromusclar activity,
The purpose of the present paper is to apply time series analysis methods
-361- _,.$ !.'vr_rr: ...... ,,,., • ....
1982005792-358
https://ntrs.nasa.gov/search.jsp?R=19820005825 2020-03-21T15:59:10+00:00Z
rto study the input-output behavior of the neuromuscular system. The time series
method is very parsimonious in the use of parameters to represent the model
structure. Normalized residual criterion (NRC) will be used to estimate the
model order (For details of this method see Suen and Liu, 1977; Osafo-Charles
et. al., 1980).
Our previous analysis was limited to analysls of the angular rotation data
and calculation of Joint compliance. The electromyographic (EMG) data was not
analysed due to inherent difficulties in representing this output by linear
transfer functions. The time series approach allows nonlinear representations
as long as the model is linear in parameter space.
Dufresne, Soechting and Terzuolo (1978) used pseudo-random torque pulses
to study the human forearm response. They developed a mode] of the EMG in
terms of the lime position and its derivatives in the fo]lowing form:
•" (1)e_IG(t) = A O(t - d) + B _(t - d) + C e(t - d)
where A, B, and C are constant parameters and d is the time delay. They
found that the motor output depends primarily on the angular velocity of the
joint. Tile time delay was found to be about 47 msec.
In a subsquent study, Dufresne, Soechtlng, and Terzuolo (1979) used
different time delay parameters for position and its two derivatives. The
best estimates for the time delays were found to be 86 msec for position, 25 mse¢
for velocity, and 45 msec for acceleration. The physiological processes
associated with these varying delays are not clear. Soechting and Dufresne
(1980) found that the linear model given in equation (I) predicted 80% of the
EMC response.
Our analysis of the EMG using time series shows that the autoregressive
terms of the EMG are important and cannot be ignored as was done in the Dufresne
et. al (1978) model.
-362-
1982005792-359
METHODS
These experiments were done using normal human subjects. A subject sat
in a chair with the right foot strapped to a footplate which could rotate about
a horizontal, dorsal-plantar axis through the medial malleolus. The plate could _
be rotated by a DC torque motor. A band-limited gaussian ( 0-15 Hz ) signal was
prerecorded from a noise generator. These time-varying signals were superimposed
on a biasing mean motor torque level. The subject was instructed to try to
maintain a constant mean force against the bias torque of the motor so that the
ankle joint movement was nearly symmetrical with respect to the reference angle.
The input was applied for 30 sec or more and the data continuously recorded on
a digital tape.
The torque was measured hy a strain gauge bridge on the side arms of the
footplate. Angular rotation was measured by a continuous capacitive transducer.
The EMGs were recorded from disc surface electrodes taped over the bellies of
the soleus (SM) and the anterior tibial (TA) muscles. These were amplified
full-wave rectified and passed through an averaging filter (i0 msec averaging
time) b_fore recording. A computer generated the motor drive voltage at
a conversion rate of 250/s and digitized data on four input channels. The angle
and the torque signals were sampled at a rate of 250/s and filtered EMGs at
d rate of 500/s. The data analysis was done off-line using the Minitab 2
statistical software package on an IBM 370 computer. "'
-363-
1982005792-360
JThe Normalized Residual Criterion ,
Time-serles analysis can be extended to obtain discrete linear transfer
functions of systems having an input x(t) and output y(t). By x(t) and y(t) we
mean pairs of observations that are available at equispaced intervals of time. /
The behavior of the dynamic system can be adequately represented by the present
and past responses and the current and past inputs of the systems. We denote this
process as transfer function (TF) models (n,m) and write its equation as
y(t) = aO + alY(t - 1) _ "'. + anY(t - n) + 80x(t)+...+fl x(t-n)_v(t)
n (2)
In (2) the parameters to be estimated are a0,...,Sn, 80 ....,Bm, n, and
m. The time series v(t) is a random term measuring the difference between the
response y(t) and the variables used to expalin the time-series data. The
parameter a0 measures the mean output.
Equation (2) reduces to an autoregressive model (AR(n)) if x(t) is omitted
from the model, and reduces to a moving average model (}_(m)) if lags of y are
omitted. The following assumptions will be made concerning _(t) for a given
ouput time sequence y(t), t = [O,T],
I) El_(t)] = 0
2
2) E[,,(i),(j)] = ov61j (3)
where
{_ for i=j_ij = for i_j "_
3) T_'n.
rrom (2), - n
v(t) * y(t)-aO- Eaty(t - i) -E _jx(t- j),
t=l j-O
Define t - I, 2, "'.,T-n. (4)
1-n T-.
. ;jvj; and -llYII (5)
t-I t'l
-364-
1982005792-361
r: .L !
T
Note that in the discussion below V and Y are vectors such that
v(1) y(1)
(2) y()
V= Y=
v (T-n) y (T-n) (6)
Squaring (4) and normalizing by the total sum of squares, we have
= - - _ = E(n,m, (7)
IIYII 2 IIYII 2
and therefore
Since y(t), the data series, is deterministic, (8) can be rewritten as
From (6) we have
_:hich by assumption 2) in (3) reduces to
Substitution of (11) in (9) we have
(T- n) 7_
, - • (12)
L J I!YIi 2
and by assumption 3), (12) becomes
$ E _(n,m, T) : ToC (13)
lIYll_
The quantity _(n,m,T} depends cn n, m, and T and is proportlonal to the
normalized variance of the regression for a given n and m. If this ratio is
mininized over n and m, then the data fit as measured by the correlation coeffi-
cient p will be maximized. Note that
:_. -365-
L
1982005792-362
_= [1- ;iv];2]if2ii_i12 (m4_
or
where T, being a constant for the data, is omitted in the optimization procedure,
and _(n,m) is the minimum value for c(n,m). This optimization technique is the
so called the Normalized Residual Criterion.
-366-
1982005792-363
RESULTS
The mathematical modeling problem was considered in two separate parts.
For the first model the applied torque is the input and the resulting angular
rotation of the joint is the output of the system. For the second model, the
angular rotation (and its derivatives) is considered as the input to the
system and the resulting stretch reflex electromyographic activity is consi-
dered as the output. It should be emphasized that the angular rotation is
the net result of two torque inputs applied at the joint; one by the external
motor torque and the other muscle forces produced by the stretch reflex mech-
anism. These mechanisms are also responsible for a signi¢icant contribution
to the joint viscous and elastic properties. Figure I shows a sample of the
data at 4 msec sampling interval. The velocity was obtained by digital
differentiation.
Angle-Torque Model
Although the data was recorded for 30 seconds at each input (Agarwal and
Gottlieb, 1977b), this method does not require such long data records which
w_uld also use too much computer time. The time series analysis was done using
only two-seconds of the data record. (The first two-seconds of the data were not
used to a11ow the turn-on transients to die out).
The values of [(n,m'Pwere computed for a given data record and then plotted
against different values of n (see Figure 2). T;le data sampling interval in
this case is 4 msec. Thi_ analysis clearly indicates that n = 2 and m • 0 is
adequate to model this data. The same data was analyzed again using the sampling
intervals of 12, 20, 40, and 60 msec. Figure 3 shows the c(n,m) values for the
sampling interval of 20 msec. Note that the minimum value ol the normalized
residual is about 60 times of that in the first case. For 40 and 60 msec sampling
c(n,m) did not reach as}_ptotic values even for model order of (8,8). The norm-
I alized resldual_ at 12 and 20 msec sampling implied a model order of (3.1).
-367-
1982005792-364
Figure 4 shows the actual angular rotation data (2 to 4 sec interval used
in this analysis ), the regression fit and the predicted output using 4 msec
sampling and model order of (3,!). The regression fit is obtained by using the
equation
• 0(t) = a 0 + a I 0(t - 1) + a 2 0(t - 2) + a 3 0(t - 3) + 80T(t) + 81_(t - '"
(16)
The error between the actual data and the regression fit is nearly zero.
The correlattop coefficient is p • 0.999. However, when this mode, _ sed to
predict the output using the first three data output values as the inxtial con-
ditions, the predicted output is a poor approximation of the actual data (see
Figure 4). Figure 5 shows the observed angle and the predicted model values for
model ordtrs of (3,1), (7,1), (9,1), and (14,1). Even the fourteenth order
model is not able to adequately reproduce the data sequence. These models are
not able to capture the steady state (or long term) behavior of the system.
Osafo-Charles, et al., (1980) showed ti_at to develor ke_ter predictive models,
the TF(n,m) models must be constrained to incorporate the steady state response
of the system.
Constrained Model
Consider the estimated model given by equation (16). Under conditions
of equllibrlum
o(t) - o(t - 1) = o(t - 2) - o(t - 3_ = 0e
,uld
T(t) ,, T(t - 1) - Te
_'here O and T are the steady state response and input respectlvelv. At physicale e "
and statistical equilibrium, with 0(t) = 0 e and T(t) = Te, equaition (16) became
_e = al 0e + a2 0e + a3 0e + _0 Te + dl Tc (17)
or
Oe t_0 "_ t31
• --_. g (18)
T e 1 - a 1- a 2 - a3
-:'8-
1982005792-365
where (18) expres:;es the steady state gain in terms of the parameters of the model.
The value of g was approximated by the slope of the curve of torque vs. angular
rotation in the relaxed ankle during slnusoldal oscillatlon at 0.i H (Gottlleb
z
and Agar_a!, 1978).
• For 8e = g Te to be true, we must have
80 " 8(1 - al - a2 - a3)- 81 (19)
From equctlons (19) and (16), we get
+ a3 (t - 3) - g T(t
+ 61 _(t - i) - T (t_ (20)
Regression analysis is used again to estimate the parameters al, a2,
a3 and 61 for a given value of gain 8. 60 is then obtained using (19). Figure
6 shows the output angle and predicted model response for a constrained model
with gains of g - 6.5, 7.5, and 8.5. The gain value of 8.5 was considered to
provide the best fit in terms of the minimum estimated standard devlatlon of the
regression.
The transfer function for the unconstrained model is:
-l
H(z) - 0.00239 - 0.00024 z
-) -_ -3 (21)
1 - 2.678 z + 2.399 z - 0.7191 z
For the constrained model with a slope of 8.5, the transfer function is:
-1
tt(z) ," 0.00238 + 0.00017 z
; (22)
-1 -2 -3
. 1 - 2.731 z + 2.503 z - 0.7717 z
-tF
,e ,-
-y
Z
-369-
I
1982005792-366
EMG Model
Our efforts to model EHG as a function of either the angular rotation
or the velocity or a combination of both were not successful. As was noted I
by Dufresne et al. (1978), the velocity of rotation is the most significant
input due to splndle properties (Matthews, 1972), However only those
components of velocity which stretch the spindle contribute to the EHG of ,_
the stretched muscle. (The splndle is silent during shortening). Therefore, J!o
J
a new velocity slgnal represnetlng only the stretching velocities was
d_fined as:
, • J
ed (t) = 0 (t) if 0 _> 0 _j
• != O if 0 < 0 (23)
The normalized residual analysis indicated a model order of (4,1) using
soleus E,_IGas the ouput and 0d as the input signal. The predicted output of
the unconstrained model and its comparison with the actual EMG slgnal is
shown in Figure 7. Since the EMG signal is a full-wave rectified and filtered
fusing an averaging filter) slgnal, it has only negative values (because of
negative filter gain). The predicted value of EMG ks a poor approxi_Jtlon of
the data.
A const_ai,,ed model was developed u_ing a similar approach as outlined
earlier. Figure (8) sho_'s tle predicted EHG and the actual data at three
values _" the gl, tn parameter. The gain of -0.005 was considered to give the
most appropriate fit.
For the constrained model with a slope of -0.005, the transfer f'anction
is:
-1
H(z) = DIG - -.004164 + .00314 z (24)
_d 1 - 1.19 z-I + ,6685 z-: - .2947 z-3 + .02104 z"_
-3?O-
1982005792-367
.......... T
CONCLUSIONS
The time series approach is a powerful and versatile technique in
developing time domain models from a given input-output data sequence. Norm-
alized residual criterion allows effective prediction of the model order.
Models developed in this manner may be satisfactory, but may not be good
predic*Ive models. It is recommended that constrained parameter modeling
which allows incorporating the steady-state behavior be used to obtain
better predictive models.
Acknowled_ nt
This work was supported in part by National Science Foundation grant
ENG-7608754 and grants from the National Institutes of Health NS-12877 and
NS-00196.
L
-311-
1982005792-368
REFERENCES
i. G. C. Agarwal and G. L. Gottlieb. Oscillation of the human ankle joint
in response to appled sinsoidal torque on the foot. .1. Physiology
(London), 268:151-176, 1977a.
L
2. G. C. Agarwal and G. L. Gottlieb. Compliance of the human ankle joint.
, ASME J. of Biomechanical Engineering, 99:166-170, 1977b.
3. J. R. Dufrense, J. F. Soechting, and C. A. Terzuolo. Electromyographic
response to pseudo-random torque disturbances of human forearm position. {
Neuroscience, 3:1213-1216, 1978.
4. J. R. Dufrense, J. F. Soechting, and C. A. Terzuolo. Reflex motor out-
put to torque pulses in man: Identification of short - and long- latency
loops with individual feedback parameters. Neuroscience, 4:1493-1500,
1979.
5. G. L. Gottliph, G. C. Agarwal, and R. Penn. Sinusoidal o_cillation of the
ankle as a means of evaluating the spastic patient. J. Neurology,
Neurosurgery, and Psychiatry, 41: 32-39, 1978. _"
6. G. L. Gottlieb and G. C. Agarwal. Dependence of human ankle compliance
on joint angle. J. Biomechanics, 11:177-181, 1978.
7. P. B. C. Matthews. Mammalian Mascle Receptors and Their Central Actions,
Williams & Wilkins Comp., Baltimore, 1972.
8. F. Osafo-Charles, G. C. Agarwal, W. D. O'Neill, and G. L. Gottlieb.
Application of time-series modeling to human operator dynamics. IEEE
Trans. Systems, Man_ and Cybernetics, SMC-10:849-860, 1980.
9, J. F. Soechting and J. R. Dufresne. An eva_uatiun of nonlinearities in
the motor output response to applied torque perturbations in man.
Biol. Cybernetics, 36:63-71, ]980.
I0. L. C. Suen and R. Liu. A normalized residual criterion for the deter-
mination of structure of linear systems. Policy Analysis & Information
S_'stems, 1:1-22, 1977. •
-372=
1982005792-369
_ _-'--- | _ |' L ''i&.| "i ,J Ill | '
_8 - !.
-16
I J_INT .N(;L.E $.I'tPLED .T 4'HS 114[ i _ _" i" e _ I I i i/
-t0 -
I _,_so,,s_,,_LEo,T, _.s ]
0.0 ........
-2.04
/ t'ELOCITY S.I'IPLED _'T 4 MS 1400[- ' r-T--F--_-_-F T _ i ,_
-400L._. t i I I [ I 1 I __ ":
2.o =A 2.e s._ 3.e 4.o
•:,, TIME(SEC)
Figure 1
i,, -373-
1982005792-370
SRIIPLEORT 4 IqS
0.016 1 I i I I i I I I
(/) 0.014
,,I - ,0
<
:) 0.012 _ m .0 TO 9 _
a_(/)
lU 0.010-
Q 0.008
N
_ o.ooe
a, 0.004
o
Z 0.002
o I I i i I I I4 6 8
NUIqBEROFLRGS-OUTPUT
Figure 2
Z2[ SRMPLEORTZOHS
f_l 0 I I I I I I | I I I I I I I I
_ 0.18
j 0 16
CC 0 14
0 12
_N
o,o i'_ 0 08-ZOo 06
0 04 f I I I I ,| I I I I I I I I, I
2 3 4 5 6 ? O
HUNGEROFLRGS-OUTPUT
t'igure 3
-3_4-
1982005792-371
....
j ll'w,j. _l lilt I ,jl,lL o.
•f / _IP
_ - - -
_ o
_,- _p
2"ii' ' I O,
I I !
i (D3(]) _rl_)NV J,NIOI"
!_ -375-
1982005792-372
, ., . ._ ........ _.,,,._,,.m,=_=,_.,._,_,,,_'_ "_'''*=vg'_m_¢_'_'_" _--
¢
[ COHSTRRINEORODEL(SLOPE'6.5) (Nz3,M,I} ] :"
l co,sT,,,._D.OOEL_S,O,-S.5)-- J
,!L , ', _", " , - ,-i
2.0 2.4 2.8 3.2 3._ 4.0
TIME(SEC)
Figure 6
:g2.0
-4.0
2.0 2.4 2.8 3.2 3.6 4.0
TIME(SEC)
Figure 7
-376-
1982005792-373
QA
a
COHSTR_IHEOflOOEL(K=-_S) N-4,fl=l
_ -2.0 -
-4,o - I:
W
L
CONSTRAINEOMODEL(K--.O06) N-4,nsl
-4.o , , _ i i i i i
2.0 2.4 2.8 3.2 3.6 4.0
! TIMEISEC)
_, Figure 8
r ;
i. -377-
1982005792-374


Computational problems in autoregressive moving average (ARMA) models

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server