487
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-22, NO. 6, NOVEMBER 1975 Select (-)
Sampled ECG
N-point
transform',
transform
com-
using the variance criterion
ponents
A
Storage edium
and human ECG's have similar characteristics, it is plausible that the method presented in this paper could also be used to study data compression of human ECG's. ACKNOWLEDGMENT The authors wish to express their indebtedness to Dr. W. W.
(a)
Koepsel and Mr. T. Natarajan of the Department of Electrical Engineering, Kansas State University, for their assistance in processing the ECG data used in the preparation of this paper. REFERENCES [11 H. C. Andrews, Introduction to Mathematical Techniques in Pat-
(b) Fig. 4. Pertaining to an m :1 data compression.
(b) Retrieval.
(a) Storage.
Using orthogonal transforms, we store N/m transform components per ECG. Each component is stored as a word in storage. Thus the total number of words required is NK/m, which implies that the storage requirements are reduced by a factor m. VI. CONCLUSIONS The experimental results presented in the last section demonstrate that it is feasible to use the variance criterion to secure data compression of canine ECG data. Since canine
tern Recognition. New York: John Wiley & Sons, 1972. [2] H. C. Andrews, "Multidimensional rotations in feature selection," IEEE Trans. on Computers, vol. C-20, pp. 1045-1051, Sept. 1971. [31 N. Ahmed and K. R. Rao, Orthogonal Transforms for Digital Signal Processing. New York/Berlin/Heidelberg: Springer Verlag, (in press). [41 J. Pearl, "Basis restricted transformations and performance measures for spectral representation," IEEE Trans. Info. Theory, vol. IT-17, pp. 751-752, 1971. [51 P. J. Milne, "Orthogonal transform processing of electrocardiographic data," Ph.D. Dissertation, 1973, Kansas State University, Manhattan, Kansas. [61 N. Ahmed et al., "Discrete cosine transform," IEEE Trans. on Computers, vol. C-23, pp. 90-93, Jan. 1974. [71 C. A. Careres and L. S. Dreifus, Clinical Electrocardiography and Computers. New York: Academic Press, 1970. [81 R. C. Balda, "Computer assisted ECG interpretation," Measuring for Medicine, vol. 7, May-Aug. 1972, Hewlett-Packard, Waltham, Mass.
Statistically Constrained Inverse Electrocardiography RICHARD
0.
MARTIN, MEMBER, IEEE, T. C. PILKINGTON, MEMBER, IEEE,
Abstract-This paper examines the feasbility of utilizing statistical constraints on the inverse potential model to determine the potential distribution over a 4 cm sphere surrounding the heart from perturbed torso potentials. These perturbed torso potentials reflect instrumentation, quadrature, electrode placement, and heart position uncertainties. This work is an extension of the authors' previous work which concluded that it is not feasible to determine this same potential distribution using unconstrained solutions. However, the results of the present work indicate that with the use of approximate signal and noise covariance matrices, it is possible to achieve estimates of this potential distribution with an average sum squared error of twenty-five percent. Further, the estimation of the signal and noise covariance matrices can be accomplished with a knowledge of heart geometry, torso geometry,
Manuscript received May 26, 1972; revised November 11, 1974, and March 24, 1975. This paper was supported in part by USPHS Grants HL 05716, HL 05372, and HL 11307. R. 0. Martin is with the Department of Electrical Engineering, Christian Brothers College, Memphis, Tenn. 38104. T. C. Pilkington and M. N. Morrow are with the Department of Biomedical Engineering, Duke University, Durham, N.C. 27706.
AND
MARY N. MORROW
the approximate measurement exror, and a rough estimate of the time an average section of myocardium is depolarized, but without an a prori specification of the activation sequence.
I. INTRODUCTION
ATTEMPTS at obtaining a physiologically meaningful solution to the inverse problem in electrocardiography have been numerous. Proposed models include dipoles [1]- [4], multipoles [5]-[9], and epicardial potentials [10]. This paper
is an extension of the epicardial potential model examined in our previous paper [10]. The same field theory developments are utilized to relate the torso potentials to the potential at any point on a 4 cm sphere surrounding the heart. The geometrical heart and torso data and the activation data used in the simulations are the same dog data previously described [11], [121 and utilized in subsequent work [101-[131. This paper extends the previous work by injecting an assumed
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, NOVEMBER 1975
488
knowledge of the signal and noise statistics [14]-[16] into the formulation of the solution vectors. The statistical constraint approach to inverse electrocardiography differs conceptually from other constrained solution approaches [4], [12], [17] in one major way. Rather than assuming deterministic relationships which imply the functional form of the solution vectors, we assume a knowledge of the average, the variance, and the covariance of the solution vectors. II. THEORETICAL DEVELOPMENTS
In the previous paper [10] we showed torso potentials were related to the potentials over a 4 cm sphere surrounding the heart ("epicardial potentials") by the equation
y =Ax +
(1)
can be utilized when m < n, but the overdetermined case of m > n is of sufficient interest to mention the following alternate expression for the B matrix which was developed by Strand and Westwater [15], [16]. B = (S-1 +A TN-'A)-'A TN-1 (8) Equation (8) requires the inversion of an n * n matrix when m > n, the inversion of a typical m * m diagonal N-matrix being trivial. III. MODELS AND DATA
The sets of geometry and torso potential data utilized in this study were obtained from a single dog, affectionately named Coriolis. A complete description of the origin of the Coriolis data is given by Rogers [13].
where y (m dimensional vector) represents the torso potentials, x (n dimensional vector) represents the "epicardial potentials," A. Geometry Data and q (m dimensional vector), the noise. The matrix A was Two geometrical models were used in this study. The first derived using potential theory and is an m * n transfer matrix. model consists of a 49 segment torso and a 32 segment epiUsing an average mean square error criterion, the infeasibility cardial surface while the second model is comprised of a 30 of estimating x from a knowledge of y and A in the presence segment torso and a 72 segment epicardial surface. The torso of realistic t? was demonstrated. data for both models represent subsets of the original 268 The present work again assumes that the basic form of the Coriolis torso segments. In each case the majority of the torso system of equations is given by (1), but the approach utilized segnents are those located closest to the heart with the rein this paper also assumes a knowledge of the ensemble of mainder chosen to provide continuity over the surface. "signal" and "noise" vectors. Thus we define In order to utilize an epicardial surface that would both bear some relationship to the real epicardial surface and also be S = ((x (x)) (x (x))T) (2) mathematically tractable, a 4 cm sphere enclosing the Coriolis and heart surface and contained within the torso boundary was selected. The sphere was segmented into rectilinear spherical surface elements; 32 segments for the first model and 72 segwhere S is an n * n variance-covariance matrix of the source ments for the second model. vector, N is an m * m variance-covariance matrix of the noise vector, and the brackets represent the ensemble average. B. Torso and Epicardial Potentials The following theorem due to Foster [14] forms the basis of The potential data used in this study were computed directly a statistical estimation scheme by defining an optimum from the activation data utilizing the infinite media transfer estimator. matrix. The relationship used to compute these potentials is Theorem: The optimum inverse B of the m * n matrix A is given in (9). (4) 76 r..- f.) B=SAT(ASAT+N)-l if ASiTO(i). = 2.693 (13 VJ(T) (9) while the optimum estimate of the solution vector x is i=l
le=By +b
(5)
b = (x) - BA (x)
V,(T) is the potential at the jth point for the Tth time instant;
rij is the vector from the ith heart segment to the jth observa-
where
(6)
under the following criterion of optimality (B, b) = ((x - X)T G(x - x)) = miniimum (7) where G is any arbitrary metric [17] . The function may be physically interpreted as the mean square distance between the true solution x and its estimate, k, and as such is a reasonable function to minimize. Equation (4) requires the inversion of an m * m matrix which limits the number of observations which can be easily manipulated when m > n. Statistical estimation techniques
tion point; fti is the outward unit normal to the ith heart segment; ASi is the corresponding surface area; 5 TQ) is a step function that is zero when the ith heart segment dipole is inactive and one otherwise. The turn-on and turn-off times for each 6 (i) are deduced from the measured activation data and correspond to depolarization and repolarization of the heart segnent. The index j runs to 49 for the torso and 32 for the spherical epicardial surface for the first model and to 30 for
the torso and 72 for the spherical epicardial surface for the second model. The constant, 2.693, is an empirical weighting factor which was chosen so that the activation data would fit the infinite media torso data [18].
489
MARTIN et al.: INVERSE ELECTROCARDIOGRAPHY
C. Activation Data The activation data (the ST(i) for each of the 34 time instants) are consistent with the aforementioned geometry and are derived from a limited number of electrode measurements recorded from Coriolis. Gaps in the experimental data were filled by educated judgment. These data are grossly consistent with the published results of Sodi-Pallares [20], [21], Scher [22]-[24], and Durrer [25]- [27]. Normal activation data were used to compute compatible sets of epicardial and torso potentials using the infinite medium transfer coefficients. IV. STATISTICAL ESTIMATION The following sections will present the methods used to estimate the statistical matrices.
A. StatisticalMatrices and Error Estimates Integral to the technique of statistical estimation is the construction of two variance-covariance matrices as defined in (2) and (3). The ensemble average denoted by ( ) implies that one should time align many sets of data and then average over these sets for each time instant. Thus, the S and N matrices obtained would be different for each time instant. This is physically reasonable and intuitively appealing since the potential distribution over the 4 cm spherical surface, hereafter called the epicardial surface, is highly time dependent. However, there are not sufficient data presently available to perform this ensemble averaging. Further, the authors believe that while constrained inverse electrocardiographic techniques must restrict the solution space, useful techniques should be simple and not overly restrictive. Therefore, we decided to test some presumptuous methods for estimating these two variance-covariance matrices. B. Methods of Estimating the S Matrix Two methods for estimating an S matrix are presented. The first method, called the "time average estimate" assumes that the ensemble average can be replaced by a time average. Epicardial potentials are generated from the activation data by using (9). After this is done for all 34 time instants, the S matrix is generated by approximating the ensemble average with temporal summation. That is, each Sii term of the S matrix is generated using the relation
Sij=
134
1: (xi (T) (xi>))(xi (T) (xi)) -
T=1
-
(I10)
where 1
34
(xi)=3 E Xi£(T). T=1
(1
Results suggested that ensemble averaging can be accurately approximated by temporal summation, but this is a very limited technique since the generation of the S matrix requires specification of the activation sequence. However, if the "time average estimate" had not produced an S matrix sufficiently accurate to estimate time dependent epicardial poten-
tials, then it would be very doubtful that ensemble averaging could be replaced by any form of temporal averaging. The second method used to estimate an S matrix is called the "Monte-Carlo Method." It is assumed that each of the 76 dipoles has the same preassigned probability P of being on, and the probability that each dipole is off is 1 - P. A particular activation state is obtained by generating an independent random number Ni(O
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-22, NO. 6, NOVEMBER 1975 Select (-)
Sampled ECG
N-point
transform',
transform
com-
using the variance criterion
ponents
A
Storage edium
and human ECG's have similar characteristics, it is plausible that the method presented in this paper could also be used to study data compression of human ECG's. ACKNOWLEDGMENT The authors wish to express their indebtedness to Dr. W. W.
(a)
Koepsel and Mr. T. Natarajan of the Department of Electrical Engineering, Kansas State University, for their assistance in processing the ECG data used in the preparation of this paper. REFERENCES [11 H. C. Andrews, Introduction to Mathematical Techniques in Pat-
(b) Fig. 4. Pertaining to an m :1 data compression.
(b) Retrieval.
(a) Storage.
Using orthogonal transforms, we store N/m transform components per ECG. Each component is stored as a word in storage. Thus the total number of words required is NK/m, which implies that the storage requirements are reduced by a factor m. VI. CONCLUSIONS The experimental results presented in the last section demonstrate that it is feasible to use the variance criterion to secure data compression of canine ECG data. Since canine
tern Recognition. New York: John Wiley & Sons, 1972. [2] H. C. Andrews, "Multidimensional rotations in feature selection," IEEE Trans. on Computers, vol. C-20, pp. 1045-1051, Sept. 1971. [31 N. Ahmed and K. R. Rao, Orthogonal Transforms for Digital Signal Processing. New York/Berlin/Heidelberg: Springer Verlag, (in press). [41 J. Pearl, "Basis restricted transformations and performance measures for spectral representation," IEEE Trans. Info. Theory, vol. IT-17, pp. 751-752, 1971. [51 P. J. Milne, "Orthogonal transform processing of electrocardiographic data," Ph.D. Dissertation, 1973, Kansas State University, Manhattan, Kansas. [61 N. Ahmed et al., "Discrete cosine transform," IEEE Trans. on Computers, vol. C-23, pp. 90-93, Jan. 1974. [71 C. A. Careres and L. S. Dreifus, Clinical Electrocardiography and Computers. New York: Academic Press, 1970. [81 R. C. Balda, "Computer assisted ECG interpretation," Measuring for Medicine, vol. 7, May-Aug. 1972, Hewlett-Packard, Waltham, Mass.
Statistically Constrained Inverse Electrocardiography RICHARD
0.
MARTIN, MEMBER, IEEE, T. C. PILKINGTON, MEMBER, IEEE,
Abstract-This paper examines the feasbility of utilizing statistical constraints on the inverse potential model to determine the potential distribution over a 4 cm sphere surrounding the heart from perturbed torso potentials. These perturbed torso potentials reflect instrumentation, quadrature, electrode placement, and heart position uncertainties. This work is an extension of the authors' previous work which concluded that it is not feasible to determine this same potential distribution using unconstrained solutions. However, the results of the present work indicate that with the use of approximate signal and noise covariance matrices, it is possible to achieve estimates of this potential distribution with an average sum squared error of twenty-five percent. Further, the estimation of the signal and noise covariance matrices can be accomplished with a knowledge of heart geometry, torso geometry,
Manuscript received May 26, 1972; revised November 11, 1974, and March 24, 1975. This paper was supported in part by USPHS Grants HL 05716, HL 05372, and HL 11307. R. 0. Martin is with the Department of Electrical Engineering, Christian Brothers College, Memphis, Tenn. 38104. T. C. Pilkington and M. N. Morrow are with the Department of Biomedical Engineering, Duke University, Durham, N.C. 27706.
AND
MARY N. MORROW
the approximate measurement exror, and a rough estimate of the time an average section of myocardium is depolarized, but without an a prori specification of the activation sequence.
I. INTRODUCTION
ATTEMPTS at obtaining a physiologically meaningful solution to the inverse problem in electrocardiography have been numerous. Proposed models include dipoles [1]- [4], multipoles [5]-[9], and epicardial potentials [10]. This paper
is an extension of the epicardial potential model examined in our previous paper [10]. The same field theory developments are utilized to relate the torso potentials to the potential at any point on a 4 cm sphere surrounding the heart. The geometrical heart and torso data and the activation data used in the simulations are the same dog data previously described [11], [121 and utilized in subsequent work [101-[131. This paper extends the previous work by injecting an assumed
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, NOVEMBER 1975
488
knowledge of the signal and noise statistics [14]-[16] into the formulation of the solution vectors. The statistical constraint approach to inverse electrocardiography differs conceptually from other constrained solution approaches [4], [12], [17] in one major way. Rather than assuming deterministic relationships which imply the functional form of the solution vectors, we assume a knowledge of the average, the variance, and the covariance of the solution vectors. II. THEORETICAL DEVELOPMENTS
In the previous paper [10] we showed torso potentials were related to the potentials over a 4 cm sphere surrounding the heart ("epicardial potentials") by the equation
y =Ax +
(1)
can be utilized when m < n, but the overdetermined case of m > n is of sufficient interest to mention the following alternate expression for the B matrix which was developed by Strand and Westwater [15], [16]. B = (S-1 +A TN-'A)-'A TN-1 (8) Equation (8) requires the inversion of an n * n matrix when m > n, the inversion of a typical m * m diagonal N-matrix being trivial. III. MODELS AND DATA
The sets of geometry and torso potential data utilized in this study were obtained from a single dog, affectionately named Coriolis. A complete description of the origin of the Coriolis data is given by Rogers [13].
where y (m dimensional vector) represents the torso potentials, x (n dimensional vector) represents the "epicardial potentials," A. Geometry Data and q (m dimensional vector), the noise. The matrix A was Two geometrical models were used in this study. The first derived using potential theory and is an m * n transfer matrix. model consists of a 49 segment torso and a 32 segment epiUsing an average mean square error criterion, the infeasibility cardial surface while the second model is comprised of a 30 of estimating x from a knowledge of y and A in the presence segment torso and a 72 segment epicardial surface. The torso of realistic t? was demonstrated. data for both models represent subsets of the original 268 The present work again assumes that the basic form of the Coriolis torso segments. In each case the majority of the torso system of equations is given by (1), but the approach utilized segnents are those located closest to the heart with the rein this paper also assumes a knowledge of the ensemble of mainder chosen to provide continuity over the surface. "signal" and "noise" vectors. Thus we define In order to utilize an epicardial surface that would both bear some relationship to the real epicardial surface and also be S = ((x (x)) (x (x))T) (2) mathematically tractable, a 4 cm sphere enclosing the Coriolis and heart surface and contained within the torso boundary was selected. The sphere was segmented into rectilinear spherical surface elements; 32 segments for the first model and 72 segwhere S is an n * n variance-covariance matrix of the source ments for the second model. vector, N is an m * m variance-covariance matrix of the noise vector, and the brackets represent the ensemble average. B. Torso and Epicardial Potentials The following theorem due to Foster [14] forms the basis of The potential data used in this study were computed directly a statistical estimation scheme by defining an optimum from the activation data utilizing the infinite media transfer estimator. matrix. The relationship used to compute these potentials is Theorem: The optimum inverse B of the m * n matrix A is given in (9). (4) 76 r..- f.) B=SAT(ASAT+N)-l if ASiTO(i). = 2.693 (13 VJ(T) (9) while the optimum estimate of the solution vector x is i=l
le=By +b
(5)
b = (x) - BA (x)
V,(T) is the potential at the jth point for the Tth time instant;
rij is the vector from the ith heart segment to the jth observa-
where
(6)
under the following criterion of optimality (B, b) = ((x - X)T G(x - x)) = miniimum (7) where G is any arbitrary metric [17] . The function may be physically interpreted as the mean square distance between the true solution x and its estimate, k, and as such is a reasonable function to minimize. Equation (4) requires the inversion of an m * m matrix which limits the number of observations which can be easily manipulated when m > n. Statistical estimation techniques
tion point; fti is the outward unit normal to the ith heart segment; ASi is the corresponding surface area; 5 TQ) is a step function that is zero when the ith heart segment dipole is inactive and one otherwise. The turn-on and turn-off times for each 6 (i) are deduced from the measured activation data and correspond to depolarization and repolarization of the heart segnent. The index j runs to 49 for the torso and 32 for the spherical epicardial surface for the first model and to 30 for
the torso and 72 for the spherical epicardial surface for the second model. The constant, 2.693, is an empirical weighting factor which was chosen so that the activation data would fit the infinite media torso data [18].
489
MARTIN et al.: INVERSE ELECTROCARDIOGRAPHY
C. Activation Data The activation data (the ST(i) for each of the 34 time instants) are consistent with the aforementioned geometry and are derived from a limited number of electrode measurements recorded from Coriolis. Gaps in the experimental data were filled by educated judgment. These data are grossly consistent with the published results of Sodi-Pallares [20], [21], Scher [22]-[24], and Durrer [25]- [27]. Normal activation data were used to compute compatible sets of epicardial and torso potentials using the infinite medium transfer coefficients. IV. STATISTICAL ESTIMATION The following sections will present the methods used to estimate the statistical matrices.
A. StatisticalMatrices and Error Estimates Integral to the technique of statistical estimation is the construction of two variance-covariance matrices as defined in (2) and (3). The ensemble average denoted by ( ) implies that one should time align many sets of data and then average over these sets for each time instant. Thus, the S and N matrices obtained would be different for each time instant. This is physically reasonable and intuitively appealing since the potential distribution over the 4 cm spherical surface, hereafter called the epicardial surface, is highly time dependent. However, there are not sufficient data presently available to perform this ensemble averaging. Further, the authors believe that while constrained inverse electrocardiographic techniques must restrict the solution space, useful techniques should be simple and not overly restrictive. Therefore, we decided to test some presumptuous methods for estimating these two variance-covariance matrices. B. Methods of Estimating the S Matrix Two methods for estimating an S matrix are presented. The first method, called the "time average estimate" assumes that the ensemble average can be replaced by a time average. Epicardial potentials are generated from the activation data by using (9). After this is done for all 34 time instants, the S matrix is generated by approximating the ensemble average with temporal summation. That is, each Sii term of the S matrix is generated using the relation
Sij=
134
1: (xi (T) (xi>))(xi (T) (xi)) -
T=1
-
(I10)
where 1
34
(xi)=3 E Xi£(T). T=1
(1
Results suggested that ensemble averaging can be accurately approximated by temporal summation, but this is a very limited technique since the generation of the S matrix requires specification of the activation sequence. However, if the "time average estimate" had not produced an S matrix sufficiently accurate to estimate time dependent epicardial poten-
tials, then it would be very doubtful that ensemble averaging could be replaced by any form of temporal averaging. The second method used to estimate an S matrix is called the "Monte-Carlo Method." It is assumed that each of the 76 dipoles has the same preassigned probability P of being on, and the probability that each dipole is off is 1 - P. A particular activation state is obtained by generating an independent random number Ni(O
