Functional Separation of EMG Signals via ARMA Identification Methods for Prosthesis Control Purposes

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS, VOL.

252

(21 [3] [4]

[5] [6] [7]

[81

D. T. McRuer, D. Graham, E. Krendel, and W. Reisener, "Human pilot dynamics in compensatory systems," AFFDL-TR65-15, July, 1965. R. W. Roig, "A comparison between human operator and optimum linear controller RMS-error performance," IEEE Trans. Hum. Factors Electron., vol. HFE-3, pp. 18-21, Mar. 1962. R. E. Kalman, "When is a linear control system optimal?" J. Basic Eng., ser. D, vol. 86, pp. 61-60, 1964. J. D. Dillow, "The 'paper pilot'-A digital computer program to predict pilot rating for the hover task," AFFDL-TR-70-40, Mar. 1971. R. L. Stapleford, "Application of parameter optimization to stability augmentation design," in Proc. 6th Ann. Conf. Manual Control, AD-709-381, Apr. 1970. T. L. Hollis, "Optimal selection of stability augmentation parameters to reduce the pilot ratings for the pitch tracking task," M.Sc. thesis, Air Force Institute of Technology, Ohio, June 1971. D. L. Kleinman, S. Baron, and W. H. Levison, "An optimal control model of human response," Automatica, vol. 6, pp.

357-369,1970.

[9] R. W. Obermayer and F. A. Muckler, "On the inverse optimal control problem in manual control systems," in IEEE Int. Conv. Rec., Part 6, pp. 153-165, 1965. [10] D. T. McRuer and E. S. Krendel, "The human operator as a

[11] [12] [13]

[14] [15] [16] [17] [18]

SMC-5, NO. 2, MARCH 1975

servo system element," J. Franklin Inst., vol. 267, pp. 511-536, June 1959. R. D. Murphy and K. S. Narendra, "Design of helicopter stabilization systems using optimal control theory," J. Aircraft, vol. 6, pp. 129-136, Mar. 1969. A. W. Starr and Y. C. Ho, "Non zero-sum differential games," J. Optimiz. Theory Appl., vol. 3, no. 3, pp. 184-206, 1969. W. S. Levine and M. Athans, "On the determination of the optimal constant output feedback gains for linear multivariable systems," IEEE Trans. Automat. Contr., vol. AC-15, pp. 44 48, Feb. 1970. B. D. 0. Anderson and J. B. Moore, Linear Optimal Control. New York: Prentice-Hall, 1971. J. C. Willems, "Least squares stationary optimal control and the algebraic Riccati equation," IEEE Trans. Automat. Contr., vol. AC-16, pp. 621-634, Dec. 1971. K. W. Anderson, "The numerical solution of the time invariant matrix Riccati equation," Dep. Elec. Eng., Univ. Queensland, Aust., CGR72/2, Oct. 1972. D. Rothschild and A. Jameson, "Comparison of four numerical algorithms for solving the Lyapunov Matrix equation," Int. J. Contr., vol. 11, no. 2, pp. 181-198, 1970. E. W. Vinje and D. P. Miller, "An analysis of pilot adaptation in a simulated VTOL hovering task," in Proc. 4th NASA Conf. Manual Control, NASA-SP-192, pp. 95-118, Mar. 1968.

Functional Separation of BMG Signals via ARMA Identification Methods for Prosthesis Control Purposes DANIEL GRAUPE, MEMBER, IEEE, AND WILLIAM K. CLINE, STUDENT MEMBER, IEEE Abstract-Multifunctional control of artificial limbs via electromyo- other interacting EMG signals and of recognizing and graphic (EMG) actuation requires means for reliably recognizing or separating the EMG signals related to different limb distinguishing between the various functions on the basis of the recorded functions. For adequate and reliable filtering and recogniEMG data. Furthermore, constraints of weight, cost, and computation time on practical prosthesis application must be satisfied. An approach to tion, a rigorous statistical analysis of EMG signals IS the aforementioned recognition problem is given in terms of deriving essential. However, the EMG treatment that is currently a fast parametric-recognition algorithm whereby the autoregressive- available in the literature is mostly restricted to ad hoc moving-average (ARMA) parameters and the Kalman filter parameters of the EMG time series are identified. It is shown that the resulting identified parameters yield sufficient information to discriminate between a small number of upper extremity functions. Problems involved in practical prosthesis control via the present approach and problems of hardware realization are discussed to illustrate the validity of the ap-

methods of spectral density and correlation evaluations

[9]

[1]-tm and therefore doesnot make optimal or at least optimal-linear use of the information content of the signals. l

Hence EMG prosthesis control has been essentially of single-function nature, and the resulting multifunctional control has been of ad hoc nature, employing some kind of proach. amplitude-level coding [5]. This in turn has usually required I. INTRODUCTION level training of the amputee to actuate multifunctional A AJOR problem in multifunctional prosthesis . control to affect reliability and amputee's comfort. Methods control via electromyographic (EMG) signals is that of pattern recognition of EMG signals were recently u or a [1 T of the filtering of the EMG signals from noise and from po s electrodes and complex memory and computation regarding Manuscript received March 21, 1974; revised August 27, 1974. This ad hoc and usually predetermined tasks to again affect work was supported in part by the Veterans Administration, Depart- .. ment of Medicine and Surgery, Prosthetic and Sensory Aids Service, reliability, simplicity, comfort, and speed. Furthermore, these again used only low frequency information from the under contract V1O1(134)P.i 11. MAJOR

problem

in

multifunctionl

piaMGsgasndwrthsotpimlrevn D. Graupe iS with the Department of Electrical Engineering, 'E Gsgasadwr hsntotmlo vnotml Colorado State University, Fort Collins, Colo. 80521. W. K. Cline was with the Department of Electrical Engineering, linear in terms of information utilization. Noting the recent Colorado State University, Fort Collins, Colo. He is now with Bell adncsiaplngrooutmeeisaayisehiqs avne napyn gru lesre nlsstcnqe Laboratories, Naperville, Ill.

253

GRAUPE AND CLINE: PROSTHESIS CONTROL

to electroencephalogram (EEG) analysis [11]-[13], a similar but non-ad-hoc philosophy has been proposed by the present author for application to EMG signals that utilizes the information in the signal in an optimal-linear fashion. The latter advances coupled with the enormous advances in microcomputer and microprocessor hardware, form the basis for this work and for its realizability with today's hardware. The present approach to EMG analysis recognizes the nonstationary and nonlinear nature of the EMG signals. Furthermore, it recognizes the practical constraints in any realization of filtering and identification or recognition, when employed in conjunction with a realistic prosthesis, the constraints being those of time available for processing and of cost and weight of computational hardware to be carried by the amputee (say, in his pocket). These constraints must and do certainly lead to compromises and to simplifying assumptions in any analysis, that may detract from the rigor of the analysis. However, by using microcomputer hardware, an adequate and reliable solution based on rigorous analysis is still realizable as long as the number of degrees of freedom considered is small. (An incorporation of the latter system with limb control as in [18]-[19] thus leads to considerable prosthesis maneuverability via EMG actuation, while causing hardly any mental strain in terms of an amputee's need to learn complex actuation combinations, etc.) The present analysis has been based on the nature of time series in that their pattern can be parameterized into a finite set of parameters of a linear signal model [14] that form a reduced minimal set as compared with the (almost) infinitum of values of the pattern itself. These parameters must not all be stationary or unique per each function for function separability. However, if at least some of these parameters are such that their range of variation for a given limb function does not overlap with that of other functions, then such separation is possible, as is shown later to be the case in all our surface-electrode EMG recordings (over 250

records). Specifically, and noting the preceding constraints, we have decided to examine minimum-order autoregressivemoving-average (ARMA) models of stationary time series [15], [20]-[24]. Although the EMG signal is not a stationary one, we have shown (see Section V) that the signal is sufficiently stationary per each of the limb functions considered to yield ARMA parameters whose range of variation with time is adequately small to facilitate discrimination of these limb functions. The choice of a minimum-parameter linear ARMA model becomes obvious since it yields the required function discrimination via a minimal-parameter

structure. Also, the linearity of our approach is justified noting that the identification of a nonlinear or of a nonstationary model is impossible in practice (especially since for practical actuation of limb functions, the identification and the recognition of a function must be completed within about 0.1 s). Furthermore, the error between the aforementioned model output and the actual data converges to a white-noise sequence that contains no information in its

first and second moments, such that all linearly retrievable information is contained in the parameters of the model. We comment that the parameters just considered need not be similar for similar functions for different patients. To overcome this, complete offline identification for establishing parameter ranges must be made for each patient prior to connecting his prosthesis controller to the identification and filtering microcomputer that is to execute the recognition. II. EXPERIMENTS For the present analysis, EMG signals have been recorded at three different locations (Li,L2,L3) along the biceps, triceps, and forearm of three volunteers. The signals relate to different limb functions: Fl, complete rest; F2, load of 1 lb on hand during arm bending; F3, grasp action; and F4, same as F2 but 5-lb load (see tables in Section IV). More than 250 records in total were considered. We comment that although an amputee cannot perform these functions (without a prosthesis), he can produce or can be trained to produce muscular tensions and flexions (above the amputation) that are related to these functions at appropriate electrode locations such that the preceding signals are available. The sampling rate for the EMG when fed to the computer for our experiments was 2000 samples/s, which is about the upper sampling rate of interest for EMG analysis [7]. However, no hardware or software limitations will hinder considering even a higher sampling rate (say, 10 000 samples/s). 1. PURE AUTOREGRESIVE (AR) AND ARMA III. EMG SIGNALS ANALYSIS FOR EMG SIGNALS Our analysis is concerned with identifying the parameters of the EMG patterns, as previously discussed. Hence, our analysis becomes that of identifying a finite number of parameters in order to detect some parameters that are stationary or that vary inside narrow ranges. Our resultant computations have shown that identification can be complete in less than 0.1 s, namely within an adequate time interval such that no unreasonable delays are introduced into the actuation of limb movements. The aforementioned computation time per pattern, namely the time required prior to limb function determination, also implies that we can cope with the lowest relevant EMG frequencies, which are of about 30 Hz [7]. The constraints mentioned in Section I clearly indicate that if a linear analysis can be justified, then it should be used, since otherwise computation will be too lengthy and complex to enable practical prosthesis applications. We note that a linear model is optimal only for a Gaussian process, which is not our case. However, stationary time series can always be modeled in a linear-optimal fashion as linear processes input with white noise [25], such that the identification error converges to a white-noise process, that is, contains no information in its first two moments. Since we are interested in piecewise-stationary features, such a linear model can be derived. Furthermore, if it will FoR

254

IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS, MARCH 1975

1.0

r-_

__

_

The preceding autocorrelation of 0k has been computed whiteness, as illustrated in Fig. 1. This check, and computations of var (0k), namely of the variance of the prediction error, were employed to establish the adequacy of the finite order used for the AR model, as was a check of the change in the estimates a1 to an. of a1 to an for assuming s = n and s = r + I for various r. Once the AR parameter estimates a1 as were available, the corresponding ARMA parameters were derived, noting the minimum-parameter character of the ARMA model, which is given by to check its

*.8

-------------t

-----

.6-

t 7

.4

4-I

.2-A

0.0-

-

X t t t t t tt i -

-~ ~~ ~ wi- - -

~

~

-

where

~

~

~

~

(B)yk= O(B)Ok

(5)

(~~B) 0(B) = a(B)

(6)

a(B) being the pure AR parameters of (1) and (2), and 2 - - - - - - - - ;(B),0(B) denoting the AR and the moving-average (MA) -.2 0 42 6 8 10 12 14 16 18 20 parameters of the mixed ARMA model, respectively, 4(B) i (DELAY) -.being of order n and 0(B) of order m. Observing (6), properties of polynomial division [24] Fig. 1. Residual autocorrelation function. were used to determine the ARMA parameters and order without any a priori knowledge and without any recursive eventually (and it actually does) yield separation of limb computation. The latter derivation, which is based on functions, then there is no point in going to a nonlinear solving an algebraic matrix equation [24], is extremely analysis, which is too lengthy anyway and computationally fast and is consistent for consistent AR parameter estimates too complex to be applicable to real life prostheses. Our (or otherwise bounded), by theorem 2.3.3. of Lukacs [26]. approach is thus in terms of identifying a pure (theoretically It is noted that the derivation of the ARMA parameters is infinite) AR model for the EMG sequence, from which the necessary only if a Kalman filter model is to be obtained, minimum-parameter ARMA model, and subsequently the for the purpose of filtering non-EMG noises. Otherwise, this derivation need not be performed, since wild order Kalman filter model, can be obtained, as follows, Let Yk denote the measured EMG time series (k = variations are not likely, and the first m + n AR param,l,2, * , denoting the interval). Noting [14], the whitening eters can be shown [22] to yield all the information on the filter for Yk of any distribution is a linear pure AR filter pattern. (This derivation must in fact be made once only, which is in theory, of infinite order, but whose parameters as a calibration for determining m and n.) form a convergent sequence for bounded Yk Since this is The preceding procedure is applied to all EMG recording our case, it is sufficient to use three to ten AR terms to locations such that for location [L(u = 1,2,- ,4), Yk of obtain a white residual 'k, as follows: (l)-(3) becomes Yk(0); 0(B), 0(B), a(B) become sO(B), ag(B), respectively; and Wk becomes Wk The + Wk Yk = (1) aiYk-i complete EMG model thus becomes a multi-output model, i= 1 whose outputs are spatially distributed. The treatment of a,being estimates of the AR parameters, s being the ordereach output is, however, via a single-input single-output of AR filter, and Wk being the white-noise residual. Equation model, for fast computation and convergence. We note that (1) can be written in operator form as a multi-input multi-output treatment, the theory of which is (2) outlined in [23], is certainly possible, but it is far lengthier. a(B)yk = Wk beiga,By = Yk- . The param- It is therefore not recommended as long as the single-input being operator, namely A d eters ai are identifiable as in [14, ch. 12] and in [24], using single-output model applied to all y(") facilitates adequate resolution of functions, and noting the constraints of Section a recursive least squares algorithm such that I. We shall return to this point in Section V. For results and A [al,,,islk their discussion, see Section IV. * ak = ak_ 1 + PkYk(yk - ak-I Yk) Yk-sII (3) IV. DERIVATION OF KALMAN FILTER MODELS FOR Yk _ [Yk-l1, lyT 4) EMG SIGNALS kI = 'k-1 - k-1 kG + Yk PY1) l Yk p (In Section III the identification of the ARMA model for f()()cnan The latter identification can be shown to be consistent FGsgaswsdsusd oee,Y (convergent in probability) for the correct AR order [15],EGsGaswsdsusdHwvr expressed as ^ , otherwise . 1 . an error that can be bounded' ~~~~~~~~non-EMGmeasurement noise. This can7 bef()()cnan [24] and to yield as required by appropriate selection of s [21]. Yk =Xk +flk (7) -

g

4,(B),

255

GRAUPE AND CLINE: PROSTHESIS CONTROL

TABLE I SEVERAL IDENTIFIED AR PARAMETERS FOR EMG SIGNAL FOR LI, F3 Autoregressive Parameters

a3

a2

a

a5

a4

a7

a6

a8

a1

a9

1.97 -1.03

.041

s=5

1.90 -.623

-.447

-.013

.166

s=10

1.88 -.569

-.431

-.055

.035

.175

-.076

.026

.090

.084

s=15

1.88 -.569

-.422

-.063

.037

.175

-.089

.033

.084

-.057

.

.

.

.

.

Li, Fl Ll, LI, F2 F3 Fl L2, F1

-. 6750

-1.772

Ll,

-1.922 -1.922

L2, F2

-

L2, F3

-

.7686

.6009

L2, -1 . 60 76 F3

- .2919 + .8235 + .9505 ++ .2009 - .2764

.9505

-.676 +1.10 .210

¢3

--

--

01

-.1117

-.320 -.3200

3

+.1505 +.1150

-.3937 +.2071

+.20021 +002170 +.0170

+.0298

+.3174

+.0357 +.0357

0217 +.0217

-3 _3

o2

1 x ¢12w

--

4.452

--

5.436

2.570 2574

Li, F3 Li, L2, Fl

7.365

L2, F3 L2, F2

1.250

--

1

6.839

-.0210

+.083 +.0083

Li,

Fl

F2

L,F

- .5296 -1.794 -1.901 -1.490 140 --.1660 .7925

TABLE IV

MODEL OF (3); ¢2

Li, Fl Li, F2 Li, F3 L2, Fl

L2, F2 L2, F3

- .8923 -1.801 -1.906 -1.377 - .8812 -1.787

-.0809 +.8517 +.9354 +.3937 -.0248 +.9486

p3. --

--

--

--

-.0816 -.1168

s

.. -.006

--

-.073

-.083

.104

.033

TABLE III MODEL OF (3); s = 5 02

-- +0434 0343 --+.i592

a15

a4

--..--

TABLE II MODEL OF (3); s = 3 ¢1 ¢2

a3

.

s=3

.-.

a12

a 11

3

-.4325 +.8445 +.9308

+.5040 +54

-.0630 +.8145

----

+.2991 -.1313 -.4386 +-1942 +.00073 +.3082 -.5775 +-0659

TABLE V

8I

e2

-.0700 -.4635 -.0303 -.4764 +.0506 -.0152

-.0817 +.2023

+.3111

+.0698

+.0412 -.4562

2

^

e3 Cw --

----

+.2299 -.2703

x

-

A

$1

10 3

4.388 1.213 2.152

5.145 6.485 6.991

nk being non-EMG measurement noise. Due to the nearly deterministic form of electrocardiographs (ECG), the ECG interference can be filtered out by currently available ad hoc filters and is considered not to be white noise. The assumption of nk being white noise can sometimes be relaxed. However, it is a common and convenient assumption for most measurement noises. Consequently, and since the results available at this time are related to this assumption, we shall restrict the present discussion to nk being white noise. Hence, substituting Ykfrom (7) into (5), we derive [16] (8) b(B)xk -= lfr(B)u E[uknl] -O V,V k,l 1] , o(B)xk= tk(B)uk ,kil (8) Uk being white noise, qk(B) being of order m < n, such that (1) and (5) yield (a2 denoting variances of the respective

Li, Fl Li, F2 Li, F3

L2, Fl L2, F2 L2, F3

- .8641 -1.807 -1.906 -1.502 - .9123 -1.716

$2 -.1074 +.8566 +.9358

+.5152

-.0240 +.8720

w

xlo4.452 1.240 2.232 5.315

--

----

-57 .69 _ .1 ~ 6.568 -.1313 +.1407 +.0835 -.0783 +.1132 -.3616 +.2015 -.2303 7.071

MODEL OF (3);

= 10

3

2

3

103

s

=

15

'3

^

61

--

-.0428

-__

--

-.0517 -.0934

^

^

62

-.0868 -.4834 +.1973 -.0318 +.3100 -.6185 +.0707 +.0192 +.0126 +.0597 -.4019

2

^

63 awx --

----

+.2340 -.2592

10

-3

4.365 1.183 2.144 5.088 6.459 6.915

space form required for Kalman filtering is straightforward [14], and all the required parameters and covariances of the Kalman filter are available. Hence, linear-optimal (and in the Gaussian case, optimal) filtering is possible. However, since we are concerned with function separation rather than with filtering, results as in Section V will indicate that Kalman filtering is not essential for our problem.

EMG data recorded at various arm locations (biceps, triceps, forearm, etc.) have been employed to derive the models of (2) and (5). The identified models obtained are tabulated in Tables I-V for various assumed AR orders s. In these tables L1,L2,** denote locations 1,2