I262
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. 12. DECEMBER 1991
Communications Deconvolution of Tracer and Dilution Data Using the Wiener Filter J . H . T . Bates Abstract-In the study of living systems it is often necessary to inject or infuse a substance into the peripheral circulation and monitor its subsequent concentration in the plasma with time. Examples abound in the pharmacokinetic study of drugs and in the use of the indicator dilution technique for measuring blood flow. Furthermore, it is often necessary to deconvolve one such measured, and hence noisy, data set with another. One of the standard methods for deconvolving noisy signals is the Wiener filter, which is generally derived as a real window in the frequency domain such that the mean squared error between the estimated deconvolved function and the truth, on average, is minimized. Application of the Wiener filter requires some (often crude) model of the noise-to-signal power ratio as a function of frequency. In the pharmacokinetic and indicator dilution situations, however, one invariably has a good model of the actual function to be deconvolved in the form of a sum of decaying exponential functions. Such a model maj be employed to calculate the signal-to-noise power ratio for use in the Wiener filter, or alternatively may be directly deconvolved itself. It is shown that better results are achieved with the Wiener filter if the model of the signal is not particularly accurate, whereas with a very accurate model it is better to deconvolve the model itself. The point at which the two deconvolution approaches perform comparably occurs when the error in the model is of a similar magnitude to the noise.
INTRODUCTION In the pharmacokinetic study of living systems it is common to measure the clearance of a tracer substance from the peripheral circulation following a bolus injection. The time-course of this clearance, h(t), may be regarded as the impulse response of the system to the tracer. Now let the same tracer be deposited in the circulation continuously from some other source by an amount given by the function x ( t ) . The time-course of tracer clearance from the circulation can again be measured, this time of by y(t). Because the system behaves linearly with respect to tracers the three functions h, x , and y are related through the convolution integral [ 11 e1
y(t) = J 0 x ( 7 ) h ( t -
7)
dr.
The convolution integral also arises in indicator dilution measurements of blood flow, where y ( t ) is the concentration of indicator at a downstream measurement site, x ( t ) is the concentration at the upstream site, and h(t) is the impulse response of the vascular bed in between [2]. While the functions y ( t ) and h(t) are readily measured experimentally, it is often the case that x ( t ) is not. For example, x ( t ) may represent the conversion of one substance to another or the release of tracer from an organ such as the liver. One therefore needs to be able to invert the integral in (1) to obtain h(t) given x ( t ) and y(r)-a process known as deconvolution. This would be straight-
forward if the measurements of x ( t ) and y ( t ) were completely free of noise, which of course is never actually the case. In practice, the process of deconvolution tends to greatly amplify measurement noise. One of the classical techniques for deconvolving a noisy signal is the Wiener filter, which seeks to find the solution that is closest to the truth in the least squares sense [3]-[6]. The Wiener filter was originally derived for giving the best reconstruction, on average, of an ensemble of signals o r images, each with the same deteministic component but with different noise [4], [5]. In the pharmacokinetic situation we consider here, however, one has only a single y(t) to deconvolve, and is interested in doing so with the greatest accuracy in that one particular case. Furthermore, precise application of the Wiener filter is actually impossible because it requires complete knowledge of the signal-to-noise power ratio as a function of frequency. In practice, therefore, one must employ an approximation to the ideal filter using some model of the signal-tonoise power ratio 161. The main purpose of the study reported in this communication was to examine the application of the Wiener filter to the deconvolution of tracer-type signals (by which is meant any smooth concentration signals of biological origin). This led to a rather more fundamental investigation of the Wiener filter itself, because of the fact that tracer-type signals are invariably well described by parametric curves consisting of polynomials o r sums of exponential functions [2], [7], [8]. Such curves can be used to calculate the signal-to-noise power spectra for use in the Wiener filter. O n the other hand, if the parametric curves are perfect representations of the noise-free signals themselves, then obviously the deconvolution can be performed on the curves, making application of the Wiener filter superfluous.
THEORY The derivation of the Wiener filter, as commonly applied, is well known [4]-[6]. It is necessary, however, to reiterate here its derivation from first principles since its definition will be extended to include three alternative forms, one of which is the standard one. Let y ( t ) be the output of a linear system with input ~ ( tand ) impulse response h(t). Let n(t) be additive noise at the output giving rise to the measured output signal z(r) where z(t) = y(t) = x(r)
+ n(t)
(2)
and 0 denotes the operation of convolution. The problem is to find an estimate of ~ ( tgiven ) z ( t ) and h(t). This is most conveniently addressed in the Fourier domain where (2) becomes, by the convolution theorem for Fourier transforms [9], Z(w) = Y(w) =
Manuscript received August 23, 1990; revised June 4, 1991. The author is with Meakins-Christie Laboratories and Biomedical Engineering Unit, McGill University, Montreal, P.Q. H2X 2P2, Canada. IEEE Log Number 9103959.
+ n(t) 0 h(t)
+ N(w)
X(w)H(w)
+ N(w)
(3)
where capital letters denote the Fourier transforms of corresponding small letters and w is angular frequency. We can now define “simple-minded” deconvolution to be the
0018-9294/91$01.00 0 1991 IEEE
1
I
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
straightforward division of Z(w) by H(w) to give an estimate of X(w) equal to
d(w)
d(w) = Z(w)/H(w) =
x ( ~+)N ( ~ ) / H ( ~ )
sr
Id(w) - X(w)12dw
-m
(4)
(6)
and W(w) is the Wiener window whose function is to suppress those frequencies dominated by the noise. (The Wiener filter is simply the Wiener window divided by H(w).) Note that minimizing J is equivalent to requiring that the time-domain estimate T(t) of x ( t ) be optimal in the least squares sense, as a result of Rayleigh’s theorem for Fourier transforms [9]. Since the integrand in (5) is positive definite it suffices to find the W(w) that minimizes the integrand at every U . The integrand is
I
=
[W(Y + N ) / H - X]*[W(Y + N ) / H - XI
(7)
where * denotes the complex conjugate of the preceding quantity and the explicit w-dependence has been dropped for clarity. The most general solution for W is found by letting W be a complex function of w.Differentiating I in (7) with respect to W , setting the result equal to zero and solving for W then yields
W, = Y / ( Y
+N)
+
+
IYI2 ( Y ) *+ I N ( 2 ’
(10)
W , is the least constrained of the three windows. Substitution of
(8) into ( 6 ) shows that it yields a perfect reconstruction of the desired signal. That is,
X, = Z W , / H =
x.
(1 1)
This is not a particularly useful result for the practicing deconvolver, however, since it requires precise knowledge of N , with which there would of course be no problem. The implementation of W , thus requires estimates of models of Y and N of N . Having one, say, automatically provides the other since
N=z-P.
I
N=N-D
(14)
XI = Z(Y =X
+ D)/[H(Y+ N)]
+ D/H.
(15)
Note that this is precisely the same result as would be obtained by bypassing any consideration of the Wiener window altogether, and simply performing “simple-minded” deconvolution [see (4)] using P instead of Z. The error in the reconstructed signal d , is
E , = Xj - X =
D/H.
(16)
Now consider the reconstruction of X provided by W 2 , with and N substituted for Y and N , respectively, (9) to give
W2 = -
-
1912 + (f*N
P
+ m*)/2
( P + NI2
+ [ ( Y + D)*(N - D ) + ( Y + D ) ( N - D ) * ] / 2 IY + N I 2 [(Y + N)*(N - 0) + ( Y + N ) ( N - D ) * ] / 2 IY + NI2
J Y + DI2
= I -
where R{} denotes the real part of the quantity enclosed in the braces, so that, from ( 6 ) ,
d2 = Z W 2 / H
= x +R { ( Y + N ) D * } + 5 { ( Y + N ) N * } H(Y
(9)
This is the second form of the Wiener window in its extended definition. If we further assume that the noise is uncorrelated with itself and with the signal then, on average, the noise cross terms in (9) are zero and we have the third (and usual) form of the Wiener window thus
W, =
(13)
(8)
which is the first form of the Wiener window in its extended definition. Now we constrain W in (7) to be a real function of w and solve in the same manner to get
1YI2 + (Y*N YN*)/2 W, = IY N I 2
F=Y+D
where D is the (unknown and hopefully small) error in the model of Y. Substitution of f and N for Y and N , respectively, in (8), followed by substitution of the result into ( 6 ) , yields
(5)
where
x(w) = Z(w) W(w)/H(w)
Let
so that
which is often an extremely poor reconstruction due to the amplification of noise at those frequencies where H(w) is negligible when N(w) is not. The aim of the Wiener filter is to find the 8 ( w ) that minimizes the cost function
J =
1263
(12)
+ N)*
(18)
where 5 { } denotes the imaginary part of the quantity enclosed in the braces. The error in d2 is thus
E2 = X 2 -
-
-
R{(Y
X
+ N)*D} + 5 { ( Y + N)*N} H(Y + N)*
An interesting question now arises: “can
(19)
Z2 ever be better than
d, (in the least squares sense) when they are both calculated using the same modes P a n d hf.Intuition might suggest not, in view of the fact that W, is more constrained (i.e., to be real) than W , . It turns out, however, that X2 can indeed be better than d ,. The condition when this occurs is found by expressing E , [see (16)] in a similar form to E2 in (19); thus
E,
=
D/H
+ N ) * / [ H ( Y + N)*] - R { ( Y + N)*D} + 5 { ( Y + N)*D} H(Y + N)* =
D(Y
(20)
Now, X2 is a better estimate than 8,at a given value of w if (E2I2 is smaller than IEll2. By Rayleigh’s theorem [9], T2(t) is a better
1264
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
estimate of x ( t ) than is f , ( t ) in the least-squares sense if the integral of lE2I2 over all w is less than that of (Ell2.Inspection of (19) and (20) reveals that this occurs when
and vice versa. This result is of practical interest because it suggests that one should deconvolve using W ,if the models F and I? are sufficiently accurate, in which case D on the right-hand side of inequality (21) is small compared to N in the corresponding position on the left-hand side. (Note that the other quantities on both sides are identical.) Indeed, in such a case the Wiener filter is superfluous because using W ,is equivalent to “simple-minded” deconvolution using the model F. Conversely, if the models P and N are poor representatives of Y and N , respectively, then D is large compared to N and one should deconvolve using W2. In general, therefore, we can say that and & are comparable estimates of X when the error in the model (D)is comparable in magnitude to the noise ( N ) . It is difficult to be more specific than this, however, because it only takes D to be negligible in a very small range of w where the denominator of above integral is nonnegligible to make the integral itself huge, with a correspondingly noisy reconstruction of x(t). Now consider the reconstruction of X using W3, again with the models f and N . From (10) we have that
which gives
X3 = Z W 3 / H
= x + y H 1
1
+ IN
-
DI2/IY+ DI2
The error in X3 is then
E3
= X3 -
x
X 1
(1
+ ( N - D I 2 / ( ( Y+ DI2 1
We can now ask a similar question as above: “can .f3 ever be better than either or both of 2, and X 2 , and if so then under what condiis better than both X, and tions?” We can see immediately that X2 in the limit as D tends to infinity (i.e., the models are infinitely bad) because E , [see (16)] and E2 [see (19)] both tend to infinity while E3 tends to the finite limit of
x3
(N/H
-
x)/2.
On the other hand, when the models are perfect, and D is zero, E, has the value
which clearly has larger magnitude than the corresponding value of E , which is zero. It is also bigger in magnitude than E2 since X2 is the best representation of X that can be obtained when W is a real function of w . Therefore, as D increases, there must be a point at which X3 becomes a better estimate of X than does X2. To summarize the above theoretical results, we have shown that the three Wiener windows W,,W2, and W,perform relatively differently depending on the accuracies of the models of the signal and noise used to implement them. In the next section this is illustrated in a numerical example.
COMPUTATIONAL RESULTS The performance of the three Wiener windows derived above were tested on simulated data having the characteristics of tracertype data. In this example noise was added only to y ( t ) since the application of the Wiener filter makes the assumptions that there is no noise in the input to the convolution integral [see (l)]. While this is by no means always the case, it does occur when ~ ( t is) a system impulse response that has been well characterized by averaging data obtained from an experiment performed many times, while y ( t ) is a noisy data set from a single experiment. Fig. l(a) shows an example single-exponential x ( t ) (dashed line) together with a single-exponential h(t) (solid line). The convolution y ( t ) of these two signals is shown in Fig. l(b) as the dashed line together with the addition of noise (Gaussian, zero mean, standard deviation 4 % of signal maximum) to produce z ( t ) (solid line). Fig. l(c) shows the result of deconvolution z(t) with h(t) using “simple-minded” deconvolution [see (4)]. The result is clearly disastrous. Fig. 2 shows the results of deconvolving z(t) with h(t) using the three Wiener windows described above, all implemented with perfect models of the noise, that is, with P = Y and N = N . Fig. l(a) shows y ( t ) and J ( t ) together; they are obviously superimposable. Fig. 2(b) shows .a?,(t) provided by W , [the inverse Fourier transform of (1 I)], which is a perfect reconstruction of x ( t ) as predicted when the models are perfect. The reconstructions provided by W2and W3 are shown (together with x ( t ) itself as a dashed line) in Fig. 2(c) and (d), respectively. The sum of the squared deviations between x ( t ) and its reconstruction in each case are given in Table I. The usual form of the Wiener window, W3, gives the poorest reconstruction in this case, while the most general form of the Window ( W , ) gives the best one (it is, in fact, perfect). As pointed out in the “Theory,” in practice one never has perfect models of Y and N with which to implement the Wiener windows. A good approximation can be obtained in the present example, however, by fitting a biexponential function to z(t). The result is compared in Fig. 3(a) (solid line) to y ( t ) (dashed line). The restructions provided by the three Wiener windows implemented with this model are shown in Fig. 3(b)-(d). That provided by W ,[Fig. 3(b)] is still extremely good, while the worst result is again provided by W, (Table I). This bears out the theoretical development above, in which it was noted that an extremely good model of Y [or, equivalently, y ( t ) ] yields the best reconstructions of x ( t ) with W,. Finally, in Fig. 4 the results of implementing the Wiener windows with a poorer model y ( t ) are shown. The model itself is a sixth-order polynomial function fitted to z(t) and is shown in Fig. 4(a) (solid line) together with y ( t ) (dashed line). The reconstruction provided by W , is shown in Fig. 4(b) and is seen to have significant systematic deviations from x(t) that reflect those in the model. The reconstructions provided by W, and W3 are also poorer than with the previous more accurate models (Figs. 2 and 3). However, the
I265
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
-- y It1 - 5 (tl
Fig. 1. (a) The true deconvolved function x ( t ) and the “point spread function” h(t). (b) The convolved function y ( t ) and the noisy convolved function z ( f ) .(c) The result of deconvolving z ( t ) with h ( f ) using the “simpleminded” method.
b L
i i , (tl
\
\
Fig. 3 . (a) The noiseless convolved function y ( t ) (dashed line) and the biexponental model y(r) (solid line). (b) The results of deconvolving the noisy convolved function z ( f ) using the Wiener window W , with the biexponential model. (c) The result of deconvolving z(t) using W,. (d) The result of deconvolving z(r) using W,. In (b)-(d) the true deconvolved function x ( f ) is shown as a dashed line for comparison.
L y It1 - It)
--
Y
Fig. 2. (a) The noiseless convolved function y(r) and the perfect model y ( t ) (the two functions are identical). (b) The results of deconvolving the noisy convolved function z(t) using the Wiener window W , with the perfect model. (c) The result of deconvolving z ( t ) using W,. (d) The result of deconvolving z ( t ) using W,. In (b)-(d) the true deconvolved function X ( t ) is shown as a dashed line for comparison.
i 1 (tl
TABLE I RESULTS OF DECONVOLVING THE NOISY SIGNAL SHOWN I N FIG. l(b) WITH THE THREE WIENER WINDOWS W , , W,, A N D W,. THRRE MODELS OF THE NOISYSIGNAL WERE USED.THERESULTS OF THE PERFECT MODELARE SHOWN IN THE FIRSTCOLUMN OF THE TABLEA N D I N FIG. 2. THERESULTS OF THE BIEXPONENTIAL MODELARE SHOWN I N THE SECOND COLUMN AND IN FIG. 3. THERESULTS OF THE SIXTH-ORDER POLYNOMIAL MODEL ARE SHOWN I N THE THIRD COLUMN AND IN FIG.4 Wiener Window Wl
W, W3
Perfect Model
Biexponential Model
Sixth-Order Polynomial Model
0 0.395 0.656
0.010 0.450 0.648
1.637 1.152 1.102
Fig. 4. (a) The noiseless convolved function y ( f ) (dashed line) and the sixth-order polynomial model y(r) (solid line). (b) The results of deconvolving the noisy convolved function z ( f ) using the Wiener window W , with the sixth-order polynomial model. (c) The result of deconvolving z(f) using W,. (d) The result of deconvolving z ( f )using W,. In (b)-(d) the true deconvolved function x ( f ) is shown as a dashed line for comparison.
DISCUSSION A N D CONCLUSIONS order of performance of the windows has now been reversed, with W, performing the best and W , the worst (Table I). This again exemplifies the theoretical result that a sufficiently poor model of Y requires the use of W, for best results.
1-
The advantage of the Wiener filter over some other methods proposed for deconvolving pharmacokinetic or indicator diluation data [7], [8], [lo] is that it makes no a priori assumption about the functional form of the data except in so far as is required for an
1266
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12. DECEMBER 1991
estimate of the Wiener window. However, the typical form of such data means that good empirical parametric models are invariably available, so that accurate estimates of the Wiener window can be made. Of course, one might expect that if the models are sufficiently accurate then the deconvolution can be performed on the models themselves, thereby bypassing any consideration of the Wiener filter. In this communication it was shown that such a crossover point in model accuracy does exist, one side of which the best results are obtained with the standard Wiener filter ( W 3 )and the other side of which it is better to perform a straightforward deconvolution of the models themselves using, for example, the method of serial division [9]. So far, it has not been demonstrated, either theoretically or by example, that W, can ever perform better than both W ,and W3. Indeed, experience seems to indicate that, while W, may certainly perform better than W , , under those conditions when it does then W3seems to perform better than either. The choice of Wiener window is thus reduced to that of W , (which is tantamount to using the model y(t) itself with no window at all) or W3 (the usual form). Although it appears difficult to make specific statements about when to use either, the results of this study indicate that the availability of a good parameteric model for y ( t ) (such as the biexponential model of the above example) would allow the use of W ,. The definition of a good model here means one whose residuals (the differences between the model and the signal being modeled) are noiselike in character. A poorer model, whose residuals would exhibit systematic deviations about zero (indicating insufficient degrees of freedom or inappropriate structure, such as the parametric but highly inappropriate polynomial model in Fig. 3) would require the use of W,. In conclusion, then, it is possible to suggest the following procedure for deciding on how to use the Wiener filter in the deconvolution of tracer-type data. First, fit an appropriate parametric model to the data to be deconvolved. Next, examine the residuals between the data and the fitted model. If they are essentially noiselike in character then one should deconvolve the model itself. If the residuals have clear deterministic features (i.e., systematic deviations about zero) then one should use the model to implement the standard Wiener filter ( W 3 ) . Finally, although this study has focused on the deconvolution of biological data, its results apply equally well to any data whose deterministic part is smooth and well representable by a parametric model. ACKNOWLEDGMENT The helpful comments of Dr. R. Keamey and Dr. I. Hunter are gratefully acknowledged. REFERENCES [I] D. M. Foster, D. G. Covell and M. Berman, “Applications of a general method for deconvolution using compartmental analysis,” Compur. B i d . Med., vol. 18, pp. 253-266, 1988.
[2] J. Gamel, W. F. Rousseau, C. R. Katholi, and E. Mesel, “Pitfalls in digital computation of the impulse response of vascular beds from indicator-dilution curves,” Circ. Res., vol. XXXII, pp. 516-523, 1973. [3] N. Wiener, Extrapolation, Interpolation and Smoothing of Stationary Time Series. New York: Wiley, 1949. [4] R. N. Bracewell, “Restoration in the presence of errors,” Proc. IRE, vol. 46, pp. 106-111, 1958. [5] J. W. Brault and 0. R. White, “The analysis and restoration of astronomical data via the fast Fourier transform,” Asrron. Astrophys., vol. 13, pp. 169-189, 1971.
W. H. Press, B. P. Flannery, S. A . Teukolsky, and W. T. Vetterlin, London: Cambridge Univ. Press, 1986, pp. 4 17-4 19. D.J. Cutler, “Numerical deconvolution by least squares: Use of prescribed input functions,” J . Pharm. Biopharm., vol. 6, pp. 227-241, 1978. -, “Numerical deconvolution by least squares: Use of polynomials to represent the input function,” J . Pharm. Biopharm., vol. 6, pp. 243-263, 1978. R. N . Bracewell, The Fourier Transform and its Applications, 2nd ed. New York: McGraw-Hill, 1978. P. V. Pedersen, “Novel deconvolution method for linear pharmacokinetic systems with polyexponential impulse response,” J . Pharm. Sci., vol. 69, pp. 312-318, 1980. Numerical Recipes.
Numerical Implementation of Sealed-End Boundary Conditions in Cable Theory Emst Niebur and Dagmar Niebur
Abstract-We show that a frequently used numerical implementation of von Neumann boundary conditions (zero inflowing current) in cable theory is incorrect. Correct implementations are given and it is shown that they yield results in good agreement with known analytical solutions.
INTRODUCTION Cable theory is the standard tool for modeling voltage distributions in spatially extended neurons or other excitable cells, like cells found in cardiac tissue. The theory is based on the observation that the intracellular electrical potential vanes much more along a long nerve fiber than between points inside the fiber in a plane perpendicular to the fiber axis. This facilitates greatly the mathematical analysis since the spatial dimension of the differential equations for the intracellular voltage is reduced from three to one. The advent of digital computers made the numerical solution of these equations possible in cases in which no analytical solutions have been found. Any method of solution, analytical or numerical, must take into account the initial conditions and the boundary conditions which make the solution of the differential equations unique. We found several articles in the literature [ l ] , [ 2 ] , [5], [6], [9], [ l l ] , [13], [16] in which the boundary conditions were incorrectly implemented, which led to a wrong solution. This communication deals with the correct implementation of the boundary conditions. Manuscript received June 10, 1991; revised August 30, 1991. The work of E. Niebur was supported by the Swiss National Science Foundation Grants No. 2000-5.295 and 8220-25941, the Air Force Office of Scientific Research, and the James S . McDonnell Foundation. It was started at the Swiss Federal Institute of Technology (EPFL), supported by EPFL, and was completed at the Jet Propulsion Laboratory, California Institute of Technology, supported by the U.S. Department of Energy through an agreement with the National Aeronautics and Space Administration. E. Niebur is with the Institute of Theoretical Physics, University of Lausanne, CH-1015 Lausanne, Switzerland, and Computation and Neural Systems 216-76, California Institute of Technology, Pasadena, CA 91125. D. Niebur is with the Department of Electrical Engineering, Swiss Federal Institute of Technology, CH-1015 Lausanne, Switzerland and Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91 109. IEEE Log Number 9103960.
0018-9294/91$01.00 O 1991 IEEE
- 1
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. 12. DECEMBER 1991
Communications Deconvolution of Tracer and Dilution Data Using the Wiener Filter J . H . T . Bates Abstract-In the study of living systems it is often necessary to inject or infuse a substance into the peripheral circulation and monitor its subsequent concentration in the plasma with time. Examples abound in the pharmacokinetic study of drugs and in the use of the indicator dilution technique for measuring blood flow. Furthermore, it is often necessary to deconvolve one such measured, and hence noisy, data set with another. One of the standard methods for deconvolving noisy signals is the Wiener filter, which is generally derived as a real window in the frequency domain such that the mean squared error between the estimated deconvolved function and the truth, on average, is minimized. Application of the Wiener filter requires some (often crude) model of the noise-to-signal power ratio as a function of frequency. In the pharmacokinetic and indicator dilution situations, however, one invariably has a good model of the actual function to be deconvolved in the form of a sum of decaying exponential functions. Such a model maj be employed to calculate the signal-to-noise power ratio for use in the Wiener filter, or alternatively may be directly deconvolved itself. It is shown that better results are achieved with the Wiener filter if the model of the signal is not particularly accurate, whereas with a very accurate model it is better to deconvolve the model itself. The point at which the two deconvolution approaches perform comparably occurs when the error in the model is of a similar magnitude to the noise.
INTRODUCTION In the pharmacokinetic study of living systems it is common to measure the clearance of a tracer substance from the peripheral circulation following a bolus injection. The time-course of this clearance, h(t), may be regarded as the impulse response of the system to the tracer. Now let the same tracer be deposited in the circulation continuously from some other source by an amount given by the function x ( t ) . The time-course of tracer clearance from the circulation can again be measured, this time of by y(t). Because the system behaves linearly with respect to tracers the three functions h, x , and y are related through the convolution integral [ 11 e1
y(t) = J 0 x ( 7 ) h ( t -
7)
dr.
The convolution integral also arises in indicator dilution measurements of blood flow, where y ( t ) is the concentration of indicator at a downstream measurement site, x ( t ) is the concentration at the upstream site, and h(t) is the impulse response of the vascular bed in between [2]. While the functions y ( t ) and h(t) are readily measured experimentally, it is often the case that x ( t ) is not. For example, x ( t ) may represent the conversion of one substance to another or the release of tracer from an organ such as the liver. One therefore needs to be able to invert the integral in (1) to obtain h(t) given x ( t ) and y(r)-a process known as deconvolution. This would be straight-
forward if the measurements of x ( t ) and y ( t ) were completely free of noise, which of course is never actually the case. In practice, the process of deconvolution tends to greatly amplify measurement noise. One of the classical techniques for deconvolving a noisy signal is the Wiener filter, which seeks to find the solution that is closest to the truth in the least squares sense [3]-[6]. The Wiener filter was originally derived for giving the best reconstruction, on average, of an ensemble of signals o r images, each with the same deteministic component but with different noise [4], [5]. In the pharmacokinetic situation we consider here, however, one has only a single y(t) to deconvolve, and is interested in doing so with the greatest accuracy in that one particular case. Furthermore, precise application of the Wiener filter is actually impossible because it requires complete knowledge of the signal-to-noise power ratio as a function of frequency. In practice, therefore, one must employ an approximation to the ideal filter using some model of the signal-tonoise power ratio 161. The main purpose of the study reported in this communication was to examine the application of the Wiener filter to the deconvolution of tracer-type signals (by which is meant any smooth concentration signals of biological origin). This led to a rather more fundamental investigation of the Wiener filter itself, because of the fact that tracer-type signals are invariably well described by parametric curves consisting of polynomials o r sums of exponential functions [2], [7], [8]. Such curves can be used to calculate the signal-to-noise power spectra for use in the Wiener filter. O n the other hand, if the parametric curves are perfect representations of the noise-free signals themselves, then obviously the deconvolution can be performed on the curves, making application of the Wiener filter superfluous.
THEORY The derivation of the Wiener filter, as commonly applied, is well known [4]-[6]. It is necessary, however, to reiterate here its derivation from first principles since its definition will be extended to include three alternative forms, one of which is the standard one. Let y ( t ) be the output of a linear system with input ~ ( tand ) impulse response h(t). Let n(t) be additive noise at the output giving rise to the measured output signal z(r) where z(t) = y(t) = x(r)
+ n(t)
(2)
and 0 denotes the operation of convolution. The problem is to find an estimate of ~ ( tgiven ) z ( t ) and h(t). This is most conveniently addressed in the Fourier domain where (2) becomes, by the convolution theorem for Fourier transforms [9], Z(w) = Y(w) =
Manuscript received August 23, 1990; revised June 4, 1991. The author is with Meakins-Christie Laboratories and Biomedical Engineering Unit, McGill University, Montreal, P.Q. H2X 2P2, Canada. IEEE Log Number 9103959.
+ n(t) 0 h(t)
+ N(w)
X(w)H(w)
+ N(w)
(3)
where capital letters denote the Fourier transforms of corresponding small letters and w is angular frequency. We can now define “simple-minded” deconvolution to be the
0018-9294/91$01.00 0 1991 IEEE
1
I
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
straightforward division of Z(w) by H(w) to give an estimate of X(w) equal to
d(w)
d(w) = Z(w)/H(w) =
x ( ~+)N ( ~ ) / H ( ~ )
sr
Id(w) - X(w)12dw
-m
(4)
(6)
and W(w) is the Wiener window whose function is to suppress those frequencies dominated by the noise. (The Wiener filter is simply the Wiener window divided by H(w).) Note that minimizing J is equivalent to requiring that the time-domain estimate T(t) of x ( t ) be optimal in the least squares sense, as a result of Rayleigh’s theorem for Fourier transforms [9]. Since the integrand in (5) is positive definite it suffices to find the W(w) that minimizes the integrand at every U . The integrand is
I
=
[W(Y + N ) / H - X]*[W(Y + N ) / H - XI
(7)
where * denotes the complex conjugate of the preceding quantity and the explicit w-dependence has been dropped for clarity. The most general solution for W is found by letting W be a complex function of w.Differentiating I in (7) with respect to W , setting the result equal to zero and solving for W then yields
W, = Y / ( Y
+N)
+
+
IYI2 ( Y ) *+ I N ( 2 ’
(10)
W , is the least constrained of the three windows. Substitution of
(8) into ( 6 ) shows that it yields a perfect reconstruction of the desired signal. That is,
X, = Z W , / H =
x.
(1 1)
This is not a particularly useful result for the practicing deconvolver, however, since it requires precise knowledge of N , with which there would of course be no problem. The implementation of W , thus requires estimates of models of Y and N of N . Having one, say, automatically provides the other since
N=z-P.
I
N=N-D
(14)
XI = Z(Y =X
+ D)/[H(Y+ N)]
+ D/H.
(15)
Note that this is precisely the same result as would be obtained by bypassing any consideration of the Wiener window altogether, and simply performing “simple-minded” deconvolution [see (4)] using P instead of Z. The error in the reconstructed signal d , is
E , = Xj - X =
D/H.
(16)
Now consider the reconstruction of X provided by W 2 , with and N substituted for Y and N , respectively, (9) to give
W2 = -
-
1912 + (f*N
P
+ m*)/2
( P + NI2
+ [ ( Y + D)*(N - D ) + ( Y + D ) ( N - D ) * ] / 2 IY + N I 2 [(Y + N)*(N - 0) + ( Y + N ) ( N - D ) * ] / 2 IY + NI2
J Y + DI2
= I -
where R{} denotes the real part of the quantity enclosed in the braces, so that, from ( 6 ) ,
d2 = Z W 2 / H
= x +R { ( Y + N ) D * } + 5 { ( Y + N ) N * } H(Y
(9)
This is the second form of the Wiener window in its extended definition. If we further assume that the noise is uncorrelated with itself and with the signal then, on average, the noise cross terms in (9) are zero and we have the third (and usual) form of the Wiener window thus
W, =
(13)
(8)
which is the first form of the Wiener window in its extended definition. Now we constrain W in (7) to be a real function of w and solve in the same manner to get
1YI2 + (Y*N YN*)/2 W, = IY N I 2
F=Y+D
where D is the (unknown and hopefully small) error in the model of Y. Substitution of f and N for Y and N , respectively, in (8), followed by substitution of the result into ( 6 ) , yields
(5)
where
x(w) = Z(w) W(w)/H(w)
Let
so that
which is often an extremely poor reconstruction due to the amplification of noise at those frequencies where H(w) is negligible when N(w) is not. The aim of the Wiener filter is to find the 8 ( w ) that minimizes the cost function
J =
1263
(12)
+ N)*
(18)
where 5 { } denotes the imaginary part of the quantity enclosed in the braces. The error in d2 is thus
E2 = X 2 -
-
-
R{(Y
X
+ N)*D} + 5 { ( Y + N)*N} H(Y + N)*
An interesting question now arises: “can
(19)
Z2 ever be better than
d, (in the least squares sense) when they are both calculated using the same modes P a n d hf.Intuition might suggest not, in view of the fact that W, is more constrained (i.e., to be real) than W , . It turns out, however, that X2 can indeed be better than d ,. The condition when this occurs is found by expressing E , [see (16)] in a similar form to E2 in (19); thus
E,
=
D/H
+ N ) * / [ H ( Y + N)*] - R { ( Y + N)*D} + 5 { ( Y + N)*D} H(Y + N)* =
D(Y
(20)
Now, X2 is a better estimate than 8,at a given value of w if (E2I2 is smaller than IEll2. By Rayleigh’s theorem [9], T2(t) is a better
1264
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
estimate of x ( t ) than is f , ( t ) in the least-squares sense if the integral of lE2I2 over all w is less than that of (Ell2.Inspection of (19) and (20) reveals that this occurs when
and vice versa. This result is of practical interest because it suggests that one should deconvolve using W ,if the models F and I? are sufficiently accurate, in which case D on the right-hand side of inequality (21) is small compared to N in the corresponding position on the left-hand side. (Note that the other quantities on both sides are identical.) Indeed, in such a case the Wiener filter is superfluous because using W ,is equivalent to “simple-minded” deconvolution using the model F. Conversely, if the models P and N are poor representatives of Y and N , respectively, then D is large compared to N and one should deconvolve using W2. In general, therefore, we can say that and & are comparable estimates of X when the error in the model (D)is comparable in magnitude to the noise ( N ) . It is difficult to be more specific than this, however, because it only takes D to be negligible in a very small range of w where the denominator of above integral is nonnegligible to make the integral itself huge, with a correspondingly noisy reconstruction of x(t). Now consider the reconstruction of X using W3, again with the models f and N . From (10) we have that
which gives
X3 = Z W 3 / H
= x + y H 1
1
+ IN
-
DI2/IY+ DI2
The error in X3 is then
E3
= X3 -
x
X 1
(1
+ ( N - D I 2 / ( ( Y+ DI2 1
We can now ask a similar question as above: “can .f3 ever be better than either or both of 2, and X 2 , and if so then under what condiis better than both X, and tions?” We can see immediately that X2 in the limit as D tends to infinity (i.e., the models are infinitely bad) because E , [see (16)] and E2 [see (19)] both tend to infinity while E3 tends to the finite limit of
x3
(N/H
-
x)/2.
On the other hand, when the models are perfect, and D is zero, E, has the value
which clearly has larger magnitude than the corresponding value of E , which is zero. It is also bigger in magnitude than E2 since X2 is the best representation of X that can be obtained when W is a real function of w . Therefore, as D increases, there must be a point at which X3 becomes a better estimate of X than does X2. To summarize the above theoretical results, we have shown that the three Wiener windows W,,W2, and W,perform relatively differently depending on the accuracies of the models of the signal and noise used to implement them. In the next section this is illustrated in a numerical example.
COMPUTATIONAL RESULTS The performance of the three Wiener windows derived above were tested on simulated data having the characteristics of tracertype data. In this example noise was added only to y ( t ) since the application of the Wiener filter makes the assumptions that there is no noise in the input to the convolution integral [see (l)]. While this is by no means always the case, it does occur when ~ ( t is) a system impulse response that has been well characterized by averaging data obtained from an experiment performed many times, while y ( t ) is a noisy data set from a single experiment. Fig. l(a) shows an example single-exponential x ( t ) (dashed line) together with a single-exponential h(t) (solid line). The convolution y ( t ) of these two signals is shown in Fig. l(b) as the dashed line together with the addition of noise (Gaussian, zero mean, standard deviation 4 % of signal maximum) to produce z ( t ) (solid line). Fig. l(c) shows the result of deconvolution z(t) with h(t) using “simple-minded” deconvolution [see (4)]. The result is clearly disastrous. Fig. 2 shows the results of deconvolving z(t) with h(t) using the three Wiener windows described above, all implemented with perfect models of the noise, that is, with P = Y and N = N . Fig. l(a) shows y ( t ) and J ( t ) together; they are obviously superimposable. Fig. 2(b) shows .a?,(t) provided by W , [the inverse Fourier transform of (1 I)], which is a perfect reconstruction of x ( t ) as predicted when the models are perfect. The reconstructions provided by W2and W3 are shown (together with x ( t ) itself as a dashed line) in Fig. 2(c) and (d), respectively. The sum of the squared deviations between x ( t ) and its reconstruction in each case are given in Table I. The usual form of the Wiener window, W3, gives the poorest reconstruction in this case, while the most general form of the Window ( W , ) gives the best one (it is, in fact, perfect). As pointed out in the “Theory,” in practice one never has perfect models of Y and N with which to implement the Wiener windows. A good approximation can be obtained in the present example, however, by fitting a biexponential function to z(t). The result is compared in Fig. 3(a) (solid line) to y ( t ) (dashed line). The restructions provided by the three Wiener windows implemented with this model are shown in Fig. 3(b)-(d). That provided by W ,[Fig. 3(b)] is still extremely good, while the worst result is again provided by W, (Table I). This bears out the theoretical development above, in which it was noted that an extremely good model of Y [or, equivalently, y ( t ) ] yields the best reconstructions of x ( t ) with W,. Finally, in Fig. 4 the results of implementing the Wiener windows with a poorer model y ( t ) are shown. The model itself is a sixth-order polynomial function fitted to z(t) and is shown in Fig. 4(a) (solid line) together with y ( t ) (dashed line). The reconstruction provided by W , is shown in Fig. 4(b) and is seen to have significant systematic deviations from x(t) that reflect those in the model. The reconstructions provided by W, and W3 are also poorer than with the previous more accurate models (Figs. 2 and 3). However, the
I265
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991
-- y It1 - 5 (tl
Fig. 1. (a) The true deconvolved function x ( t ) and the “point spread function” h(t). (b) The convolved function y ( t ) and the noisy convolved function z ( f ) .(c) The result of deconvolving z ( t ) with h ( f ) using the “simpleminded” method.
b L
i i , (tl
\
\
Fig. 3 . (a) The noiseless convolved function y ( t ) (dashed line) and the biexponental model y(r) (solid line). (b) The results of deconvolving the noisy convolved function z ( f ) using the Wiener window W , with the biexponential model. (c) The result of deconvolving z(t) using W,. (d) The result of deconvolving z(r) using W,. In (b)-(d) the true deconvolved function x ( f ) is shown as a dashed line for comparison.
L y It1 - It)
--
Y
Fig. 2. (a) The noiseless convolved function y(r) and the perfect model y ( t ) (the two functions are identical). (b) The results of deconvolving the noisy convolved function z(t) using the Wiener window W , with the perfect model. (c) The result of deconvolving z ( t ) using W,. (d) The result of deconvolving z ( t ) using W,. In (b)-(d) the true deconvolved function X ( t ) is shown as a dashed line for comparison.
i 1 (tl
TABLE I RESULTS OF DECONVOLVING THE NOISY SIGNAL SHOWN I N FIG. l(b) WITH THE THREE WIENER WINDOWS W , , W,, A N D W,. THRRE MODELS OF THE NOISYSIGNAL WERE USED.THERESULTS OF THE PERFECT MODELARE SHOWN IN THE FIRSTCOLUMN OF THE TABLEA N D I N FIG. 2. THERESULTS OF THE BIEXPONENTIAL MODELARE SHOWN I N THE SECOND COLUMN AND IN FIG. 3. THERESULTS OF THE SIXTH-ORDER POLYNOMIAL MODEL ARE SHOWN I N THE THIRD COLUMN AND IN FIG.4 Wiener Window Wl
W, W3
Perfect Model
Biexponential Model
Sixth-Order Polynomial Model
0 0.395 0.656
0.010 0.450 0.648
1.637 1.152 1.102
Fig. 4. (a) The noiseless convolved function y ( f ) (dashed line) and the sixth-order polynomial model y(r) (solid line). (b) The results of deconvolving the noisy convolved function z ( f ) using the Wiener window W , with the sixth-order polynomial model. (c) The result of deconvolving z(f) using W,. (d) The result of deconvolving z ( f )using W,. In (b)-(d) the true deconvolved function x ( f ) is shown as a dashed line for comparison.
DISCUSSION A N D CONCLUSIONS order of performance of the windows has now been reversed, with W, performing the best and W , the worst (Table I). This again exemplifies the theoretical result that a sufficiently poor model of Y requires the use of W, for best results.
1-
The advantage of the Wiener filter over some other methods proposed for deconvolving pharmacokinetic or indicator diluation data [7], [8], [lo] is that it makes no a priori assumption about the functional form of the data except in so far as is required for an
1266
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12. DECEMBER 1991
estimate of the Wiener window. However, the typical form of such data means that good empirical parametric models are invariably available, so that accurate estimates of the Wiener window can be made. Of course, one might expect that if the models are sufficiently accurate then the deconvolution can be performed on the models themselves, thereby bypassing any consideration of the Wiener filter. In this communication it was shown that such a crossover point in model accuracy does exist, one side of which the best results are obtained with the standard Wiener filter ( W 3 )and the other side of which it is better to perform a straightforward deconvolution of the models themselves using, for example, the method of serial division [9]. So far, it has not been demonstrated, either theoretically or by example, that W, can ever perform better than both W ,and W3. Indeed, experience seems to indicate that, while W, may certainly perform better than W , , under those conditions when it does then W3seems to perform better than either. The choice of Wiener window is thus reduced to that of W , (which is tantamount to using the model y(t) itself with no window at all) or W3 (the usual form). Although it appears difficult to make specific statements about when to use either, the results of this study indicate that the availability of a good parameteric model for y ( t ) (such as the biexponential model of the above example) would allow the use of W ,. The definition of a good model here means one whose residuals (the differences between the model and the signal being modeled) are noiselike in character. A poorer model, whose residuals would exhibit systematic deviations about zero (indicating insufficient degrees of freedom or inappropriate structure, such as the parametric but highly inappropriate polynomial model in Fig. 3) would require the use of W,. In conclusion, then, it is possible to suggest the following procedure for deciding on how to use the Wiener filter in the deconvolution of tracer-type data. First, fit an appropriate parametric model to the data to be deconvolved. Next, examine the residuals between the data and the fitted model. If they are essentially noiselike in character then one should deconvolve the model itself. If the residuals have clear deterministic features (i.e., systematic deviations about zero) then one should use the model to implement the standard Wiener filter ( W 3 ) . Finally, although this study has focused on the deconvolution of biological data, its results apply equally well to any data whose deterministic part is smooth and well representable by a parametric model. ACKNOWLEDGMENT The helpful comments of Dr. R. Keamey and Dr. I. Hunter are gratefully acknowledged. REFERENCES [I] D. M. Foster, D. G. Covell and M. Berman, “Applications of a general method for deconvolution using compartmental analysis,” Compur. B i d . Med., vol. 18, pp. 253-266, 1988.
[2] J. Gamel, W. F. Rousseau, C. R. Katholi, and E. Mesel, “Pitfalls in digital computation of the impulse response of vascular beds from indicator-dilution curves,” Circ. Res., vol. XXXII, pp. 516-523, 1973. [3] N. Wiener, Extrapolation, Interpolation and Smoothing of Stationary Time Series. New York: Wiley, 1949. [4] R. N. Bracewell, “Restoration in the presence of errors,” Proc. IRE, vol. 46, pp. 106-111, 1958. [5] J. W. Brault and 0. R. White, “The analysis and restoration of astronomical data via the fast Fourier transform,” Asrron. Astrophys., vol. 13, pp. 169-189, 1971.
W. H. Press, B. P. Flannery, S. A . Teukolsky, and W. T. Vetterlin, London: Cambridge Univ. Press, 1986, pp. 4 17-4 19. D.J. Cutler, “Numerical deconvolution by least squares: Use of prescribed input functions,” J . Pharm. Biopharm., vol. 6, pp. 227-241, 1978. -, “Numerical deconvolution by least squares: Use of polynomials to represent the input function,” J . Pharm. Biopharm., vol. 6, pp. 243-263, 1978. R. N . Bracewell, The Fourier Transform and its Applications, 2nd ed. New York: McGraw-Hill, 1978. P. V. Pedersen, “Novel deconvolution method for linear pharmacokinetic systems with polyexponential impulse response,” J . Pharm. Sci., vol. 69, pp. 312-318, 1980. Numerical Recipes.
Numerical Implementation of Sealed-End Boundary Conditions in Cable Theory Emst Niebur and Dagmar Niebur
Abstract-We show that a frequently used numerical implementation of von Neumann boundary conditions (zero inflowing current) in cable theory is incorrect. Correct implementations are given and it is shown that they yield results in good agreement with known analytical solutions.
INTRODUCTION Cable theory is the standard tool for modeling voltage distributions in spatially extended neurons or other excitable cells, like cells found in cardiac tissue. The theory is based on the observation that the intracellular electrical potential vanes much more along a long nerve fiber than between points inside the fiber in a plane perpendicular to the fiber axis. This facilitates greatly the mathematical analysis since the spatial dimension of the differential equations for the intracellular voltage is reduced from three to one. The advent of digital computers made the numerical solution of these equations possible in cases in which no analytical solutions have been found. Any method of solution, analytical or numerical, must take into account the initial conditions and the boundary conditions which make the solution of the differential equations unique. We found several articles in the literature [ l ] , [ 2 ] , [5], [6], [9], [ l l ] , [13], [16] in which the boundary conditions were incorrectly implemented, which led to a wrong solution. This communication deals with the correct implementation of the boundary conditions. Manuscript received June 10, 1991; revised August 30, 1991. The work of E. Niebur was supported by the Swiss National Science Foundation Grants No. 2000-5.295 and 8220-25941, the Air Force Office of Scientific Research, and the James S . McDonnell Foundation. It was started at the Swiss Federal Institute of Technology (EPFL), supported by EPFL, and was completed at the Jet Propulsion Laboratory, California Institute of Technology, supported by the U.S. Department of Energy through an agreement with the National Aeronautics and Space Administration. E. Niebur is with the Institute of Theoretical Physics, University of Lausanne, CH-1015 Lausanne, Switzerland, and Computation and Neural Systems 216-76, California Institute of Technology, Pasadena, CA 91125. D. Niebur is with the Department of Electrical Engineering, Swiss Federal Institute of Technology, CH-1015 Lausanne, Switzerland and Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91 109. IEEE Log Number 9103960.
0018-9294/91$01.00 O 1991 IEEE
- 1
