- Email: [email protected]
Comment on “fit to experimental data with exponential functions using the fast fourier transform”
NUCLEAR
INSTRUMENTS
AND
METHODS
II 4 (I974) 171-173; ©
NORTH-HOLLAND
PUBLISHING
CO.
C O M M E N T O N " F I T T O E X P E R I M E N T A L DATA W I T H E X P O N E N T I A L F U N C T I O N S
USING THE FAST FOURIER TRANSFORM" l) MICHAEL R. SMITH Department of Physics, The University of Calgary and the Department of Pathology, Foothills Hospital, Calgary, Alberta, Canada T2N 1N4
and SORIN COHN-SFETCU Department of Electrical Engineering, The University of Calgar), Calgary, Alberta, Canada T2N IN4
Received 23 August 1973 The problems of applying the FFT algorithm to the method of Gardner et al. for analysing multicomponent exponential decay curves are discussed. Improvements on Schlesinger's approach are indicated. Schlesinger 1) has shown that the Fast Fourier Transform (FFT) algorithm 2) can be used to modernize the numerical calculations involved in the analysis of multicomponent exponential decay curves using the method of Gardner et al.3). This comment presents a more rigorous mathematical treatment of this method, and draws attention to the problems and literature regarding its implementation via the finite Discrete Fourier Transform (DFT)4). Many research fields (nuclear physics, biology, medicine, magnetic resonance and electrical engineering) deal with experimental observations which may be described by
f(t)= ~ Aiexp(-2tt);
t>0,
(1)
i=1
where the values of A i, 2i and n are to be determined. Mathematically, eq. (1) is a particular case of the more general integral equation f(t)=
g(2) e x p ( - 2 t ) d 2 ;
t>0,
(2)
-oc
[s(x)l = s ( ~ ) s(x) e x p ( - i2 rr/~x)dx.
=
The functional relationship between g(e -x) and x is equivalent to that between #(2)/2 and ln(2). Neither Gardner et al. nor Schlesinger recognized that this solution implies a mapping of integral eq. (2) into the convolution integral equation if(x) = gg(x)®kk(x) =
gg(z) k k ( x - z ) d z ,
(6)
which may be deconvolved for gg(x) by applying Fourier transform techniques: gg(x) = f f - '
{~.~[ff(x)]/~[kk(x)]}.
(7)
The mapping of Gardner et al. is expressed by the relationships gg (x) = g (e- x), k k ( x ) = eXexp(-e~),
provided
(5)
-oc
(8)
f f ( x ) = eXf(e~).
g(2) = ~ A15(2-2~),
(3)
i=1
where 6(2) is the Dirac delta function. Gardner et al. have shown that a solution of eq. (2) is g(e-~) = ~--1 {o~[e~f(e~)]/~-[e ~ exp(_e~)]},
(4)
where f f is the Fourier transform operator defined by 171
The numerical solution of convolution integral equations has been extensively analysed in the literatureS-8), particularly since the advent of the F F T algorithm. The success of the numerical solution is limited by 1) the fact that, in general, 99(x) is not a wellbehaved function, 2) the necessity of approximating the infinite contin-
172
M. R. S M I T H A N D S. C O H N - S F E T C U
uous Fourier transform by a finite discrete Fourier transform, and 3) the use of the F F T algorithm implying the approximation of the aperiodical convolution [eq. (6)] by a periodic convolution defined through the DFT4).
Fig. 1 illustrates the advantage of using a Gaussian rather than a square cutoff filter. The dispersion parameter /~D has been chosen to give the same resolution as Schlesinger's filter. The success of approximating the infinite deconvolution by a finite deconvolution between the limits x = +__X0 is tied to the condition thatff(x), kk(x) and their derivatives tend to zero as x tends to X0 4). This requires X o to be large2). The application of the FFT algorithm to solve eq. (6) means that the resolution of ##(x) is determined by sampling the interval Ax for if(x). A small Ax and large X 0 requires that if(x) be known at a large number of points N, where N = 2Xo/ Ax. Fig. 2 shows the result of decomposing a function of six exponentials. Comparison with fig. 3 of Schlesinger, who was unable to fully resolve four exponentials illustrates the improvement obtained by using a small Ax, large Xo and Gaussian filtering. It is the intention of the authors to pursue the application of this technique to real data, analysing the effects of 1) amplitude noise, 2) time jitter in the sampling device, and 3) interpolation procedures for obtaining if(x), w h e n f ( t ) is not sampled at times t = exp (kAx), where k = 0 , +1, + 2 . . . . , + ( ½ N - l ) , +-iN.
The deconvolution process favours the high frequency components of the data. This enhances the experimental and computational noiseg), obscuring the information carried by gg(x). This phenomenon can be compensated by decreasing the high frequency components of GG(lO=,~[gg(x)]. Whilst such low-pass filtering techniques improve the signal-to-noise ratio of the output function gg(x), they also distort it, leading to a loss of resolution. Schlesinger suggested equating to zero the terms of GG(It) greater than a certain/t 0. This is equivalent to using a square cutoff low-pass filter. The expected series of delta functions for gg(x) becomes convolved with sin(2rrlZoX)/Orx) producing a series of sin (x)/x functions, whose ripples may be misinterpreted as additional peaks. The optimum processing is obtained with a Gaussian low-pass filter, i.e., multiplying GG(lt) by exp (_/~2/p~), leading to a series of Gaussian functions. The Gaussian function, which exhibits the minimum duration-bandwidth product'°), permits the best compromise to be made between loss of resolution and increase in signal-toThe calculations for this letter were performed on a noise ratio of the output signal. This compromise is controlled by the dispersion parameter /~D" A small DEC PDP-11/20 computer kindly loaned by Prof. /~D produces a good signal-to-noise ratio but poor H. A. Buckmaster. The authors wish to acknowledge resolution while a large /~o has the opposite effect. the financial support of the Medical and National I
I
12
I
g( ,l
I
1
I
,o 8 6 4
/
2 C -2 10 °
-2 I0 °
I
I
I
I 0 -r
~0 -~
I0 ~
I0"
I
1
10 -2
iO-S
i 0 -4
), rO -4
Fig. 2. The spectrum g(2)/2 resulting from the decomposition of
f(t) = 10 exp ( - 0 . 0 0 2 t ) + 4 5 exp ( - - 0 . 0 0 9 t ) + 100 exp ( - 0 . 0 2 t ) + Fig. 1. The spectrum g(~.)/2, w h e n f ( t ) = e x p ( - 0 . 0 2 t ) is analysed using square cutoff (@) and Gaussian ( 0 ) filtering. The processing parameters are N = 128, Ax = 0.2 and pD = p o = 24. The ordinate units are arbitrary.
+ 225 exp ( - 0.045 t) + 1000 exp ( - 0.2 t) + 2250 exp ( - 0.45 t). The theoretical spectrum of delta functions is included for comparison. The processing parameters are N = 256, d x = 0.1 and/~r) = 27. The ordinate units are arbitrary.
COMMENT ON FIT TO E X P E R I M E N T A L DATA Research Councils of Canada and The University of Calgary.
References 1) j. Schlesinger, Nucl. Instr. and Meth. 106 (1973) 503. 2) j. W. Cooley and J. W. Tukey, Math. Comput. 19 (1965) 297. 3) D. G. Gardner, J. C. Gardner, G. Laush and W. W. Meinke, J. Chem. Phys. 31 (1959) 978. 4) j. W. Cooley, P. A. W. Lewis and P. D. Welsh, IEEE Trans. Audio. Electroacoust. AU-15 (1967) 79.
173
5) H. D. Helms, IEEE Trans. Audio. Electroacoust. AU-15 (1967) 85. 6) S. Cohn-Sfetcu, M. Caprini and A. Manof, Nucl. Instr. and Meth. 85 (1970) 329. 7) H. F. Silvermann and A. E. Pearson, IEEE Trans. Audio Electroacoust. AU-21 (1973) 112. 8) M. P. Ekstrom, IEEE Trans. Audio Electroacoust. AU-21 (1973) 344. 9) j. E. Welchel and D. E. Guinn, tEEE Trans. Audio. Electroacoust. AU-18 (1970) 159. lo) L. E. Franks, Signal theoo, (Prentice Hall; Englewood Cliffs, N. J., 1969) ch. 6.
INSTRUMENTS
AND
METHODS
II 4 (I974) 171-173; ©
NORTH-HOLLAND
PUBLISHING
CO.
C O M M E N T O N " F I T T O E X P E R I M E N T A L DATA W I T H E X P O N E N T I A L F U N C T I O N S
USING THE FAST FOURIER TRANSFORM" l) MICHAEL R. SMITH Department of Physics, The University of Calgary and the Department of Pathology, Foothills Hospital, Calgary, Alberta, Canada T2N 1N4
and SORIN COHN-SFETCU Department of Electrical Engineering, The University of Calgar), Calgary, Alberta, Canada T2N IN4
Received 23 August 1973 The problems of applying the FFT algorithm to the method of Gardner et al. for analysing multicomponent exponential decay curves are discussed. Improvements on Schlesinger's approach are indicated. Schlesinger 1) has shown that the Fast Fourier Transform (FFT) algorithm 2) can be used to modernize the numerical calculations involved in the analysis of multicomponent exponential decay curves using the method of Gardner et al.3). This comment presents a more rigorous mathematical treatment of this method, and draws attention to the problems and literature regarding its implementation via the finite Discrete Fourier Transform (DFT)4). Many research fields (nuclear physics, biology, medicine, magnetic resonance and electrical engineering) deal with experimental observations which may be described by
f(t)= ~ Aiexp(-2tt);
t>0,
(1)
i=1
where the values of A i, 2i and n are to be determined. Mathematically, eq. (1) is a particular case of the more general integral equation f(t)=
g(2) e x p ( - 2 t ) d 2 ;
t>0,
(2)
-oc
[s(x)l = s ( ~ ) s(x) e x p ( - i2 rr/~x)dx.
=
The functional relationship between g(e -x) and x is equivalent to that between #(2)/2 and ln(2). Neither Gardner et al. nor Schlesinger recognized that this solution implies a mapping of integral eq. (2) into the convolution integral equation if(x) = gg(x)®kk(x) =
gg(z) k k ( x - z ) d z ,
(6)
which may be deconvolved for gg(x) by applying Fourier transform techniques: gg(x) = f f - '
{~.~[ff(x)]/~[kk(x)]}.
(7)
The mapping of Gardner et al. is expressed by the relationships gg (x) = g (e- x), k k ( x ) = eXexp(-e~),
provided
(5)
-oc
(8)
f f ( x ) = eXf(e~).
g(2) = ~ A15(2-2~),
(3)
i=1
where 6(2) is the Dirac delta function. Gardner et al. have shown that a solution of eq. (2) is g(e-~) = ~--1 {o~[e~f(e~)]/~-[e ~ exp(_e~)]},
(4)
where f f is the Fourier transform operator defined by 171
The numerical solution of convolution integral equations has been extensively analysed in the literatureS-8), particularly since the advent of the F F T algorithm. The success of the numerical solution is limited by 1) the fact that, in general, 99(x) is not a wellbehaved function, 2) the necessity of approximating the infinite contin-
172
M. R. S M I T H A N D S. C O H N - S F E T C U
uous Fourier transform by a finite discrete Fourier transform, and 3) the use of the F F T algorithm implying the approximation of the aperiodical convolution [eq. (6)] by a periodic convolution defined through the DFT4).
Fig. 1 illustrates the advantage of using a Gaussian rather than a square cutoff filter. The dispersion parameter /~D has been chosen to give the same resolution as Schlesinger's filter. The success of approximating the infinite deconvolution by a finite deconvolution between the limits x = +__X0 is tied to the condition thatff(x), kk(x) and their derivatives tend to zero as x tends to X0 4). This requires X o to be large2). The application of the FFT algorithm to solve eq. (6) means that the resolution of ##(x) is determined by sampling the interval Ax for if(x). A small Ax and large X 0 requires that if(x) be known at a large number of points N, where N = 2Xo/ Ax. Fig. 2 shows the result of decomposing a function of six exponentials. Comparison with fig. 3 of Schlesinger, who was unable to fully resolve four exponentials illustrates the improvement obtained by using a small Ax, large Xo and Gaussian filtering. It is the intention of the authors to pursue the application of this technique to real data, analysing the effects of 1) amplitude noise, 2) time jitter in the sampling device, and 3) interpolation procedures for obtaining if(x), w h e n f ( t ) is not sampled at times t = exp (kAx), where k = 0 , +1, + 2 . . . . , + ( ½ N - l ) , +-iN.
The deconvolution process favours the high frequency components of the data. This enhances the experimental and computational noiseg), obscuring the information carried by gg(x). This phenomenon can be compensated by decreasing the high frequency components of GG(lO=,~[gg(x)]. Whilst such low-pass filtering techniques improve the signal-to-noise ratio of the output function gg(x), they also distort it, leading to a loss of resolution. Schlesinger suggested equating to zero the terms of GG(It) greater than a certain/t 0. This is equivalent to using a square cutoff low-pass filter. The expected series of delta functions for gg(x) becomes convolved with sin(2rrlZoX)/Orx) producing a series of sin (x)/x functions, whose ripples may be misinterpreted as additional peaks. The optimum processing is obtained with a Gaussian low-pass filter, i.e., multiplying GG(lt) by exp (_/~2/p~), leading to a series of Gaussian functions. The Gaussian function, which exhibits the minimum duration-bandwidth product'°), permits the best compromise to be made between loss of resolution and increase in signal-toThe calculations for this letter were performed on a noise ratio of the output signal. This compromise is controlled by the dispersion parameter /~D" A small DEC PDP-11/20 computer kindly loaned by Prof. /~D produces a good signal-to-noise ratio but poor H. A. Buckmaster. The authors wish to acknowledge resolution while a large /~o has the opposite effect. the financial support of the Medical and National I
I
12
I
g( ,l
I
1
I
,o 8 6 4
/
2 C -2 10 °
-2 I0 °
I
I
I
I 0 -r
~0 -~
I0 ~
I0"
I
1
10 -2
iO-S
i 0 -4
), rO -4
Fig. 2. The spectrum g(2)/2 resulting from the decomposition of
f(t) = 10 exp ( - 0 . 0 0 2 t ) + 4 5 exp ( - - 0 . 0 0 9 t ) + 100 exp ( - 0 . 0 2 t ) + Fig. 1. The spectrum g(~.)/2, w h e n f ( t ) = e x p ( - 0 . 0 2 t ) is analysed using square cutoff (@) and Gaussian ( 0 ) filtering. The processing parameters are N = 128, Ax = 0.2 and pD = p o = 24. The ordinate units are arbitrary.
+ 225 exp ( - 0.045 t) + 1000 exp ( - 0.2 t) + 2250 exp ( - 0.45 t). The theoretical spectrum of delta functions is included for comparison. The processing parameters are N = 256, d x = 0.1 and/~r) = 27. The ordinate units are arbitrary.
COMMENT ON FIT TO E X P E R I M E N T A L DATA Research Councils of Canada and The University of Calgary.
References 1) j. Schlesinger, Nucl. Instr. and Meth. 106 (1973) 503. 2) j. W. Cooley and J. W. Tukey, Math. Comput. 19 (1965) 297. 3) D. G. Gardner, J. C. Gardner, G. Laush and W. W. Meinke, J. Chem. Phys. 31 (1959) 978. 4) j. W. Cooley, P. A. W. Lewis and P. D. Welsh, IEEE Trans. Audio. Electroacoust. AU-15 (1967) 79.
173
5) H. D. Helms, IEEE Trans. Audio. Electroacoust. AU-15 (1967) 85. 6) S. Cohn-Sfetcu, M. Caprini and A. Manof, Nucl. Instr. and Meth. 85 (1970) 329. 7) H. F. Silvermann and A. E. Pearson, IEEE Trans. Audio Electroacoust. AU-21 (1973) 112. 8) M. P. Ekstrom, IEEE Trans. Audio Electroacoust. AU-21 (1973) 344. 9) j. E. Welchel and D. E. Guinn, tEEE Trans. Audio. Electroacoust. AU-18 (1970) 159. lo) L. E. Franks, Signal theoo, (Prentice Hall; Englewood Cliffs, N. J., 1969) ch. 6.