- Email: [email protected]
Estimation of respiratory parameters by the method of covariance ratios
COMPUTERS
AND
BIOMEDICAL
RESEARCH
7,21-39
(1974)
Estimation of Respiratory Parameters by the Method of Covariance Ratios ROB ROY* Systems Engineering Division, Rensselaer Polytechnic Institute, Troy, New York
SAMUELR. POWERS,JR. Albany Medical College, Albany, New York
AND WILLIAM R. KIMBALL Rensselaer Polytechnic Institute, Troy, New York and Albany Medical College, Albany, New York Received June 11,1973
INTRODUCTION
Recent evidence (2, 4, 19, 21, 22, 23, 25) has dramatically pointed out that the survival of posttraumatic patients is highly dependent on the maintainence of adequate respiratory function. Consequently, the modeling of respiratory mechanics and the determination of the parameters of these models becomes of increasing importance. There have been numerous studies of respiratory mechanics which have yielded models ranging from simple first-order linear models (10,17,27) to highly complex nonlinear models (5,24). However, little attention has been placed upon the worthiness of the numbers which have been derived by matching the model to actual physical behavior. Indeed, it makes IittIe sense in a clinical setting to present the physician with a set of numbers describing a complex model if these numbers do not contribute to his understanding of the patient’s condition. In fact, these noncontributory numbers may actually obscure the fundamental changes in respiratory function. Early attempts at respiratory parameter determination used only a simple firstorder linear model, by which dynamic compliance and inelastic airway resistance * Special Fellow in Respiratory Physiology, National (on leave). 21 Copyright 0 I974 by Academic Press, Inc. All rights of reproduction Printed in Great Britain
in any form
reserved.
Institute of General Medical
Sciences.
22
ROY,
POWERS,
JR.,
AND
KIMBALL
were determined. The equation used for these models was (6)
P=K,V+K,li,
(I?
where P = transpulmonary pressure, V = tidal volume, ri = airflow, K1 = reciprocal of dynamic compliance, and K2 = inelastic airway resistance. It should be pointed out that the “constants” K1 and K2 are really not constants. In a biological system, these numbers are in fact random variables, and the number assigned to a particular variable is the estimate of the mean of this variable. Consequently, the assigning of a number to these constants should be done on a proper statistical basis. For example, one procedure (8,18) which is often used to determine K, is to divide P by I/ when ri = 0. Not only is the procedure statistically incorrect, but it suffers from the usual problems associated with the determination of accurate zero crossings in the presence of noise. However, much work was performed on this simple model by these early investigators using analog computers (7, 22, 13, 26). The use of analog computation had the advantage that it was easy, inexpensive, and provided an on-line parameter determination. Later, the gradual replacement of analog computers by digital computers provided a much more powerful tool for the modeling of biological systems (IZ). As a result, the nonlinear aspects of respiratory mechanics were modeled, and the number flood appeared. Since this study was associated with a clinical unit where the numbers from a respiratory model were to be used to determine therapeutic maneuvers, it was of importance to use the simplest possible model on both a physiological and computational basis. Therefore, it was important to determine whether actual patient data contained sufficient information to warrant higher order models. This question was answered by a comparison of least-squares analysis on a set of proposed models (6) which were of the form
P = 5 Kij‘(tj). i=l
(2)
Appendix A presents the fundamentals of least-squares analysis. To some readers, this may appear as another reinvention of the wheel; however, the presentation made serves a dual function. First, the pictorial approach to least squares points out what can and cannot be accomplished by this procedure. Second, it shows that the procedure that is recommended, that of covariance ratios, can be viewed as a least-squares procedure. The data used in this analysis consisted of the original raw data plus data sets which were formed by Fourier decomposition and selected reconstruction. This enabled a comparison of the effects of filtering on the parameter values. In addition, it provided a means for examining the models on a frequency by frequency basis. An analysis of parameter determination on the basis of single frequency information is presented, indicating the limitations and possible applications.
ESTIMATION OF RESPIRATORY PARAMETERS
23
Finally, the new idea, that of covariance ratios, is described. Although simple in description, it relies on the key point that the covariance between a real bandlimited signal and its derivative is identically zero. This point is fully developed in Appendix B. The advantages of this technique are fully discussed in the section covering this topic. However, it can be summarized by stating that the use of covariance ratios provides a simple, statistically correct procedure for parameter determination that can easily be programmed digitally, on an analog computer, or even in special-purpose hardware. Rather than overtax this paper with tables and graphs from a sequence of patients, data from a single patient is presented for illustrative purposes. It is hoped that this approach will be accepted in the spirit in which it is offered, that of clarification. This paper is basically mathematical, designed to place conceptual ideas on a firm mathematical basis. One might view it as being half tutorial and half new technique. Therefore, the data presented is meant to serve these purposes and not overwhelm them. DATA
COLLECTION
The data used in this study was obtained from patients who were being treated in the Trauma Unit at Albany Medical Center Hospital. This unit is a special one-bed section which was set up to study the pulmonary function aspects of severe trauma. A joint medical-engineering team consisting of personnel from both the Albany Medical Center and Rensselaer Polytechnic Institute are responsible for the operation and direction of this unit. In general, the complex consists of instrumentation requisite to maintain ventilatory support and measure pertinent physiological parameters. The interrelationships between the measured parameters and the patients treatment and progress are determined by analysis of the raw data by a direct coupled PDP-15 digital computer. In particular, the measurement of respiratory flow, volume, and pressure is performed as follows. Respiratory gas flow is measured with a bag-box (9) connected to a wedge spirometer. This combination allows the patient to breath air from a flexible bag and back out into the box surrounding the bag. Inspired and expired air do not mix, and yet the volume remains constant. Volume changes due to each breath actuate the wedge spirometer which contains a transducer, which measures both volume and flow signals. Transpulmonary pressure is measured by means of an esophageal balloon (to measure intrapleural pressure) and a differential strain gauge at the mouth of an endotracheal tube which is connected to the bag-box. These three measurements, pressure, flow, and volume, were used to determine the model of respiratory mechanics. Although most of the patients in this unit receive ventilatory assistance from an Ohio respirator, the data used in this study was from spontaneously breathing patients. This enabled a comparison between an inspiratory model and an expiratory model.
24
ROY, POWERS, JR., AND KIMBALL
The sampling rate used was 40 samples/set for all data channels. Since the data has 99 % of its energy located below 1 Hz, this sampling rate is more than adequate to prevent data aliasing. Spectral Energy Distribution
Both pressure and volume waveforms were decomposed into spectral components by the use of the discrete (3) Fourier Transform (DFT) N;l
F(kc0,) = &J’(n) e-j(2n;N)nk,
where o,, = fundamental angular frequency = 27c/N, N = number of discrete sample points, k = index of frequency component, n = index of sample number, F(kq,) = kth complex Fourier component, f(n) = nth discrete sample point. A Fast Fourier Transform (FFT) subroutine was used to implement this decomposition. One method of displaying the spectrum is shown in Fig. 1. This graph shows the cumulative energy distribution plotted versus the harmonic component
(,,),,
,,,,,,,,,,,,,,,,,,, 10
20 HARMONlC
FIG. 1. Energy distribution
COMPONENT
30 (40/1024
40 HERTZI
of pressure and volume waveforms.
index. The total energy of the set of samples used for the FFT was taken to be unity, and the energy of each frequency component was correspondingly scaled and cumulatively added, starting at the first harmonic. The cumulative plot is shown as a percentage of the total energy in the waveform. Any zero frequency component was removed prior to the analysis. The plots show that most of the energy in both pressure and volume waveforms exists between the 12th and 16th harmonic, with the 14th harmonic containing the
ESTIMATION
OF RESPIRATORY
PARAMETERS
25
largest energy component. The 14th harmonic of the pressure waveform contains 60% of the total pressure waveform energy, and the 14th harmonic of the volume waveform contains 71 ‘A of the total volume waveform energy. Note that the term “harmonic” is in reference to the total number of sample points chosen. For this study, 1024 points were selected. Thus, each harmonic is a multiple of (40/1024) Hz, since the original data was sampled at 40 Hz. The 14th harmonic represents a frequency of (14)(40)/1024, or 0.546 Hz. Least-Squares
Analysis
The least-squares procedure detailed in Appendix A was used to find the Ki coefficients for the following models : (1) P=KlV+K,li+K,(v)2+Kqii
(4a)
(2) P=K,
V+K,
ViK@Y
(4b)
(3) P=K,
Vv+K,
3
(4C)
The data used for analysis were the original P and V waveforms and a similar set which was reconstructed from a truncated discrete Fourier Transform. In effect, the use of this reduced data is equivalent to filtering the higher frequency content of the individual waveforms. The basic reconstructed waveforms used the first 17 discrete components of the DFT, each discrete component being spaced (40/1024) Hz apart. For the models where a (p)’ term was used, the pressure waveform was reconstructed using 34 discrete components, since the squared term will produce all discrete components up to twice the upper frequency component of the band-limited signal. The original data was also analyzed for each model as a comparison. The data was analyzed on a breath by breath basis, dividing each breath into an inspiratory and an expiratory set. Start of inspiration was defined as the point .where li = 0, v > 0. Start of expiration was defined as the point where p = 0, P < 0. Due to the noise in the original data, a separate subroutine was used for testing and interpolation in the vicinity of zero flow to determine the start of either inspiration or expiration. Analysis was performed on a cumulative basis, Each new breath cycle was added to the data from the previous breath cycles (separated into inspiration and expiration), and a running cumulative least-squares analysis was performed. Essentially, this is a “growing memory” form of model matching. The results of this analysis are condensed into Table 1 for comparison purposes. The following discussion is in reference to this table. A comparison of Models 1 and 2 points out that (i) the use of the v (or inertia) term does not lower the RMS error and (ii) the use of additional frequency components to account for the theoretical double frequency spectrum associated with the (P)” (or turbulence) term increases the RMS error. If recovery of the K,, Kz coefficients is the primary goal of the analysis, then certainly Model 2 using P/34,
26
ROY,
POWERS,
JR., AND
KIMBALL
V/l7 will provide the same numbers, with the same RMS error, as will Mode1 I using P/34, V/17. A comparison of Models 2 and 3 reveals that the addition of the turbulence (p)’ term, using a P/17, V/17 data base, neither improves the RMS error nor provides a different set of numbers for the K,, Kz coefficient. Based upon the primary goal of obtaining the KI. K2 coefficients, Model 3 is the most appropriate choice. As seen in a later section, use of a model with first-order dynamics allows the use of covariance techniques for the measurement of the desired respiratory parameters. TABLE LEAST-SQUARES
1
ESTIMATES OF RESPIRATORY 10 Breath
--~---
~~ -~~ Imp.
K,
-~~~~ Exp. Insp.
K2
~~~.~_~~ Exp. Imp.
PAUAMETERS
Cycles RMS
K4
k
Exp.
Model I P/34 v/17 P/17 v/17 Original
26.8 18.8 20.3
13.3 19.3 15.5
6.6 4.8 6.4
3.0 4.7 5.0
-5.7 -0.71 -4.7
-4.2 -0.49 -6.5
Model 2 P/34 v/17 P/17 v/17
27.0 19.7
12.9 19.7
6.6 4.8
3.0 4.7
-5.7 -0.78
-4.2 0.46
Model 3 P/17 v/17 Original
19.7 23.3
19.7 18.0
4.5 4.2
4.5 4.6
Insp. -.017 -.091 -0.14
Exp. -.033 -.038 -0.11
error
Imp.
Exp.
0.62 0.34 0.13
0.78 0.35 0.86
0.62 0.34
0.78 0.35
0.34 0.75
0.35 0.88
A point of philosophy should be injected into the discussion at this stage. The use of least-squares analysis does not “prove” that a model is the correct one, nor does it indicate anything concerning the physics or physiology of the biological system under concern. This analysis simply provides a set of numbers indicating the values of coefficients which fit, in a least-squares sense, an assumed differential equation. As a measure of the “fit,” the analysis also provides the RMS error between the model assumed and the actual system. An injudicious choice of a set of models to choose between would lead to a poor model of the actual system, although it may be best among the choices listed. For this reason, the models selected were those which physiologists have considered to be “reasonable” replicas of the physiological system (6, 14, 15). From the numbers shown, it is difficult to justify the selection of one model over another solely on the basis of RMS error. A possible reason for the closeness in RMS errors is that the first two terms of the general form assumed for the model
ESTIMATION
OF RESPIRATORY
PARAMETERS
27
.
... IL “( :“‘ *. ”...
c
qc”c
w. )(L
l L .
: . . l . . . . . ‘ . . . . . . . . . . . . l
f
. *
”
C”
c: C”‘”
* . . * . . . . .
lb,?9
: / : : * . : . r* . . . . : . : : c*
L
*
t”
:
7” El : L* t* C” :. -aa L. .
FIG. 2. Pressure waveform after harmonic filtering. measured filtered data.
C= computed from model; M=
actual
28
ROY, POWERS, JR., AND KIMBALL
.. .. :. .
FIG. 3. Pressure waveform, unfiltered. waveform, unfiltered.
C = computed
from model:
BY = actual measured
ESTIMATION
OF RESPIRATORY
PARAMETERS
29
are the dominant terms and that the remaining terms represent “noise variables” which contribute little to the data fitting. As indicated in Appendix A, this would be illustrated by additional vectors which do not lie in the space which generates the fitted variable (II). Consequently, partly on the basis of simplicity, Model 3 is selected as the model of choice. Figures 2 and 3 illustrate how Model 3 follows the actual data. The actual pressure waveform is shown as the sequence of points marked “M”, and the computed, or model, pressure waveform is the sequence of points marked “c”. Where both waveforms occupy the same point, a “period” is shown. Figure 2 shows the matching for the reduced or filtered waveform, and Fig. 3 shows the matching for the raw data. dc values have been removed from both waveforms. From a visual standpoint, this can be considered to be good matching. Consequently, it can be assumed that this is not the best of a bad set of models but is indeed a reasonably good model. Single-Frequency
Analysis
Several investigators (I, 16, 20) have used the numbers associated with the principal or fundamental frequency of the pressure-volume waveform to determine the respiratory parameters of the first-order model P=K,V+K,J?
(5)
The following analysis shows the conditions under which this procedure can be considered valid. Consider, then, a real band-limited waveform which can be represented by a finite Fourier series
x(t) = k=NF
Ckuejkoot.
(6)
Assume that the true volume signal is the sum of x(t) and a wide spectrum noise signal n,(t) c(t) = x(t) + n,(t). (7) Similarly,
Substitution
the pressure waveform can be assumed to be of the form
of Eqs. (6), (7), and (8) into Eq. (5) yields eikwot + jK, kw,
2 C,, ejkwot + n,(t) - rip(t) k=-N
= (KI +jKz km,)
5 C,, ejkwOt f n,(t) - n,(t). k=-N
(9?
The question to be resolved is what are the conditions under which knowledge of only one pair of complex Fourier coefficients can yield both unknowns KI and Kz ?
30
ROY,
POWERS,
JR.,
AND
Klh4BALL
Assume that the coefficients of the mth harmonic are known. The complex frequencies ejkwot form a orthonormal set of basis vectors, such that
where ( , ) represents the scalar (or inner) product
Consequently, elmDot yields
taking the inner product of both sides of Eq. (9) with respect to c,, = c,,(K,
+,iK, ItiC&) -t i@uot, (n,(t) - f’$,(t))::,.
(11)
It is clear that the accuracy of the determination of Ki, K2 is highly dependent upon the component of noise in the specific frequency moo. If the noise energy is small with respect to the signal energy at the frequency of measurement, then Kl and K2 can be determined as follows : and
(12)
where Re( } = Real Part, Im { ) = Imaginary
Part.
In an actual data analysis, the noise terms no(t) and n,(t) are included in the complex Fourier coefficients. Therefore, the use of Eq. (12) always implies an error in the calculation. To minimize this error, the frequency component with the largest signal energy should be used for the calculation. Use of harmonic components with small signal energy can lead to meaningless numbers for K,, K2. For the illustrative data run shown in this paper, the values of K,, K2 obtained by the single frequency analysis, using the 14th harmonic (0.546 Hz), are K1 = 19.7 and K2 = 4.4. Other investigators have used a single-frequency analysis to determine these parameters but have done so with a highly error-prone technique, that of examining zero crossings. Similarly zero crossings can be directly applied to Eq. (5) to determine K,, K,. However, the use of zero crossings in a noisy signal can lead to erroneous results. Even the phase shift inherent in the filter can distort the time relationships between pressure and volume signals. This time distortion is also present in linear phase filters since pressure and volume signals do not have identical frequency spectrums.
ESTIMATION
OF RESPIRATORY
31
PARAMETERS
Parameter Determination by Variance Ratios Assume that the differential equation governing the relationship and volume is given by Eq. (1) :
between pressure
P=K,V+KJ?
(1)
Assuming that the dc terms have been removed, the covariance between pressure and volume is o&, = E{PV} = E{K, Vz + K2 liV}
(13)
= Kl E{ V’} + K2 E{ J+‘) = Kl o;, + K2 L&. However, the covariance between a real band-limited identically zero (see Appendix B). Consequently, K 1 d&&
(14) signal and its derivative is
covariance between pressure and volume signal mean square value of volume signal
(15a)
covariance between pressure and flow signal mean square value of flow signal
(15b)
In a similar manner,
,g&=
o;p
There are several advantages in using covariance ratios to estimate respiratory parameters. First, since both Kl and K2 are determined by dividing a covariance by a mean-squared value (which is the total signal power), the calculations are inherently normalized to signal strength. Second, filtering is not required. The use of mean values provides a measure of noise immunity provided the noise is uncorrelated with the signal. Third, zero crossings are not required. As stated earlier, zero-crossing techniques are at best treacherous and should be avoided if possible. Fourth, the calculations are simple. Although the method of single-frequency analysis of the spectral component with maximum energy and the method of leastsquares analysis will provide reliable numbers, there is a certain minimum amount of calculation which must be performed for each method. The single-frequency analysis requires a Fourier decomposition followed by maximum detection and then simple ratios. The least-squares procedure involves data storage, matrix multiplication, and matrix inversion. With judicious programming, these techniques can produce numbers quite rapidly; however, a digital computation facility is required. Alternatively, although digital computation is desirable for the covariance ratio method, it is not necessary. The calculations required can be performed by analog computation, providing an on-line running estimate of the parameters K,, K2. To illustrate this point, examine the block diagram shown in Fig. 4. The signals P and V are passed through an R-C circuit which removes the mean .T
32
ROY,
POWERS,
JR., AND
KIMBALL
or dc value from the signals. The signal pihas zero mean. Alternatively, an integrator and a subtractor could be used to remove the mean from each of the signals. The zero mean signals V’ and P are then squared and cross-multiplied by P to form V2, PV, Pp, and (lj)l. Each of these signals are then passed through “leaky integrators”. Basically, these are operational amplifier integrators with a time constant adjusted so that past data can “leak off” exponentially and new data incorporated into the integrator. The output of the integrators provides a running measure of the integral of the input without saturation of the amplifier. The formation of K,, K2 requires the ratios of expected value which are produced by integration with respect to time and then division by the integration time. Since C v
” 5
-II
(
2 a
I2
s+a
R
A DIVIDE
PV
MULTIPLY
K,
T
0
s+a
A
p-it P;
0
s+a
3 UIVIDE
T FIG.
I
7
4. Analog computation
vh2
a
K2
?
s+a
of respiratory parameters.
only ratios are formed, division of each integral by the time of observation is not required since the time baseappearsin both numerator and denominator and therefore cancel. The output of the division circuits are the desired numbers Kl, Kz. If K,, Kz vary with respect to time, this arrangement of exponential averaging will allow a continuous tracking of the changing values. This circuit represents a method of obtaining on-line measurementof the respiratory parameters Kl and Kz. As an example, a running estimate of parameter K1 at the end of successive breaths is listed in Table 2. Also shown is the sameestimate as provided by least squares, separated into inspiration and expiration. These estimates were obtained on the raw unfiltered data. At the conclusion of ten breaths, there is only a 3 % difference between the values obtained by covariance ratios and that obtained by least squares.This is to be expected, for as shown in Appendix A, the covariance ratio technique is actually a least-squaresprocedure simplified by the fact that the covariance of a real band-limited signal and its time derivative is identically zero. The small differences shown are due to both the separation into inspiration and expiration and a nonzero value of COV{XJ?} over a finite time interval. As shown
33
ESTIMATION OF RESPIRATORY PARAMETERS
in Table 2, this difference is small, showing that the covariance ratio technique is a viable method for obtaining respiratory parameters. TABLE CUMULATIVE
End of breath
1 2 3 4 5 6 7 8 9 10
2
ESTIMATE OF PARAMETER
Covariance
ratios
19.3 20.4 20.2 20.1 20.2 20.7 20.9 20.8 20.4 20.3
K,
Least squares Insp. Exp. 19.4 17.8 18.5 19.2 19.5 19.3 19.5 19.7 19.7 19.7
18.2 18.1 19.1 19.9 19.8 19.9 19.9 19.9 19.8 19.7
CONCLUSIONS
A new technique has been presented for the estimation of the respiratory parameters associated with a first-order linear model of respiratory mechanics. This method, that of covariance ratios, is shown to be a modification of least-squares analysis. However, this modification allows for a considerably simplified (but mathematically correct) estimation procedure which can be programmed with a small analog computer. Consequently a real-time on-line least-squaresestimation of respiratory parameters can be provided without the necessitydigital computation and data storage. APPENDIX
A: LEAST-SQUARES
ANALYSIS
(3)
Consider a systemwith a single scalar output y which is to be approximated by a linear combination of a set of linearly independent variables xl, x2,. . ., x,. Call the approximation j. Thus L’~~=k~x~+k,x,+.~.+k,x,. (A-1) Assume that A4 successivemeasurementsof both y and the variables xi are made. A set of equations could then be formed : ~‘(1)~=:(1)=k,x,(l)+k,x,(l)+~~~+k,x,(l) y(2) z j(2) = kl x,(2) + k2 x2(2) + . . . + kN x,,(2) .... ..... ..... ..... ..... .... ...... ... . ~‘(M)~~(~)=k~x,(M)+k,x,(M)+..-+kk,x,(M)
(A-2)
34
ROY,
POWERS,
JR.,
AND
KIMBALL
To simplify the notation, let each column variable be designated as a vector such that (A-3)
y~:=kklx,$-k,x,+.‘.+li,~x,.
Alternatively,
a matrix X can be formed consisting of the column vectors x=
[XlX2...Xh]
(A-4)
Y=Xk.
r A-5)
such that
The error between the true value y and the approximation vector E:
y is given by the error
g=y-f=y-Xk.
(A-6)
The problem is to determine a vector k which minimizes the error vector in a least-squares sense. A pictorial appreciation of the procedure involved in leastsquares approximation is provided in Fig. A-l. For the purpose of simplicity, only two variables, x1 and x2, will be used, thus only two vectors x1 and x2 are shown.
FIG.
A-l.
Orthogonal projection:
the principle of least-squares fitting.
These two vectors x1 and x2 form a plane. The measurement vector y is assumed to lie outside of the plane formed by x1 and x2. The approximation to y is formed by a linear combination of x1 and x2 : f=klxl+k,x,.
t A-7) This approximation must lie in the plane formed by x1 and x2. From basic geometric considerations, the best approximation to y is given by an orthogonal projection of y onto the x1, x2 plane. This is the orthogonal projection theorem, which forms the basis for least-squares approximation. If J; is formed as shown, then 9 is a leastquares approximation to y. Note that this projection requires that the error vector
ESTIMATION
OF RESPIRATORY
35
PARAMETERS
E is orthogonal to the plane x1, x2 and, thus, is orthogonal to both x1 and x2. This concept is easily extended to higher dimensions. Furthermore, note that perfect matching of y is possible only when the measurement y lies within the same space subtended by the vectors xi. Therefore, it is extremely important that the independent variables x1, x,, . . ., X~ form as complete a set of measurements as possible and that a sufficient number of measurements be made (M > N). The presence of measurement noise will of course add an additional dimension to the space such that y cannot possible lie in the subspace subtended by the xi vectors. Continuing with the analysis, observe that E must be simultaneously orthogonal to all of the xi vectors. This requirement can be written as XiT&=O
for all i,
or XTE = 0.
(A-8)
From Eq. (A-6), XTz=O=XTy-X=Xk or
(A-9)
k = (X’ X)-l XT y.
Equation (A-9) is the basic equation of least-squares analysis. It is instructive to note the dimensions of each of the terms of this equation : k=Nx
g=Mx
N variables, M successive measurements,
1, 1,
X=MxN, (XT X)-l = N x N.
For the case illustrated, where N = 2, the vector k is of order (2 x 1). Since the inverse of a 2 x 2 matrix is given by
A-’ =
(A-10)
the vector k is given by igl
x22@?
C ~~(9~49- 2 x1(9
& x12W C x2W.W
x2(i)
- C x1(9x2(i) C xlG).W
C xl”(i) 2 x2’(i) - CCxl(~)x2W”.
It is important
I-
C x2G>v(i)
(A-l 1)
to note that if the cross product term is zero,
jtl xl(i) x2(i)= 0,
(A-12)
ROY,POWERS,JR.,AND
36
then the values fork,
KlMBALL
and k, are simply
(A-13a)
M 2 x,(i)y(i k, = k&e-~-
j
‘;; x,2(i) ,Tl
Equation (A-l 3) will be used to show that the covariance ratio technique is identical to a least-squares procedure, provided that the covariance between x, and .x2 is identically zero. APPENDIX The cross correlation
B: COVARIANCE RATIO ANALYSIS (3)
function between a signal and its derivative is given by
(B-1)
R&l, t2)= Wt,) qaj. where E{ } represents the expectation operator. In terms of the definition derivative, [B-l] can be rewritten as
x(t, + 6)- x(fJ -.---irm--L--
of a
I)
(B-2)
If x(t) is a stationary signal, the correlation function is dependent only on the difference between tl and t2 and not on their absolute values. Letting t1 - f2 = z. the correlation function can be stated as R,*(z) = -aR,,(r)/aT.
Let x(t) be a real band-limited series
(B-3)
signal which can be represented by a finite Fourier
Ckejkwot.
(B-4)
ESTIMATION OF RESPIRATORY PARAMETERS
37
The correlation function R,,(z) is given by R,,(T) = E
5 Ck &k%t + C, ejj.-o(t+r) AYN ck&N 1
=E
5 2 ck Ciei(k+E.)wotejl.wor . (B-5) 1,,-N k=-N 1 Due to the orthogonality of the exponential functions, the above expectation is zero except when k + 1. = 0. Therefore, R,,(T) = 5 C, Cmk(ejkwor + e-jkoor). k=-N
(B-6)
Since x(t) is real, the Fourier coefficients are symmetric about k = 0, thus Ck = C-,. R,,(z) = C,” f 2 : C,” cos kw,. k=l
(B-7):
From [B-3],
R,JT) =-v
=2 ? Ck2(kCo,sin ko, 2). k:l
(B-8)
The zero value of R,,(z) is identically zero: R,,(O) = 0.
(B-9)
The zero value of a correlation function represents a mean-square value: R,,(t) = &x(t) x(t + T)}.
(B-10)
R,,(O) = E{x2(t)) = x’(t).
(B-l I)
For T = 0, For a signal whose mean value is zero, x2(t) = mean-square value = oh = variance of x.
(B-12)
In summary, a:, = c,’ + 2 ;: ck2, kT1
CT:, = 0.
(B-13)
ACKNOWLEDGMENTS This work was performed in conjunction with the Departments of Surgery and Physiology, Albany Medical College of Union University, Albany, NY and the Systems Engineering Division of Rensselaer Polytechnic Institute, Troy, NY. This work was supported by Grant No. GM 15426 from the National Institute of General
38
ROY,
POWERS,
JR.,
AND
KIMBALL
Medical Sciences, General Clinical Research Center Grant No. 5 MO I-RR 00094m, National Institute of General Medical Sciences Special Fellowship No. 5 F03 GM 53156-02, and Project Themis Contract No. DAAB07-69-C-0365.
REFERENCES 1. ALBRIGHT, C. D., AND BONDURANT, S. Some effects of respiratory frequency on pulmonary mechanics. J. Clin. Invest. 44(S), 1362-1370 (1965). 2. ASHBAUGH, D. G. Continuous positive pressure breathing (CPPB) in adult respiratory distress syndrome. J. Thoruc. Curdio. Surg. 57,31-41 (1969). 3. DERUSSO, P. M., ROY, R. J., AND CLOSE, C. M. “State Variables for Engineers.” John Wiley and Sons, Inc., New York, 1965. 4. FALKE, K. .I. Ventilation with end-expiratory pressure in acute lung disease. J. Clin. Invest. 51(9),
2315-2323
(1972).
5. FEINBERG, B. N. Parameter estimation: A diagnostic aid for lung diseases. Instrumentation Technol. 4Cb46 (1970). 6. FENN, W. O., AND H. RAHN (Ed.), Handbook of Physiology, Section 3 : Respiration, Vol. 1, pp. 363-376,411428,451462. Washington, D.C.: American Physiological Society, 1964. 7. FLETCHER,G., AND BELLVILLE, J. W. On line computation of pulmonary compliance and work of breathing. J. A&. Physiol. 21(4), 1321-1327 (1966). 8. FRY, D. L., AND HYATT, R. E. A unified analysis of the relationship between pressure, volume and gasflow in the lungs of normal and diseased human subjects. Amer. J. Med. 672-689 (1960). 9. GISSER,D. G., NEWELL, J. C., AND POWERS,S. R., JR. A refined bag-box spirometer. J. .4ss. Advan. Med. Znstr. 6, 167 (1972). 10. JODAT, R. W. Simulation of respiratory mechanics. BiophE’s. J. 6,773-785 (1966). 11. LEWIS, F. J. Analysis of respiration by an on-line digital computer system. Proc. Amer. Surg. Ass. 211-221 (1966). 12. MCILROY, M. B. A new method for measurement of compliance and resistance of lungs and thorax. J. Appl. Physiol. 18(2), 424-427 (1963). 13. MCWILLIAMS, R., AND ADAMS, A. H. Analogue computer for study of breathing mechanics. Med. Electron. Biol. Eng. 1,353 (1963). 14. MEAD, J., AND WH~~TENBERGER,J. L. Physical properties of human lungs measured during spontaneous respiration. J. Appl. Physiol. 5,779 (1953). 15. MEAD, J. Mechanical properties of lungs. Physiol. Rev. 41,281 (1964). 16. MEDITCH, J. S., AND STOLL, P. J. On respiratory system parameter estimation. 1970 Joint Automatic Control Conference, Paper 16-D, pp. 383-390. 17. MILLER, C. W. “A Least Squares Fit Method of Measuring Lung and Thorax Compliance in a Continuous Patient Monitoring Systems.” 22ACEMB, Paper 10.13, p. 124,197O. 18. MINDLIN, S. “On-Line Computation of Respiratory Characteristics.” EAI Applications Bulletin 953055, 1969. 19. MONACO, V., BURDGE,R., NEWELL, J., LEATHER, R., AND POWERS,S. R., JK. Clinical significance of the functional residual capacity in the post-injury state. Sarg. Forum 22,42-43 (1971). 20. OSTRANDER,L., AND CHESTER,E. H. “Frequency Analysis of Pulmonary Pressure-Flow Data.’ 23ACEMB, Paper 16.11, p. 208,197O. 21. PETERS,R. M. Work of breathing. J. Trauma 8(5), 915-921 (1968). 22. PONTOPPIDAN, M., GEFFIN, B., AND LOWENSTEIN, E. Acute respiratory failure in the adult. (Three Parts) New Eng. J. Med. 287: 14,15,16,690-698,743-752,799-806,1972. 23. PRICE, A., AND POWERS, S. R., JR. “The Effect of a Period of Acute Hemorrhagic Shock in the Respiratory Compliance of the Dog,” 21st ACEMB, Paper 42.1, 1968.
ESTIMATION
OF RESPIRATORY
PARAMETERS
39
24. ROSSING, R. G. “A Digital Computer Program for the Estimation of Respiratory Mechanical Parameters,” Proceedings of Symposium on Biomedical Engineering, Marquette University, pp. 134-136,1966. 25. SEALY, W. C. The lung in hemorrhagic shock. J. Trauma 8,774-781(1968). 26. STACY, R. W., AND PETERS, R. M. Computations of respiratory mechanical parameters. “Computers in Biomedical Research.” Chapter 11, pp. 269-288, Academic Press, New York, 1965. 27. WALD, A. A computer system for respiratory parameters. Comp. Biomed. Res. 2,41 l-429,1969.
AND
BIOMEDICAL
RESEARCH
7,21-39
(1974)
Estimation of Respiratory Parameters by the Method of Covariance Ratios ROB ROY* Systems Engineering Division, Rensselaer Polytechnic Institute, Troy, New York
SAMUELR. POWERS,JR. Albany Medical College, Albany, New York
AND WILLIAM R. KIMBALL Rensselaer Polytechnic Institute, Troy, New York and Albany Medical College, Albany, New York Received June 11,1973
INTRODUCTION
Recent evidence (2, 4, 19, 21, 22, 23, 25) has dramatically pointed out that the survival of posttraumatic patients is highly dependent on the maintainence of adequate respiratory function. Consequently, the modeling of respiratory mechanics and the determination of the parameters of these models becomes of increasing importance. There have been numerous studies of respiratory mechanics which have yielded models ranging from simple first-order linear models (10,17,27) to highly complex nonlinear models (5,24). However, little attention has been placed upon the worthiness of the numbers which have been derived by matching the model to actual physical behavior. Indeed, it makes IittIe sense in a clinical setting to present the physician with a set of numbers describing a complex model if these numbers do not contribute to his understanding of the patient’s condition. In fact, these noncontributory numbers may actually obscure the fundamental changes in respiratory function. Early attempts at respiratory parameter determination used only a simple firstorder linear model, by which dynamic compliance and inelastic airway resistance * Special Fellow in Respiratory Physiology, National (on leave). 21 Copyright 0 I974 by Academic Press, Inc. All rights of reproduction Printed in Great Britain
in any form
reserved.
Institute of General Medical
Sciences.
22
ROY,
POWERS,
JR.,
AND
KIMBALL
were determined. The equation used for these models was (6)
P=K,V+K,li,
(I?
where P = transpulmonary pressure, V = tidal volume, ri = airflow, K1 = reciprocal of dynamic compliance, and K2 = inelastic airway resistance. It should be pointed out that the “constants” K1 and K2 are really not constants. In a biological system, these numbers are in fact random variables, and the number assigned to a particular variable is the estimate of the mean of this variable. Consequently, the assigning of a number to these constants should be done on a proper statistical basis. For example, one procedure (8,18) which is often used to determine K, is to divide P by I/ when ri = 0. Not only is the procedure statistically incorrect, but it suffers from the usual problems associated with the determination of accurate zero crossings in the presence of noise. However, much work was performed on this simple model by these early investigators using analog computers (7, 22, 13, 26). The use of analog computation had the advantage that it was easy, inexpensive, and provided an on-line parameter determination. Later, the gradual replacement of analog computers by digital computers provided a much more powerful tool for the modeling of biological systems (IZ). As a result, the nonlinear aspects of respiratory mechanics were modeled, and the number flood appeared. Since this study was associated with a clinical unit where the numbers from a respiratory model were to be used to determine therapeutic maneuvers, it was of importance to use the simplest possible model on both a physiological and computational basis. Therefore, it was important to determine whether actual patient data contained sufficient information to warrant higher order models. This question was answered by a comparison of least-squares analysis on a set of proposed models (6) which were of the form
P = 5 Kij‘(tj). i=l
(2)
Appendix A presents the fundamentals of least-squares analysis. To some readers, this may appear as another reinvention of the wheel; however, the presentation made serves a dual function. First, the pictorial approach to least squares points out what can and cannot be accomplished by this procedure. Second, it shows that the procedure that is recommended, that of covariance ratios, can be viewed as a least-squares procedure. The data used in this analysis consisted of the original raw data plus data sets which were formed by Fourier decomposition and selected reconstruction. This enabled a comparison of the effects of filtering on the parameter values. In addition, it provided a means for examining the models on a frequency by frequency basis. An analysis of parameter determination on the basis of single frequency information is presented, indicating the limitations and possible applications.
ESTIMATION OF RESPIRATORY PARAMETERS
23
Finally, the new idea, that of covariance ratios, is described. Although simple in description, it relies on the key point that the covariance between a real bandlimited signal and its derivative is identically zero. This point is fully developed in Appendix B. The advantages of this technique are fully discussed in the section covering this topic. However, it can be summarized by stating that the use of covariance ratios provides a simple, statistically correct procedure for parameter determination that can easily be programmed digitally, on an analog computer, or even in special-purpose hardware. Rather than overtax this paper with tables and graphs from a sequence of patients, data from a single patient is presented for illustrative purposes. It is hoped that this approach will be accepted in the spirit in which it is offered, that of clarification. This paper is basically mathematical, designed to place conceptual ideas on a firm mathematical basis. One might view it as being half tutorial and half new technique. Therefore, the data presented is meant to serve these purposes and not overwhelm them. DATA
COLLECTION
The data used in this study was obtained from patients who were being treated in the Trauma Unit at Albany Medical Center Hospital. This unit is a special one-bed section which was set up to study the pulmonary function aspects of severe trauma. A joint medical-engineering team consisting of personnel from both the Albany Medical Center and Rensselaer Polytechnic Institute are responsible for the operation and direction of this unit. In general, the complex consists of instrumentation requisite to maintain ventilatory support and measure pertinent physiological parameters. The interrelationships between the measured parameters and the patients treatment and progress are determined by analysis of the raw data by a direct coupled PDP-15 digital computer. In particular, the measurement of respiratory flow, volume, and pressure is performed as follows. Respiratory gas flow is measured with a bag-box (9) connected to a wedge spirometer. This combination allows the patient to breath air from a flexible bag and back out into the box surrounding the bag. Inspired and expired air do not mix, and yet the volume remains constant. Volume changes due to each breath actuate the wedge spirometer which contains a transducer, which measures both volume and flow signals. Transpulmonary pressure is measured by means of an esophageal balloon (to measure intrapleural pressure) and a differential strain gauge at the mouth of an endotracheal tube which is connected to the bag-box. These three measurements, pressure, flow, and volume, were used to determine the model of respiratory mechanics. Although most of the patients in this unit receive ventilatory assistance from an Ohio respirator, the data used in this study was from spontaneously breathing patients. This enabled a comparison between an inspiratory model and an expiratory model.
24
ROY, POWERS, JR., AND KIMBALL
The sampling rate used was 40 samples/set for all data channels. Since the data has 99 % of its energy located below 1 Hz, this sampling rate is more than adequate to prevent data aliasing. Spectral Energy Distribution
Both pressure and volume waveforms were decomposed into spectral components by the use of the discrete (3) Fourier Transform (DFT) N;l
F(kc0,) = &J’(n) e-j(2n;N)nk,
where o,, = fundamental angular frequency = 27c/N, N = number of discrete sample points, k = index of frequency component, n = index of sample number, F(kq,) = kth complex Fourier component, f(n) = nth discrete sample point. A Fast Fourier Transform (FFT) subroutine was used to implement this decomposition. One method of displaying the spectrum is shown in Fig. 1. This graph shows the cumulative energy distribution plotted versus the harmonic component
(,,),,
,,,,,,,,,,,,,,,,,,, 10
20 HARMONlC
FIG. 1. Energy distribution
COMPONENT
30 (40/1024
40 HERTZI
of pressure and volume waveforms.
index. The total energy of the set of samples used for the FFT was taken to be unity, and the energy of each frequency component was correspondingly scaled and cumulatively added, starting at the first harmonic. The cumulative plot is shown as a percentage of the total energy in the waveform. Any zero frequency component was removed prior to the analysis. The plots show that most of the energy in both pressure and volume waveforms exists between the 12th and 16th harmonic, with the 14th harmonic containing the
ESTIMATION
OF RESPIRATORY
PARAMETERS
25
largest energy component. The 14th harmonic of the pressure waveform contains 60% of the total pressure waveform energy, and the 14th harmonic of the volume waveform contains 71 ‘A of the total volume waveform energy. Note that the term “harmonic” is in reference to the total number of sample points chosen. For this study, 1024 points were selected. Thus, each harmonic is a multiple of (40/1024) Hz, since the original data was sampled at 40 Hz. The 14th harmonic represents a frequency of (14)(40)/1024, or 0.546 Hz. Least-Squares
Analysis
The least-squares procedure detailed in Appendix A was used to find the Ki coefficients for the following models : (1) P=KlV+K,li+K,(v)2+Kqii
(4a)
(2) P=K,
V+K,
ViK@Y
(4b)
(3) P=K,
Vv+K,
3
(4C)
The data used for analysis were the original P and V waveforms and a similar set which was reconstructed from a truncated discrete Fourier Transform. In effect, the use of this reduced data is equivalent to filtering the higher frequency content of the individual waveforms. The basic reconstructed waveforms used the first 17 discrete components of the DFT, each discrete component being spaced (40/1024) Hz apart. For the models where a (p)’ term was used, the pressure waveform was reconstructed using 34 discrete components, since the squared term will produce all discrete components up to twice the upper frequency component of the band-limited signal. The original data was also analyzed for each model as a comparison. The data was analyzed on a breath by breath basis, dividing each breath into an inspiratory and an expiratory set. Start of inspiration was defined as the point .where li = 0, v > 0. Start of expiration was defined as the point where p = 0, P < 0. Due to the noise in the original data, a separate subroutine was used for testing and interpolation in the vicinity of zero flow to determine the start of either inspiration or expiration. Analysis was performed on a cumulative basis, Each new breath cycle was added to the data from the previous breath cycles (separated into inspiration and expiration), and a running cumulative least-squares analysis was performed. Essentially, this is a “growing memory” form of model matching. The results of this analysis are condensed into Table 1 for comparison purposes. The following discussion is in reference to this table. A comparison of Models 1 and 2 points out that (i) the use of the v (or inertia) term does not lower the RMS error and (ii) the use of additional frequency components to account for the theoretical double frequency spectrum associated with the (P)” (or turbulence) term increases the RMS error. If recovery of the K,, Kz coefficients is the primary goal of the analysis, then certainly Model 2 using P/34,
26
ROY,
POWERS,
JR., AND
KIMBALL
V/l7 will provide the same numbers, with the same RMS error, as will Mode1 I using P/34, V/17. A comparison of Models 2 and 3 reveals that the addition of the turbulence (p)’ term, using a P/17, V/17 data base, neither improves the RMS error nor provides a different set of numbers for the K,, Kz coefficient. Based upon the primary goal of obtaining the KI. K2 coefficients, Model 3 is the most appropriate choice. As seen in a later section, use of a model with first-order dynamics allows the use of covariance techniques for the measurement of the desired respiratory parameters. TABLE LEAST-SQUARES
1
ESTIMATES OF RESPIRATORY 10 Breath
--~---
~~ -~~ Imp.
K,
-~~~~ Exp. Insp.
K2
~~~.~_~~ Exp. Imp.
PAUAMETERS
Cycles RMS
K4
k
Exp.
Model I P/34 v/17 P/17 v/17 Original
26.8 18.8 20.3
13.3 19.3 15.5
6.6 4.8 6.4
3.0 4.7 5.0
-5.7 -0.71 -4.7
-4.2 -0.49 -6.5
Model 2 P/34 v/17 P/17 v/17
27.0 19.7
12.9 19.7
6.6 4.8
3.0 4.7
-5.7 -0.78
-4.2 0.46
Model 3 P/17 v/17 Original
19.7 23.3
19.7 18.0
4.5 4.2
4.5 4.6
Insp. -.017 -.091 -0.14
Exp. -.033 -.038 -0.11
error
Imp.
Exp.
0.62 0.34 0.13
0.78 0.35 0.86
0.62 0.34
0.78 0.35
0.34 0.75
0.35 0.88
A point of philosophy should be injected into the discussion at this stage. The use of least-squares analysis does not “prove” that a model is the correct one, nor does it indicate anything concerning the physics or physiology of the biological system under concern. This analysis simply provides a set of numbers indicating the values of coefficients which fit, in a least-squares sense, an assumed differential equation. As a measure of the “fit,” the analysis also provides the RMS error between the model assumed and the actual system. An injudicious choice of a set of models to choose between would lead to a poor model of the actual system, although it may be best among the choices listed. For this reason, the models selected were those which physiologists have considered to be “reasonable” replicas of the physiological system (6, 14, 15). From the numbers shown, it is difficult to justify the selection of one model over another solely on the basis of RMS error. A possible reason for the closeness in RMS errors is that the first two terms of the general form assumed for the model
ESTIMATION
OF RESPIRATORY
PARAMETERS
27
.
... IL “( :“‘ *. ”...
c
qc”c
w. )(L
l L .
: . . l . . . . . ‘ . . . . . . . . . . . . l
f
. *
”
C”
c: C”‘”
* . . * . . . . .
lb,?9
: / : : * . : . r* . . . . : . : : c*
L
*
t”
:
7” El : L* t* C” :. -aa L. .
FIG. 2. Pressure waveform after harmonic filtering. measured filtered data.
C= computed from model; M=
actual
28
ROY, POWERS, JR., AND KIMBALL
.. .. :. .
FIG. 3. Pressure waveform, unfiltered. waveform, unfiltered.
C = computed
from model:
BY = actual measured
ESTIMATION
OF RESPIRATORY
PARAMETERS
29
are the dominant terms and that the remaining terms represent “noise variables” which contribute little to the data fitting. As indicated in Appendix A, this would be illustrated by additional vectors which do not lie in the space which generates the fitted variable (II). Consequently, partly on the basis of simplicity, Model 3 is selected as the model of choice. Figures 2 and 3 illustrate how Model 3 follows the actual data. The actual pressure waveform is shown as the sequence of points marked “M”, and the computed, or model, pressure waveform is the sequence of points marked “c”. Where both waveforms occupy the same point, a “period” is shown. Figure 2 shows the matching for the reduced or filtered waveform, and Fig. 3 shows the matching for the raw data. dc values have been removed from both waveforms. From a visual standpoint, this can be considered to be good matching. Consequently, it can be assumed that this is not the best of a bad set of models but is indeed a reasonably good model. Single-Frequency
Analysis
Several investigators (I, 16, 20) have used the numbers associated with the principal or fundamental frequency of the pressure-volume waveform to determine the respiratory parameters of the first-order model P=K,V+K,J?
(5)
The following analysis shows the conditions under which this procedure can be considered valid. Consider, then, a real band-limited waveform which can be represented by a finite Fourier series
x(t) = k=NF
Ckuejkoot.
(6)
Assume that the true volume signal is the sum of x(t) and a wide spectrum noise signal n,(t) c(t) = x(t) + n,(t). (7) Similarly,
Substitution
the pressure waveform can be assumed to be of the form
of Eqs. (6), (7), and (8) into Eq. (5) yields eikwot + jK, kw,
2 C,, ejkwot + n,(t) - rip(t) k=-N
= (KI +jKz km,)
5 C,, ejkwOt f n,(t) - n,(t). k=-N
(9?
The question to be resolved is what are the conditions under which knowledge of only one pair of complex Fourier coefficients can yield both unknowns KI and Kz ?
30
ROY,
POWERS,
JR.,
AND
Klh4BALL
Assume that the coefficients of the mth harmonic are known. The complex frequencies ejkwot form a orthonormal set of basis vectors, such that
where ( , ) represents the scalar (or inner) product
Consequently, elmDot yields
taking the inner product of both sides of Eq. (9) with respect to c,, = c,,(K,
+,iK, ItiC&) -t i@uot, (n,(t) - f’$,(t))::,.
(11)
It is clear that the accuracy of the determination of Ki, K2 is highly dependent upon the component of noise in the specific frequency moo. If the noise energy is small with respect to the signal energy at the frequency of measurement, then Kl and K2 can be determined as follows : and
(12)
where Re( } = Real Part, Im { ) = Imaginary
Part.
In an actual data analysis, the noise terms no(t) and n,(t) are included in the complex Fourier coefficients. Therefore, the use of Eq. (12) always implies an error in the calculation. To minimize this error, the frequency component with the largest signal energy should be used for the calculation. Use of harmonic components with small signal energy can lead to meaningless numbers for K,, K2. For the illustrative data run shown in this paper, the values of K,, K2 obtained by the single frequency analysis, using the 14th harmonic (0.546 Hz), are K1 = 19.7 and K2 = 4.4. Other investigators have used a single-frequency analysis to determine these parameters but have done so with a highly error-prone technique, that of examining zero crossings. Similarly zero crossings can be directly applied to Eq. (5) to determine K,, K,. However, the use of zero crossings in a noisy signal can lead to erroneous results. Even the phase shift inherent in the filter can distort the time relationships between pressure and volume signals. This time distortion is also present in linear phase filters since pressure and volume signals do not have identical frequency spectrums.
ESTIMATION
OF RESPIRATORY
31
PARAMETERS
Parameter Determination by Variance Ratios Assume that the differential equation governing the relationship and volume is given by Eq. (1) :
between pressure
P=K,V+KJ?
(1)
Assuming that the dc terms have been removed, the covariance between pressure and volume is o&, = E{PV} = E{K, Vz + K2 liV}
(13)
= Kl E{ V’} + K2 E{ J+‘) = Kl o;, + K2 L&. However, the covariance between a real band-limited identically zero (see Appendix B). Consequently, K 1 d&&
(14) signal and its derivative is
covariance between pressure and volume signal mean square value of volume signal
(15a)
covariance between pressure and flow signal mean square value of flow signal
(15b)
In a similar manner,
,g&=
o;p
There are several advantages in using covariance ratios to estimate respiratory parameters. First, since both Kl and K2 are determined by dividing a covariance by a mean-squared value (which is the total signal power), the calculations are inherently normalized to signal strength. Second, filtering is not required. The use of mean values provides a measure of noise immunity provided the noise is uncorrelated with the signal. Third, zero crossings are not required. As stated earlier, zero-crossing techniques are at best treacherous and should be avoided if possible. Fourth, the calculations are simple. Although the method of single-frequency analysis of the spectral component with maximum energy and the method of leastsquares analysis will provide reliable numbers, there is a certain minimum amount of calculation which must be performed for each method. The single-frequency analysis requires a Fourier decomposition followed by maximum detection and then simple ratios. The least-squares procedure involves data storage, matrix multiplication, and matrix inversion. With judicious programming, these techniques can produce numbers quite rapidly; however, a digital computation facility is required. Alternatively, although digital computation is desirable for the covariance ratio method, it is not necessary. The calculations required can be performed by analog computation, providing an on-line running estimate of the parameters K,, K2. To illustrate this point, examine the block diagram shown in Fig. 4. The signals P and V are passed through an R-C circuit which removes the mean .T
32
ROY,
POWERS,
JR., AND
KIMBALL
or dc value from the signals. The signal pihas zero mean. Alternatively, an integrator and a subtractor could be used to remove the mean from each of the signals. The zero mean signals V’ and P are then squared and cross-multiplied by P to form V2, PV, Pp, and (lj)l. Each of these signals are then passed through “leaky integrators”. Basically, these are operational amplifier integrators with a time constant adjusted so that past data can “leak off” exponentially and new data incorporated into the integrator. The output of the integrators provides a running measure of the integral of the input without saturation of the amplifier. The formation of K,, K2 requires the ratios of expected value which are produced by integration with respect to time and then division by the integration time. Since C v
” 5
-II
(
2 a
I2
s+a
R
A DIVIDE
PV
MULTIPLY
K,
T
0
s+a
A
p-it P;
0
s+a
3 UIVIDE
T FIG.
I
7
4. Analog computation
vh2
a
K2
?
s+a
of respiratory parameters.
only ratios are formed, division of each integral by the time of observation is not required since the time baseappearsin both numerator and denominator and therefore cancel. The output of the division circuits are the desired numbers Kl, Kz. If K,, Kz vary with respect to time, this arrangement of exponential averaging will allow a continuous tracking of the changing values. This circuit represents a method of obtaining on-line measurementof the respiratory parameters Kl and Kz. As an example, a running estimate of parameter K1 at the end of successive breaths is listed in Table 2. Also shown is the sameestimate as provided by least squares, separated into inspiration and expiration. These estimates were obtained on the raw unfiltered data. At the conclusion of ten breaths, there is only a 3 % difference between the values obtained by covariance ratios and that obtained by least squares.This is to be expected, for as shown in Appendix A, the covariance ratio technique is actually a least-squaresprocedure simplified by the fact that the covariance of a real band-limited signal and its time derivative is identically zero. The small differences shown are due to both the separation into inspiration and expiration and a nonzero value of COV{XJ?} over a finite time interval. As shown
33
ESTIMATION OF RESPIRATORY PARAMETERS
in Table 2, this difference is small, showing that the covariance ratio technique is a viable method for obtaining respiratory parameters. TABLE CUMULATIVE
End of breath
1 2 3 4 5 6 7 8 9 10
2
ESTIMATE OF PARAMETER
Covariance
ratios
19.3 20.4 20.2 20.1 20.2 20.7 20.9 20.8 20.4 20.3
K,
Least squares Insp. Exp. 19.4 17.8 18.5 19.2 19.5 19.3 19.5 19.7 19.7 19.7
18.2 18.1 19.1 19.9 19.8 19.9 19.9 19.9 19.8 19.7
CONCLUSIONS
A new technique has been presented for the estimation of the respiratory parameters associated with a first-order linear model of respiratory mechanics. This method, that of covariance ratios, is shown to be a modification of least-squares analysis. However, this modification allows for a considerably simplified (but mathematically correct) estimation procedure which can be programmed with a small analog computer. Consequently a real-time on-line least-squaresestimation of respiratory parameters can be provided without the necessitydigital computation and data storage. APPENDIX
A: LEAST-SQUARES
ANALYSIS
(3)
Consider a systemwith a single scalar output y which is to be approximated by a linear combination of a set of linearly independent variables xl, x2,. . ., x,. Call the approximation j. Thus L’~~=k~x~+k,x,+.~.+k,x,. (A-1) Assume that A4 successivemeasurementsof both y and the variables xi are made. A set of equations could then be formed : ~‘(1)~=:(1)=k,x,(l)+k,x,(l)+~~~+k,x,(l) y(2) z j(2) = kl x,(2) + k2 x2(2) + . . . + kN x,,(2) .... ..... ..... ..... ..... .... ...... ... . ~‘(M)~~(~)=k~x,(M)+k,x,(M)+..-+kk,x,(M)
(A-2)
34
ROY,
POWERS,
JR.,
AND
KIMBALL
To simplify the notation, let each column variable be designated as a vector such that (A-3)
y~:=kklx,$-k,x,+.‘.+li,~x,.
Alternatively,
a matrix X can be formed consisting of the column vectors x=
[XlX2...Xh]
(A-4)
Y=Xk.
r A-5)
such that
The error between the true value y and the approximation vector E:
y is given by the error
g=y-f=y-Xk.
(A-6)
The problem is to determine a vector k which minimizes the error vector in a least-squares sense. A pictorial appreciation of the procedure involved in leastsquares approximation is provided in Fig. A-l. For the purpose of simplicity, only two variables, x1 and x2, will be used, thus only two vectors x1 and x2 are shown.
FIG.
A-l.
Orthogonal projection:
the principle of least-squares fitting.
These two vectors x1 and x2 form a plane. The measurement vector y is assumed to lie outside of the plane formed by x1 and x2. The approximation to y is formed by a linear combination of x1 and x2 : f=klxl+k,x,.
t A-7) This approximation must lie in the plane formed by x1 and x2. From basic geometric considerations, the best approximation to y is given by an orthogonal projection of y onto the x1, x2 plane. This is the orthogonal projection theorem, which forms the basis for least-squares approximation. If J; is formed as shown, then 9 is a leastquares approximation to y. Note that this projection requires that the error vector
ESTIMATION
OF RESPIRATORY
35
PARAMETERS
E is orthogonal to the plane x1, x2 and, thus, is orthogonal to both x1 and x2. This concept is easily extended to higher dimensions. Furthermore, note that perfect matching of y is possible only when the measurement y lies within the same space subtended by the vectors xi. Therefore, it is extremely important that the independent variables x1, x,, . . ., X~ form as complete a set of measurements as possible and that a sufficient number of measurements be made (M > N). The presence of measurement noise will of course add an additional dimension to the space such that y cannot possible lie in the subspace subtended by the xi vectors. Continuing with the analysis, observe that E must be simultaneously orthogonal to all of the xi vectors. This requirement can be written as XiT&=O
for all i,
or XTE = 0.
(A-8)
From Eq. (A-6), XTz=O=XTy-X=Xk or
(A-9)
k = (X’ X)-l XT y.
Equation (A-9) is the basic equation of least-squares analysis. It is instructive to note the dimensions of each of the terms of this equation : k=Nx
g=Mx
N variables, M successive measurements,
1, 1,
X=MxN, (XT X)-l = N x N.
For the case illustrated, where N = 2, the vector k is of order (2 x 1). Since the inverse of a 2 x 2 matrix is given by
A-’ =
(A-10)
the vector k is given by igl
x22@?
C ~~(9~49- 2 x1(9
& x12W C x2W.W
x2(i)
- C x1(9x2(i) C xlG).W
C xl”(i) 2 x2’(i) - CCxl(~)x2W”.
It is important
I-
C x2G>v(i)
(A-l 1)
to note that if the cross product term is zero,
jtl xl(i) x2(i)= 0,
(A-12)
ROY,POWERS,JR.,AND
36
then the values fork,
KlMBALL
and k, are simply
(A-13a)
M 2 x,(i)y(i k, = k&e-~-
j
‘;; x,2(i) ,Tl
Equation (A-l 3) will be used to show that the covariance ratio technique is identical to a least-squares procedure, provided that the covariance between x, and .x2 is identically zero. APPENDIX The cross correlation
B: COVARIANCE RATIO ANALYSIS (3)
function between a signal and its derivative is given by
(B-1)
R&l, t2)= Wt,) qaj. where E{ } represents the expectation operator. In terms of the definition derivative, [B-l] can be rewritten as
x(t, + 6)- x(fJ -.---irm--L--
of a
I)
(B-2)
If x(t) is a stationary signal, the correlation function is dependent only on the difference between tl and t2 and not on their absolute values. Letting t1 - f2 = z. the correlation function can be stated as R,*(z) = -aR,,(r)/aT.
Let x(t) be a real band-limited series
(B-3)
signal which can be represented by a finite Fourier
Ckejkwot.
(B-4)
ESTIMATION OF RESPIRATORY PARAMETERS
37
The correlation function R,,(z) is given by R,,(T) = E
5 Ck &k%t + C, ejj.-o(t+r) AYN ck&N 1
=E
5 2 ck Ciei(k+E.)wotejl.wor . (B-5) 1,,-N k=-N 1 Due to the orthogonality of the exponential functions, the above expectation is zero except when k + 1. = 0. Therefore, R,,(T) = 5 C, Cmk(ejkwor + e-jkoor). k=-N
(B-6)
Since x(t) is real, the Fourier coefficients are symmetric about k = 0, thus Ck = C-,. R,,(z) = C,” f 2 : C,” cos kw,. k=l
(B-7):
From [B-3],
R,JT) =-v
=2 ? Ck2(kCo,sin ko, 2). k:l
(B-8)
The zero value of R,,(z) is identically zero: R,,(O) = 0.
(B-9)
The zero value of a correlation function represents a mean-square value: R,,(t) = &x(t) x(t + T)}.
(B-10)
R,,(O) = E{x2(t)) = x’(t).
(B-l I)
For T = 0, For a signal whose mean value is zero, x2(t) = mean-square value = oh = variance of x.
(B-12)
In summary, a:, = c,’ + 2 ;: ck2, kT1
CT:, = 0.
(B-13)
ACKNOWLEDGMENTS This work was performed in conjunction with the Departments of Surgery and Physiology, Albany Medical College of Union University, Albany, NY and the Systems Engineering Division of Rensselaer Polytechnic Institute, Troy, NY. This work was supported by Grant No. GM 15426 from the National Institute of General
38
ROY,
POWERS,
JR.,
AND
KIMBALL
Medical Sciences, General Clinical Research Center Grant No. 5 MO I-RR 00094m, National Institute of General Medical Sciences Special Fellowship No. 5 F03 GM 53156-02, and Project Themis Contract No. DAAB07-69-C-0365.
REFERENCES 1. ALBRIGHT, C. D., AND BONDURANT, S. Some effects of respiratory frequency on pulmonary mechanics. J. Clin. Invest. 44(S), 1362-1370 (1965). 2. ASHBAUGH, D. G. Continuous positive pressure breathing (CPPB) in adult respiratory distress syndrome. J. Thoruc. Curdio. Surg. 57,31-41 (1969). 3. DERUSSO, P. M., ROY, R. J., AND CLOSE, C. M. “State Variables for Engineers.” John Wiley and Sons, Inc., New York, 1965. 4. FALKE, K. .I. Ventilation with end-expiratory pressure in acute lung disease. J. Clin. Invest. 51(9),
2315-2323
(1972).
5. FEINBERG, B. N. Parameter estimation: A diagnostic aid for lung diseases. Instrumentation Technol. 4Cb46 (1970). 6. FENN, W. O., AND H. RAHN (Ed.), Handbook of Physiology, Section 3 : Respiration, Vol. 1, pp. 363-376,411428,451462. Washington, D.C.: American Physiological Society, 1964. 7. FLETCHER,G., AND BELLVILLE, J. W. On line computation of pulmonary compliance and work of breathing. J. A&. Physiol. 21(4), 1321-1327 (1966). 8. FRY, D. L., AND HYATT, R. E. A unified analysis of the relationship between pressure, volume and gasflow in the lungs of normal and diseased human subjects. Amer. J. Med. 672-689 (1960). 9. GISSER,D. G., NEWELL, J. C., AND POWERS,S. R., JR. A refined bag-box spirometer. J. .4ss. Advan. Med. Znstr. 6, 167 (1972). 10. JODAT, R. W. Simulation of respiratory mechanics. BiophE’s. J. 6,773-785 (1966). 11. LEWIS, F. J. Analysis of respiration by an on-line digital computer system. Proc. Amer. Surg. Ass. 211-221 (1966). 12. MCILROY, M. B. A new method for measurement of compliance and resistance of lungs and thorax. J. Appl. Physiol. 18(2), 424-427 (1963). 13. MCWILLIAMS, R., AND ADAMS, A. H. Analogue computer for study of breathing mechanics. Med. Electron. Biol. Eng. 1,353 (1963). 14. MEAD, J., AND WH~~TENBERGER,J. L. Physical properties of human lungs measured during spontaneous respiration. J. Appl. Physiol. 5,779 (1953). 15. MEAD, J. Mechanical properties of lungs. Physiol. Rev. 41,281 (1964). 16. MEDITCH, J. S., AND STOLL, P. J. On respiratory system parameter estimation. 1970 Joint Automatic Control Conference, Paper 16-D, pp. 383-390. 17. MILLER, C. W. “A Least Squares Fit Method of Measuring Lung and Thorax Compliance in a Continuous Patient Monitoring Systems.” 22ACEMB, Paper 10.13, p. 124,197O. 18. MINDLIN, S. “On-Line Computation of Respiratory Characteristics.” EAI Applications Bulletin 953055, 1969. 19. MONACO, V., BURDGE,R., NEWELL, J., LEATHER, R., AND POWERS,S. R., JK. Clinical significance of the functional residual capacity in the post-injury state. Sarg. Forum 22,42-43 (1971). 20. OSTRANDER,L., AND CHESTER,E. H. “Frequency Analysis of Pulmonary Pressure-Flow Data.’ 23ACEMB, Paper 16.11, p. 208,197O. 21. PETERS,R. M. Work of breathing. J. Trauma 8(5), 915-921 (1968). 22. PONTOPPIDAN, M., GEFFIN, B., AND LOWENSTEIN, E. Acute respiratory failure in the adult. (Three Parts) New Eng. J. Med. 287: 14,15,16,690-698,743-752,799-806,1972. 23. PRICE, A., AND POWERS, S. R., JR. “The Effect of a Period of Acute Hemorrhagic Shock in the Respiratory Compliance of the Dog,” 21st ACEMB, Paper 42.1, 1968.
ESTIMATION
OF RESPIRATORY
PARAMETERS
39
24. ROSSING, R. G. “A Digital Computer Program for the Estimation of Respiratory Mechanical Parameters,” Proceedings of Symposium on Biomedical Engineering, Marquette University, pp. 134-136,1966. 25. SEALY, W. C. The lung in hemorrhagic shock. J. Trauma 8,774-781(1968). 26. STACY, R. W., AND PETERS, R. M. Computations of respiratory mechanical parameters. “Computers in Biomedical Research.” Chapter 11, pp. 269-288, Academic Press, New York, 1965. 27. WALD, A. A computer system for respiratory parameters. Comp. Biomed. Res. 2,41 l-429,1969.