IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
417
A Nonstationary Model for the Electromyogram EDWARD -SHWEDYK, MEMBER, IEEE, R. BALASUBRAMANIAN, AND R. N. SCOTT, SENIOR MEMBER, IEEE
Abstract-A theoretical model of the electromyographic (EMG) signal has been developed. In the model, the neural pulse train inputs were considered to be point processes which passed through linear, timeinvariant systems that represented the respective motor unit action potential. The outputs were then summed to produce the EMG. It was assumed, that in the production of muscle force, the controlled parameter was the number of active motor units, n (t). The model then showed that the EMG can be represented as an amplitude modulation process of the form EMG = [Kn (t)J 1/2 w(t) with the stochastic process, w(t), having the spectral and probability characteristics of the EMG during a constant contraction. Various assumptions made in the model development have been verified by experiments.
\Neural Pulse Train
Ae(t) n
Fig.
1.
n
Model for EMG signal generation. The number of active motor
units are n (t). The different motor unit action potentials are repreINTRODUCTION sented by the impulse responses h1 (t),-- hn(t). Coefficients ki represent the respective weighting constants. Characters with underCONTROL of skeletal muscle force is achieved by neural bars are boldface in text. signals, transmitted via motorneurons to the muscle, from the anterior horns of the spinal cord. The motorneuron, along with the muscle fibers it supplies, is called a motor unit and is do not achieve this ideal and design of them has been based, in the basic muscle unit under voluntary control. The muscle part, on intuition and convenience rather than on a fundamenforce produced is controlled primarily by (i) varying the num- tal knowledge of the EMG and its characteristics. This was to ber of activated motor units and (ii) varying the rate of activa- a large extent due to the lack of a suitable model for EMG gention of an individual motor unit. Though both factors are eration. Properties of a particular class of EMG processors present, (i) is generally considered to be the predominant have recently been described(9). Theoretical analysis of the one.(0) processors was based on an amplitude modulation model Muscle force generation is also accompanied by electrical of EMG signal generation. The resemblance of EMG to signal activity which, detected by surface electrodes, produces an amplitude modulation process has been pointed out the electromyogram (EMG). The EMG has a complex tempo- previously(1O0 1 1) ral pattern and is best described as a random signal which may This paper describes the development of an analytical model be modeled as a weighted sum of the action potentials of the of EMG generation. The development shows the restrictions active -motor units. Therefore, its statistical properties are de- and assumptions under which the EMG may be considered to pendent on the number of active motor units and their rate of be an amplitude modulation process. In the model, the moduactivation. lation signal is the number of active motor units while the carA functional relationship thus exists between the EMG and rier is a random signal whose statistical properties are those of muscle force as has been shown in several studies(2,3 ,4). This the EMG during an isometric contraction. Experimental evirelationship has been exploited to use the EMG for diag- dence is presented to justify the assumptions made in the nosis(5), analysis of skill(6) and gait studies"7. Further, since model development. The model and experiments bring into muscle force can be voluntarily controlled and hence so can focus the EMG nonstationarity and the difficulties in processthe EMG, the EMG has been used to control powered orthotic/ ing it. The paper concludes by suggesting a possible further diprosthetic devices(8). rection for research into EMG processors. In the above applications, processors are required which ideally transform the EMG into a noiseless output signal monoII. DEVELOPMENT OF THE MODEL tonically related to muscle contraction or equivalently to the Fig. 1 depicts the generation of EMG in a block diagram number of active motor units. Practical processors of course form. The neural pulse trains go through different linear, time-invariant systems, characterized by impulse responses Manuscript received July 17, 1975; revised July 19, 1976. This work was performed at the Bio-Engineering Institute, University of New h1 (t) h.(t), to generate the motor unit trains which are Brunswick, Fredericton, N.B., Canada, and was supported in part summed to produce the EMG. The weighting coefficients, by grants from Health and Welfare, Canada, and the National Research ki, are assumed to be random and a function of electrode Council of Canada. E. Shwedyk is with the Department of Electrical Engineering, Uni- placement. It is also assumed that the motor unit trains are versity of Manitoba, Winnipeg, Man., Canada. uncorrelated and, that each train is a renewal point process R. Balasubramanian and R. N. Scott are with the Department of Electrical Engineering and the Bio-Engineering Institute, University of New with identical interval statistics. The number of active motor units n (t) contributing to the EMG is a function of time, Brunswick, Fredericton, N.B., Canada. I.
,
...
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
418
Neural Pulse Train
W(t)
e(t)
Fig. 3. Block diagram representation of the EMG as an amplitude modulation process. The input, w(t), is a Gaussian white noise signal. It is passed through the linear filter, h(t), which shapes the spectrum of w(t) to the EMG spectrum. Characters with underbars are boldface in text. Fig. 2. Reduced model for EMG signal generation. Each motor unit action potential is assumed to have the shape h(t). Characters with underbars are boldface in text.
where 1xx(T) is the autocorrelation function of the renewal process generating the neural pulse trains. In the Appendix, kxx(T) is shown to have the following form, [eqn. A-lI,
which implies that the output signal, e(t), is a nonstationary stochastic process. If the linear systems are assumed identical, then Fig. 1 0.(-r) = ka [6(r) + 41(r)] reduces to Fig. 2. The nonstationary process e(t) then has with a corresponding power density spectrum (PDS) of a mean value given by (5) 41x(w) = ka I I + ' )xx1. C n(t) The FT of a 2 (t) is then lAe(t)'= E{e(t)=E h (t)* kj E (t - tij) F[c2(t)J E{k}} kaN(Ij.) [H(jNo) * Ii(fjA) x(c)] (6) or (1) where N(jw) = F {Minimum [n (X), n(r)] }. Me(t) = E{kj} kah(t) * n(t) The determination of cDxx(') for different interval statistics where is shown in the Appendix. Fig. A2 shows that Dxx(w) is relatively constant for frequencies above 35 Hz. Since most of is the expectation operation, E{ } the power content of H(co) lies above this frequency (see Fig. * is the convolution operation, 5) the FT for 02(t) can be reduced to 6(t - ha1) is a train of impulses occurring at random times ti', (7) Nho) [H(jo) * H(Uc)J F[ae2(t)] E {k/ } kaN( is the average firing rate of a motor unit train. ka or The Fourier transform (FT) of pj(t) is F[ue2(t)] N(j.X) = (8) (2) F[Me(t)] = E {ki} ka H(IO) N(MA) N1j-)= E kj2} ka4H(/Cj) * H(jc.)1 with H(jco) and N(jc) being the FT of h (t) and n(t) respec- Postulating that [H(jc) * H(jW)] is constant over the fretively. From equation (2) it is clear that the mean value of quency range where F[a,(t)] is appreciable then e(t) is not necessarily zero even though H(jc) represents a F[crf (t)J bandpass system with a zero at the origin. The mean is zero (9) N(jcr ) E k/ 2(0) if either E {k1} = 0 or H(jw) and N(jw) fall in different frequency bands. In practice, the second alternative would be or more likely, and thus, in subsequent derivations, it is assumed n (t) = E {e2 (t)}/K (10) that H(jw) and N(j) are mutually disjoint. Experimental evidence is presented later to support this assumption. where The mean squared value of e(t) is given by K-El{ki} kaH2 (0). oon' rrn(t) Therefore, EMG generation in terms of its mean and variance 2~(t)=E{e2(t)}=Ejh(t)e £k1 E 6(t -t) may be modeled, from eqn. (10), as e(t) = [Kn(t)] 1/2 w(t). b j~~l i=-c The stochastic process, w(t), is a stationary process with zero 1= n(t) and unit variance. Its spectral characteristics are those mean 5 .h(t) *E k, 8(t - tm) (3) EMG during a constant contraction. If w(t) is a Gaussof the m=-C* /=I ian process then only the first and second order moments are Therefore, under the assumption that individual motor units necessary to specify it. A possible block diagram representaare uncorrelated, the above reduces to tion is shown in Fig. 3.
ae(t) = ff dX dT h (t [n ()), n.( r)] f~(7
III. EXPERIMENTS The block diagram representation of EMG generation shown in Fig. 3 was developed under a variety of assumptions. For (4)
--) t - X) E{k, }Minimum
419
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
ch. 1 ch. 2
Fig. 4. Placement of the twisted bipolar electrodes for the motor unit correlation study. The spacing, d, was approximately 12 mm. Each channel was connected to one set of bipolar electrodes with the reference electrode placed on the forearm.
the sake of completeness these assumptions are listed below: (i) Individual motor unit trains are uncorrelated, (ii) k~Aco) is constant over the frequency range where H(ico) has appreciable power, (iii) H(jco) and N(jc) are disjoint, (iv) [H(j&.) * H(jco)] is constant over the frequency range where F[au(t)I is appreciable. The object of the experiments described in this Section is to assess the validity of these assumptions. A. Behavior of Motor Unit Trains Motor unit trains in the model were assumed to be renewal processes and uncorrelated. That individual motor unit trains may be considered to be renewal processes was shown by Clamman(1 2). To show that the individual motor units are uncorrelated, the following experiment was performed. Two sets of twisted bipolar electrodes were inserted intramuscularly("') into the long head of the biceps brachii muscle as shown in Fig. 4. The subject then held weights of 2.3 kgm and 6.8 kgn respectively. The autocorrelation of each EMG was found along with the crosscorrelation between the two EMG's. The electrodes sampled only the activity of motor units in their immediate vicinity. Since the territory occupied by one motor unit is small("*) the two EMG's represent essentially two different groups of motor units. The correlation coefficient between the two groups was then found from
Maximum 012 (r)I 12 |10 1I(0) 02 2 (0) 1/ Values of p obtained are shown in Table 1. From this it can be seen that the crosscorrelation terms contribute less than 10%, justifying their omission in the derivation of equation (4).
TABLE I CORRELATION COEFFICIENT BETWEEN Two SETS OF TWISTED BIPOLAR ELECTRODES 1
EXPERIMENT 2
3
2.3 Kgm
0.030
0.033
0.035
6.8 Kgm
0.105
0.049
0.055
WEIGHT
2.5 cm and (ii) an EMG recorded by surface electrodes during a moderate contraction, are shown in Fig. 5 and Fig. 6, respectively. Comparing these spectra with x(c) (Fig. A2) over the range of frequencies where they have appreciable power, it is evident that 4xx(w) is relatively constant. Thus the PDS of e(t) during an isometric contraction may be approximated by
sDee(@j) -NE[kJ
IH(iW) 12 ka.
I
(12)
Fig. A2 indicates also that the approximation used to derive eqn. (12) is valid regardless of the assumed interval density for the neural pulse trains. Therefore, equation (12) may be interpreted to say that as far as the EMG power density spectrum is concerned, the individual neural pulse trains may be considered to be Poisson point processes with parameter X = k. Equation (12) may be interpreted in a different manner as well. If the individual motor unit shapes were different, that is, if the hi(t) in Fig. 1 were not identical for all i, one readily shows that the PDS of the EMG under an isometric contraction is
4ee((.))=NE{kf}ka [I/N.. IH.(io)12J
(13)
Therefore, H(ijc) 12 in eqn. (1.2) may be considered to be the average of the different IHI(iw)12 or equivalently h(t) represents the average shape of a motor unit pulse. The average shape, h (t), may be obtained by fitting suitable analytical forms to the measured spectra shown in Figs. 5 and 6. A specific form that provides a reasonable fit to the spectra is (14) H(jo) = Kjo/( jw + a)'. The corresponding shape of an average motor unit pulse is (15) h(t) = Kte-at [1 - at/2]. Fig. 7 compares a measured motor unit pulse with the average motor unit pulse obtained from the above form.
C. Experimental Measurement of a2 (t) Two experiments were performed to measure c2 (t). In the first experiment the subject was instructed to isometrically contract, as quickly as possible, the biceps muscle to a maximum effort and then relax immediately. This contraction was repeated twenty to thirty times at a rate which was comfortable for the subject. The resultant EMG activity was a series (11) oee(e) = NE{k/} I H(jco)12 Fx(c) of bursts as shown in Fig. 8. In the second experiment the subject was instructed to proPower density spectra for, (i) a single motor unit recorded by a set of bipolar intramuscular electrodes at a spacing of duce a certain torque about the elbow joint which was mea-
B. Nature of H(jcwJ) To compare H(jco) and 1(c), the PDS of the EMG during a constant contraction was measured. Since the number of active motor units is a constant, N, the EMG is a stationary process with a PDS given by
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
420
20-/1 ~
\s
10 cli -r
\
7
6 5 4
N.
3 K11
l
118
59
176
236
295
f (h z)
354
413
Fig. 5. PDS of a single motor unit action potential recorded from the biceps brachii muscle with surface electrodes. The PDS was found by first calculating H(jo), the amplitude density spectrum and then computing H(Jo) 12. 300 200 100
-
3-
3 ._
(a)
30-
20 10
h(t) = Kte-at _
a = 800
le Approximation I
321 l
_
(2-atj
1 Division = 1 ms
(b) , 4
10
20
100
200
400
Fig. 7. (a) A typical motor unit pulse measured from the biceps brachii muscle. (b) Motor unit pulse obtained from the transfer function Kjw/(/w + a)3. Parameter, a, was found from the Bode plot in Fig. 6 and corresponds to a break frequency of 127 Hz.
Fig. 6. PDS of the EMG for a moderate contraction of the biceps brachii muscle. The EMG was recorded with surface bipolar electrodes placed parallel to the muscle fibers at a spacing of 8 cm. The PDS was calculated by a FFT algorithm using a 10 s record length.
sured by a tensiometer connected to the wrist. The subject would start from a relaxed state, contract isometrically to the desired level at a rate normal to him, maintain this contraction for three to four seconds and then relax. As in the first experiment this was repeated twenty to thirty times. Fig. 9 shows the experimental set-up along with a typical signal. For both experiments the EMG activity of the long head of the biceps brachii was detected by surface electrodes and recorded on a PI 207 FM analog tape recorder at a tape speed of 60 in/s. Before determining a,(t), the recorded signal was played back at a speed of 7 2 in/s onto a chart recorder. The
Fig. 8. Typical burst of EMG activity recorded from the biceps brachii muscle during repetitive isometric contractions. Horizontal scale: 1 div. = 125 ms.
length of each burst was determined and any burst that obviously was not typical either in terms of amplitude or duration was rejected for further processing. The bursts of EMG activity recorded for a particular subject were each considered to belong to the ensemble of functions representing the stochastic process which generated the EMG. The mean squared value, u2(t), was then found by averaging
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
421
(a)
,
A1119 _'W -1il.
2i T1^,1111sk,
L
*--
(b)
Fig. 9. (a) Experiment set-up for the generation of an ensemble of EMG signal. The subject increased muscle force until the tensiometer reading was 6.8 kgm. Horizontal scale: 1 div. = 500 ins. (b) Typical functions of the ensemble recorded from the long head of the biceps brachii muscle.
Tape Recorder Fig.10.Ino
.
e setupaor hve
I
: :
tss
of he
onyionerpole with atimeconsttFolter 6 _
Lug
Level
HP Corredtor
Signal ecovery M~~~~~~~~ode
||
T20 Analog_
Computer|Ge
rao
Fig. 10. Instrumentation set-up for the estimation of the time varying mean squared value, eer The low pass filter had only one pole with a time constant of 6.25 ins.,
across the ensemble as shown in Fig. 10. The recorded signal was squared and low pass filtered on a TR20 Analog Computer. The low pass filter was used to decrease the variance in the estimated °,(t). Its time constant, however, was chosen so that significant bias was not introduced in the estimate. The squarer output was fed to a level detector which signaled the start of a particular member of the ensemble and triggered a sync pulse generator. Ensemble averaging was then performed on the HP 3721A correlator which was set to the sig-
nal recovery mode. In this mode, every time a sync pulse was generated the correlator would digitize a specified number of samples from the signal channel and add them to a running average of the previous digitized sample blocks. The mean squared value, a',(t), found for the two experiments is shown in Fig. 11 and Fig. 12. The amplitude density spectrum of GQ(t) for experiment 1 is shown in Fig. 13 along with the PDS of [H(Ijc) * H(jco)] for an EMG signal. As can be seen [H(jco) * H(jw)] is relatively
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
422
5
-
4
-
2-
3
-
1-
2 -
4-
Subject: TAR M= 26
3-
2 50
125
1
875 1000
750
625
500
375
a) 4-
-a
Subject: ES M 24
3-
M=15
-
a)
5-
-E
4.
225 0
D
3 To
37
, 125
250
375
,
,
,
, 500
750
625
, 875
, 1000
Time (msec)-.
L,-
~~~SbUDect:. ANN r'. .1,
v
1250
-
a)
ar)
1000
750
500
_0
=
2-
Subject: AN
A K
ar
S.ubject:
2-
TAR
M 14 =
1
-
Time (msec )-m
c 0a a) a)
750
500
250
0
1000
1250
5
4
3
3
2-
Subject ES M 19
2
=
12 5
2 50
375
500
62 5
750
875 1000
Fig. 11. Calculated o2(t) for three subjects. Typical signals used to calculate a2(t) are shown in Fig. 8. M is the number of bursts used to calculate a (t).
IT
0
I
2 50
II
500
I1
750
1
1 000
I
12 50
Fig. 12. Calculated a2(t) for the typical signals shown in Fig. 9 (b). M again is the number of ensemble functions used to calculate ae (t).
over the frequency range that a,(t) has appreciable amplitude density. Since the subject was instructed to contract the muscle to a maximum effort as quickly as possible, u, (t) found in this experiment should have the maximum possible bandwidth it can assume. The approximation in eqn. (9) thus is reasonable and the amplitude modulated model for the EMG is valid.
constant
and a standard error of 10% in the estimated a2(t) is desired, then an EMG record length of 200 ms(I 7) is necessary. Use of this record length to estimate ue (t) would obscure completely its transient portion. Design of EMG processors because of the nonstationarity has proven to be difficult. In experimental situations such as gait studies where an ensemble of functions can be generated the method of ensemble averaging may be employed. Another IV. DISCUSSION AND CONCLUSIONS method that is particularly applicable to nonstationary signals An amplitude modulation model for the EMG has been de- is that of Kalman filtering. The application of this important veloped theoretically along with experimental verification of filter for EMG processing has been limited"'8). We hope that the necessary assumptions. The modulating signal was shown the model developed in this paper will facilitate and stimulate to be a function of n (t), the number of active motor units. research in this area. Since muscle force production may be modeled similarly(1 5) then the basic parameter to be estimated by an EMG processor APPENDIX is n (t), or equivalently a, (t). PDS OF A RENEWAL POINT PROCESS In estimating cQ(t), experiments 1 and 2 suggest that the Consider a series of pulses of amplitude A according in time nonstationarity of the EMG must be taken into account. The t such that AAt -+'I as At 0, i.e., a unit impulse. The autoresults of experiment 2 indicate that ce2(t) has a 10%-90% correlation function of this series of pulses for > 0 is: risetime of approximately 300 ms which implies a bandwidth of 1.4 Hz(16) for a2 (t). This bandwidth compares with that ¢(T)= Xl X2 P{X1, X2: T} found for qe2(t) from the first experiment (Fig. 13) and is of where xl is an event (i.e., pulse of amplitude A) occurring at the same order of magnitude as the muscle force bandwidth. time t; x2 is an event occurring at time t + r; P{x1, x2: r} is The difficulty in estimating ae (t) and hence muscle force is the probability of a pulse occurring at t and t + r. illustrated by considering the EMG to be quasi-stationary. If an equivalent bandwidth' of 500 Hz is assumed for the EMG 0(T) =A2P{x } P{x2/x1 ma}. r
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
423
TABLE 2 LAPLACE TRANSFORMS OF FOUR INTERSPIKE INTERVAL DENSITIES Equation
Dens i ty
100I-
\
PDS of H (jw)
*
Graph
or
Ie-ds
Rectangular
H (jw)
ds
d
(4/ s ) [e -/21 ]
Triangular 2 - CT A T6(
Gammna
10765-
Beta
43-
AT6(1
I
100
50
0
200
150
Consider now the following: (
P{event occurring in (r,
hiM
ArO
30[1-6 +.)
417
A Nonstationary Model for the Electromyogram EDWARD -SHWEDYK, MEMBER, IEEE, R. BALASUBRAMANIAN, AND R. N. SCOTT, SENIOR MEMBER, IEEE
Abstract-A theoretical model of the electromyographic (EMG) signal has been developed. In the model, the neural pulse train inputs were considered to be point processes which passed through linear, timeinvariant systems that represented the respective motor unit action potential. The outputs were then summed to produce the EMG. It was assumed, that in the production of muscle force, the controlled parameter was the number of active motor units, n (t). The model then showed that the EMG can be represented as an amplitude modulation process of the form EMG = [Kn (t)J 1/2 w(t) with the stochastic process, w(t), having the spectral and probability characteristics of the EMG during a constant contraction. Various assumptions made in the model development have been verified by experiments.
\Neural Pulse Train
Ae(t) n
Fig.
1.
n
Model for EMG signal generation. The number of active motor
units are n (t). The different motor unit action potentials are repreINTRODUCTION sented by the impulse responses h1 (t),-- hn(t). Coefficients ki represent the respective weighting constants. Characters with underCONTROL of skeletal muscle force is achieved by neural bars are boldface in text. signals, transmitted via motorneurons to the muscle, from the anterior horns of the spinal cord. The motorneuron, along with the muscle fibers it supplies, is called a motor unit and is do not achieve this ideal and design of them has been based, in the basic muscle unit under voluntary control. The muscle part, on intuition and convenience rather than on a fundamenforce produced is controlled primarily by (i) varying the num- tal knowledge of the EMG and its characteristics. This was to ber of activated motor units and (ii) varying the rate of activa- a large extent due to the lack of a suitable model for EMG gention of an individual motor unit. Though both factors are eration. Properties of a particular class of EMG processors present, (i) is generally considered to be the predominant have recently been described(9). Theoretical analysis of the one.(0) processors was based on an amplitude modulation model Muscle force generation is also accompanied by electrical of EMG signal generation. The resemblance of EMG to signal activity which, detected by surface electrodes, produces an amplitude modulation process has been pointed out the electromyogram (EMG). The EMG has a complex tempo- previously(1O0 1 1) ral pattern and is best described as a random signal which may This paper describes the development of an analytical model be modeled as a weighted sum of the action potentials of the of EMG generation. The development shows the restrictions active -motor units. Therefore, its statistical properties are de- and assumptions under which the EMG may be considered to pendent on the number of active motor units and their rate of be an amplitude modulation process. In the model, the moduactivation. lation signal is the number of active motor units while the carA functional relationship thus exists between the EMG and rier is a random signal whose statistical properties are those of muscle force as has been shown in several studies(2,3 ,4). This the EMG during an isometric contraction. Experimental evirelationship has been exploited to use the EMG for diag- dence is presented to justify the assumptions made in the nosis(5), analysis of skill(6) and gait studies"7. Further, since model development. The model and experiments bring into muscle force can be voluntarily controlled and hence so can focus the EMG nonstationarity and the difficulties in processthe EMG, the EMG has been used to control powered orthotic/ ing it. The paper concludes by suggesting a possible further diprosthetic devices(8). rection for research into EMG processors. In the above applications, processors are required which ideally transform the EMG into a noiseless output signal monoII. DEVELOPMENT OF THE MODEL tonically related to muscle contraction or equivalently to the Fig. 1 depicts the generation of EMG in a block diagram number of active motor units. Practical processors of course form. The neural pulse trains go through different linear, time-invariant systems, characterized by impulse responses Manuscript received July 17, 1975; revised July 19, 1976. This work was performed at the Bio-Engineering Institute, University of New h1 (t) h.(t), to generate the motor unit trains which are Brunswick, Fredericton, N.B., Canada, and was supported in part summed to produce the EMG. The weighting coefficients, by grants from Health and Welfare, Canada, and the National Research ki, are assumed to be random and a function of electrode Council of Canada. E. Shwedyk is with the Department of Electrical Engineering, Uni- placement. It is also assumed that the motor unit trains are versity of Manitoba, Winnipeg, Man., Canada. uncorrelated and, that each train is a renewal point process R. Balasubramanian and R. N. Scott are with the Department of Electrical Engineering and the Bio-Engineering Institute, University of New with identical interval statistics. The number of active motor units n (t) contributing to the EMG is a function of time, Brunswick, Fredericton, N.B., Canada. I.
,
...
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
418
Neural Pulse Train
W(t)
e(t)
Fig. 3. Block diagram representation of the EMG as an amplitude modulation process. The input, w(t), is a Gaussian white noise signal. It is passed through the linear filter, h(t), which shapes the spectrum of w(t) to the EMG spectrum. Characters with underbars are boldface in text. Fig. 2. Reduced model for EMG signal generation. Each motor unit action potential is assumed to have the shape h(t). Characters with underbars are boldface in text.
where 1xx(T) is the autocorrelation function of the renewal process generating the neural pulse trains. In the Appendix, kxx(T) is shown to have the following form, [eqn. A-lI,
which implies that the output signal, e(t), is a nonstationary stochastic process. If the linear systems are assumed identical, then Fig. 1 0.(-r) = ka [6(r) + 41(r)] reduces to Fig. 2. The nonstationary process e(t) then has with a corresponding power density spectrum (PDS) of a mean value given by (5) 41x(w) = ka I I + ' )xx1. C n(t) The FT of a 2 (t) is then lAe(t)'= E{e(t)=E h (t)* kj E (t - tij) F[c2(t)J E{k}} kaN(Ij.) [H(jNo) * Ii(fjA) x(c)] (6) or (1) where N(jw) = F {Minimum [n (X), n(r)] }. Me(t) = E{kj} kah(t) * n(t) The determination of cDxx(') for different interval statistics where is shown in the Appendix. Fig. A2 shows that Dxx(w) is relatively constant for frequencies above 35 Hz. Since most of is the expectation operation, E{ } the power content of H(co) lies above this frequency (see Fig. * is the convolution operation, 5) the FT for 02(t) can be reduced to 6(t - ha1) is a train of impulses occurring at random times ti', (7) Nho) [H(jo) * H(Uc)J F[ae2(t)] E {k/ } kaN( is the average firing rate of a motor unit train. ka or The Fourier transform (FT) of pj(t) is F[ue2(t)] N(j.X) = (8) (2) F[Me(t)] = E {ki} ka H(IO) N(MA) N1j-)= E kj2} ka4H(/Cj) * H(jc.)1 with H(jco) and N(jc) being the FT of h (t) and n(t) respec- Postulating that [H(jc) * H(jW)] is constant over the fretively. From equation (2) it is clear that the mean value of quency range where F[a,(t)] is appreciable then e(t) is not necessarily zero even though H(jc) represents a F[crf (t)J bandpass system with a zero at the origin. The mean is zero (9) N(jcr ) E k/ 2(0) if either E {k1} = 0 or H(jw) and N(jw) fall in different frequency bands. In practice, the second alternative would be or more likely, and thus, in subsequent derivations, it is assumed n (t) = E {e2 (t)}/K (10) that H(jw) and N(j) are mutually disjoint. Experimental evidence is presented later to support this assumption. where The mean squared value of e(t) is given by K-El{ki} kaH2 (0). oon' rrn(t) Therefore, EMG generation in terms of its mean and variance 2~(t)=E{e2(t)}=Ejh(t)e £k1 E 6(t -t) may be modeled, from eqn. (10), as e(t) = [Kn(t)] 1/2 w(t). b j~~l i=-c The stochastic process, w(t), is a stationary process with zero 1= n(t) and unit variance. Its spectral characteristics are those mean 5 .h(t) *E k, 8(t - tm) (3) EMG during a constant contraction. If w(t) is a Gaussof the m=-C* /=I ian process then only the first and second order moments are Therefore, under the assumption that individual motor units necessary to specify it. A possible block diagram representaare uncorrelated, the above reduces to tion is shown in Fig. 3.
ae(t) = ff dX dT h (t [n ()), n.( r)] f~(7
III. EXPERIMENTS The block diagram representation of EMG generation shown in Fig. 3 was developed under a variety of assumptions. For (4)
--) t - X) E{k, }Minimum
419
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
ch. 1 ch. 2
Fig. 4. Placement of the twisted bipolar electrodes for the motor unit correlation study. The spacing, d, was approximately 12 mm. Each channel was connected to one set of bipolar electrodes with the reference electrode placed on the forearm.
the sake of completeness these assumptions are listed below: (i) Individual motor unit trains are uncorrelated, (ii) k~Aco) is constant over the frequency range where H(ico) has appreciable power, (iii) H(jco) and N(jc) are disjoint, (iv) [H(j&.) * H(jco)] is constant over the frequency range where F[au(t)I is appreciable. The object of the experiments described in this Section is to assess the validity of these assumptions. A. Behavior of Motor Unit Trains Motor unit trains in the model were assumed to be renewal processes and uncorrelated. That individual motor unit trains may be considered to be renewal processes was shown by Clamman(1 2). To show that the individual motor units are uncorrelated, the following experiment was performed. Two sets of twisted bipolar electrodes were inserted intramuscularly("') into the long head of the biceps brachii muscle as shown in Fig. 4. The subject then held weights of 2.3 kgm and 6.8 kgn respectively. The autocorrelation of each EMG was found along with the crosscorrelation between the two EMG's. The electrodes sampled only the activity of motor units in their immediate vicinity. Since the territory occupied by one motor unit is small("*) the two EMG's represent essentially two different groups of motor units. The correlation coefficient between the two groups was then found from
Maximum 012 (r)I 12 |10 1I(0) 02 2 (0) 1/ Values of p obtained are shown in Table 1. From this it can be seen that the crosscorrelation terms contribute less than 10%, justifying their omission in the derivation of equation (4).
TABLE I CORRELATION COEFFICIENT BETWEEN Two SETS OF TWISTED BIPOLAR ELECTRODES 1
EXPERIMENT 2
3
2.3 Kgm
0.030
0.033
0.035
6.8 Kgm
0.105
0.049
0.055
WEIGHT
2.5 cm and (ii) an EMG recorded by surface electrodes during a moderate contraction, are shown in Fig. 5 and Fig. 6, respectively. Comparing these spectra with x(c) (Fig. A2) over the range of frequencies where they have appreciable power, it is evident that 4xx(w) is relatively constant. Thus the PDS of e(t) during an isometric contraction may be approximated by
sDee(@j) -NE[kJ
IH(iW) 12 ka.
I
(12)
Fig. A2 indicates also that the approximation used to derive eqn. (12) is valid regardless of the assumed interval density for the neural pulse trains. Therefore, equation (12) may be interpreted to say that as far as the EMG power density spectrum is concerned, the individual neural pulse trains may be considered to be Poisson point processes with parameter X = k. Equation (12) may be interpreted in a different manner as well. If the individual motor unit shapes were different, that is, if the hi(t) in Fig. 1 were not identical for all i, one readily shows that the PDS of the EMG under an isometric contraction is
4ee((.))=NE{kf}ka [I/N.. IH.(io)12J
(13)
Therefore, H(ijc) 12 in eqn. (1.2) may be considered to be the average of the different IHI(iw)12 or equivalently h(t) represents the average shape of a motor unit pulse. The average shape, h (t), may be obtained by fitting suitable analytical forms to the measured spectra shown in Figs. 5 and 6. A specific form that provides a reasonable fit to the spectra is (14) H(jo) = Kjo/( jw + a)'. The corresponding shape of an average motor unit pulse is (15) h(t) = Kte-at [1 - at/2]. Fig. 7 compares a measured motor unit pulse with the average motor unit pulse obtained from the above form.
C. Experimental Measurement of a2 (t) Two experiments were performed to measure c2 (t). In the first experiment the subject was instructed to isometrically contract, as quickly as possible, the biceps muscle to a maximum effort and then relax immediately. This contraction was repeated twenty to thirty times at a rate which was comfortable for the subject. The resultant EMG activity was a series (11) oee(e) = NE{k/} I H(jco)12 Fx(c) of bursts as shown in Fig. 8. In the second experiment the subject was instructed to proPower density spectra for, (i) a single motor unit recorded by a set of bipolar intramuscular electrodes at a spacing of duce a certain torque about the elbow joint which was mea-
B. Nature of H(jcwJ) To compare H(jco) and 1(c), the PDS of the EMG during a constant contraction was measured. Since the number of active motor units is a constant, N, the EMG is a stationary process with a PDS given by
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
420
20-/1 ~
\s
10 cli -r
\
7
6 5 4
N.
3 K11
l
118
59
176
236
295
f (h z)
354
413
Fig. 5. PDS of a single motor unit action potential recorded from the biceps brachii muscle with surface electrodes. The PDS was found by first calculating H(jo), the amplitude density spectrum and then computing H(Jo) 12. 300 200 100
-
3-
3 ._
(a)
30-
20 10
h(t) = Kte-at _
a = 800
le Approximation I
321 l
_
(2-atj
1 Division = 1 ms
(b) , 4
10
20
100
200
400
Fig. 7. (a) A typical motor unit pulse measured from the biceps brachii muscle. (b) Motor unit pulse obtained from the transfer function Kjw/(/w + a)3. Parameter, a, was found from the Bode plot in Fig. 6 and corresponds to a break frequency of 127 Hz.
Fig. 6. PDS of the EMG for a moderate contraction of the biceps brachii muscle. The EMG was recorded with surface bipolar electrodes placed parallel to the muscle fibers at a spacing of 8 cm. The PDS was calculated by a FFT algorithm using a 10 s record length.
sured by a tensiometer connected to the wrist. The subject would start from a relaxed state, contract isometrically to the desired level at a rate normal to him, maintain this contraction for three to four seconds and then relax. As in the first experiment this was repeated twenty to thirty times. Fig. 9 shows the experimental set-up along with a typical signal. For both experiments the EMG activity of the long head of the biceps brachii was detected by surface electrodes and recorded on a PI 207 FM analog tape recorder at a tape speed of 60 in/s. Before determining a,(t), the recorded signal was played back at a speed of 7 2 in/s onto a chart recorder. The
Fig. 8. Typical burst of EMG activity recorded from the biceps brachii muscle during repetitive isometric contractions. Horizontal scale: 1 div. = 125 ms.
length of each burst was determined and any burst that obviously was not typical either in terms of amplitude or duration was rejected for further processing. The bursts of EMG activity recorded for a particular subject were each considered to belong to the ensemble of functions representing the stochastic process which generated the EMG. The mean squared value, u2(t), was then found by averaging
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
421
(a)
,
A1119 _'W -1il.
2i T1^,1111sk,
L
*--
(b)
Fig. 9. (a) Experiment set-up for the generation of an ensemble of EMG signal. The subject increased muscle force until the tensiometer reading was 6.8 kgm. Horizontal scale: 1 div. = 500 ins. (b) Typical functions of the ensemble recorded from the long head of the biceps brachii muscle.
Tape Recorder Fig.10.Ino
.
e setupaor hve
I
: :
tss
of he
onyionerpole with atimeconsttFolter 6 _
Lug
Level
HP Corredtor
Signal ecovery M~~~~~~~~ode
||
T20 Analog_
Computer|Ge
rao
Fig. 10. Instrumentation set-up for the estimation of the time varying mean squared value, eer The low pass filter had only one pole with a time constant of 6.25 ins.,
across the ensemble as shown in Fig. 10. The recorded signal was squared and low pass filtered on a TR20 Analog Computer. The low pass filter was used to decrease the variance in the estimated °,(t). Its time constant, however, was chosen so that significant bias was not introduced in the estimate. The squarer output was fed to a level detector which signaled the start of a particular member of the ensemble and triggered a sync pulse generator. Ensemble averaging was then performed on the HP 3721A correlator which was set to the sig-
nal recovery mode. In this mode, every time a sync pulse was generated the correlator would digitize a specified number of samples from the signal channel and add them to a running average of the previous digitized sample blocks. The mean squared value, a',(t), found for the two experiments is shown in Fig. 11 and Fig. 12. The amplitude density spectrum of GQ(t) for experiment 1 is shown in Fig. 13 along with the PDS of [H(Ijc) * H(jco)] for an EMG signal. As can be seen [H(jco) * H(jw)] is relatively
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 5, SEPTEMBER 1977
422
5
-
4
-
2-
3
-
1-
2 -
4-
Subject: TAR M= 26
3-
2 50
125
1
875 1000
750
625
500
375
a) 4-
-a
Subject: ES M 24
3-
M=15
-
a)
5-
-E
4.
225 0
D
3 To
37
, 125
250
375
,
,
,
, 500
750
625
, 875
, 1000
Time (msec)-.
L,-
~~~SbUDect:. ANN r'. .1,
v
1250
-
a)
ar)
1000
750
500
_0
=
2-
Subject: AN
A K
ar
S.ubject:
2-
TAR
M 14 =
1
-
Time (msec )-m
c 0a a) a)
750
500
250
0
1000
1250
5
4
3
3
2-
Subject ES M 19
2
=
12 5
2 50
375
500
62 5
750
875 1000
Fig. 11. Calculated o2(t) for three subjects. Typical signals used to calculate a2(t) are shown in Fig. 8. M is the number of bursts used to calculate a (t).
IT
0
I
2 50
II
500
I1
750
1
1 000
I
12 50
Fig. 12. Calculated a2(t) for the typical signals shown in Fig. 9 (b). M again is the number of ensemble functions used to calculate ae (t).
over the frequency range that a,(t) has appreciable amplitude density. Since the subject was instructed to contract the muscle to a maximum effort as quickly as possible, u, (t) found in this experiment should have the maximum possible bandwidth it can assume. The approximation in eqn. (9) thus is reasonable and the amplitude modulated model for the EMG is valid.
constant
and a standard error of 10% in the estimated a2(t) is desired, then an EMG record length of 200 ms(I 7) is necessary. Use of this record length to estimate ue (t) would obscure completely its transient portion. Design of EMG processors because of the nonstationarity has proven to be difficult. In experimental situations such as gait studies where an ensemble of functions can be generated the method of ensemble averaging may be employed. Another IV. DISCUSSION AND CONCLUSIONS method that is particularly applicable to nonstationary signals An amplitude modulation model for the EMG has been de- is that of Kalman filtering. The application of this important veloped theoretically along with experimental verification of filter for EMG processing has been limited"'8). We hope that the necessary assumptions. The modulating signal was shown the model developed in this paper will facilitate and stimulate to be a function of n (t), the number of active motor units. research in this area. Since muscle force production may be modeled similarly(1 5) then the basic parameter to be estimated by an EMG processor APPENDIX is n (t), or equivalently a, (t). PDS OF A RENEWAL POINT PROCESS In estimating cQ(t), experiments 1 and 2 suggest that the Consider a series of pulses of amplitude A according in time nonstationarity of the EMG must be taken into account. The t such that AAt -+'I as At 0, i.e., a unit impulse. The autoresults of experiment 2 indicate that ce2(t) has a 10%-90% correlation function of this series of pulses for > 0 is: risetime of approximately 300 ms which implies a bandwidth of 1.4 Hz(16) for a2 (t). This bandwidth compares with that ¢(T)= Xl X2 P{X1, X2: T} found for qe2(t) from the first experiment (Fig. 13) and is of where xl is an event (i.e., pulse of amplitude A) occurring at the same order of magnitude as the muscle force bandwidth. time t; x2 is an event occurring at time t + r; P{x1, x2: r} is The difficulty in estimating ae (t) and hence muscle force is the probability of a pulse occurring at t and t + r. illustrated by considering the EMG to be quasi-stationary. If an equivalent bandwidth' of 500 Hz is assumed for the EMG 0(T) =A2P{x } P{x2/x1 ma}. r
SHWEDYK et al.: MODEL FOR THE ELECTROMYOGRAM
423
TABLE 2 LAPLACE TRANSFORMS OF FOUR INTERSPIKE INTERVAL DENSITIES Equation
Dens i ty
100I-
\
PDS of H (jw)
*
Graph
or
Ie-ds
Rectangular
H (jw)
ds
d
(4/ s ) [e -/21 ]
Triangular 2 - CT A T6(
Gammna
10765-
Beta
43-
AT6(1
I
100
50
0
200
150
Consider now the following: (
P{event occurring in (r,
hiM
ArO
30[1-6 +.)
