Improved Multiexponential Transient Spectroscopy by Iterative Deconvolution

—The analysis of multiexponential decays is challenging because of their complex nature. When analyzing these signals, not only the parameters, but also the orders of the models, have to be estimated. Here, we present an improved spectroscopic technique specially suited for this purpose. The proposed algorithm combines an iterative linear filter with an iterative deconvolution method. A thorough analysis of the noise effect is presented. The performance is tested with synthetic and experimental data.


Improved Multiexponential Transient Spectroscopy By Iterative Deconvolution Santiago Marco, Jordi Palacín, and Josep Samitier
Abstract-The analysis of multiexponential decays is challenging because of their complex nature.When analyzing these signals, not only the parameters, but also the orders of the models, have to be estimated.Here, we present an improved spectroscopic technique specially suited for this purpose.The proposed algorithm combines an iterative linear filter with an iterative deconvolution method.A thorough analysis of the noise effect is presented.The performance is tested with synthetic and experimental data.

I. INTRODUCTION
M ULTIEXPONENTIAL signals are found in such diverse areas as nuclear science, chemical kinetics, biomedicine, semiconductors, thermal responses, etc.Many different methods have been proposed for the analysis of these signals: nonlinear least-squares fitting [1], the method of moments [2], linear prediction [3], exponential series method [4], Gardner transform [5], the Prony algorithm [6], or the Padé-Laplace transform [7] are some examples.However, we are facing a still unsolved problem [8].The underlying difficulty is that exponential decays, in contrast to pure sinusoids, do not form an orthogonal basis of functions on the real axis.This fact makes the estimation process extremely sensitive to experimental noise and truncation of the measured signal.
The authors have proposed a spectroscopic technique for the analysis of these kinds of transient signals: multiexponential transient spectroscopy (METS) [9].It is important to note that the spectroscopic nature of the algorithm gives more information than other techniques because i) permits dealing with discrete and continuous time distributions; ii) shape of the peaks gives information about the adequacy of the model; iii) number of peaks gives information about the number of exponential terms to include in the model.Moreover, as it will soon follow, the procedure leads to the natural implementation of iterative deconvolution algorithms that can improve the separation capabilities of the technique.A final step based on constrained nonlinear least squares furnishes the final numerical estimates of the unknown parameters.
In the following sections, we will present, for the sake of completeness, the basic theory of the algorithm, an analysis of the noise influence, details of the discrete-time implementation, the introduction of an additional nonlinear deconvolution step, and some application examples with synthetic and experimental signals.The paper will end with the main conclusions.

II. THEORY
We consider the set of signals generated from a certain timeconstant distribution (1) The objective is to recover the underlying .The problem of exponential analysis is solved by taking the inverse Laplace transform of the transient.This is only feasible if the analytical expression of is known.However, this is not the case with experimental signals.The authors have proposed the following approach for the solution of (1).The first-order METS signal is defined

METS
(2) It can be deduced that this signal can be written as a convolution

METS
(3) with , , and , , is a peak-shaped kernel, and contains the information to be recovered.On most occasions this technique is applied to the analysis of decays arising from a discrete time constant distribution.In this case, the METS signal will consist of a finite set of peaks whose shape will be given by the kernel function, because the METS signal is the convolution of a set of -Dirac with the kernel function.However, due to the finite width of the peak, blurring will fuse near time constants in a single broader peak (Fig. 1).This signal analysis brings up immediate analogies with spectroscopic techniques.In this case the "instrument impulse response" can be identified with the kernel function.
The recovery of from the measured METS signal is a complex deconvolution problem.In the frequency domain, the kernel function acts as a lowpass filter, removing high frequencies present in .The Fourier transform of can be cal- At high frequencies, the noise level might be higher than the signal, and subsequent attempts to recover have to face noise amplification that can easily hide the information in the deconvolved results.Up to this point, this theory is similar to the Gardner transform method.However, in the latter, deconvolution is accomplished using the Fourier transform

METS
(5) Since is a lowpass signal, this division enhances the higher frequency components of METS extremely.This is essential for resolution increase, but it produces an unwanted increase of the noise level, as well.Usually, to regularize this deconvolution, a Gaussian lowpass filter is applied before the inverse Fourier transform is applied (6) where the parameter is used as a regularizing parameter.An example of this procedure can be found in [10].
Here, we will present an iterative linear deconvolution technique specially suited for this problem.Additionally, we propose to apply nonlinear iterative deconvolution techniques that usually improve the results achieved with linear deconvolution methods.
The iterative linear filter is given by METS METS METS (7) The cumulative application of this filter acts as a high-pass filter.Fig. 2 shows the evolution of the filter frequency response.As increases, more weight is given to the high-frequency content of the signal.
One of the most interesting features of this iterative filter is that it preserves the METS signal as a new convolution of with a sharper kernel signal.

METS (8)
where (9) is the Gamma function.All the kernels have unity area.The shape of these kernels can be observed in Fig. 3. Time resolution increases as the user computes higher order METS signals.This can be observed in the simulation example (Fig. 4), where the signals appear after normalization and time shift.The analyzed signals consist of four time constants whose amplitudes can be observed in the ninth-order METS signal.An additional advantage of the convolution preservation feature is the possibility to discard spurious peaks with reduced width that can not be due to new time constants.
As it can be observed in Fig. 3, the peak becomes higher and narrower.The maximum value is (10) Note that, although not apparent in the figure, the peak suffers a controlled shift toward higher times.This is not a problem, since the METS signal can be shifted ' ' to the left for interpretation purposes or visual representation.However, a side consequence of such a behavior is that the maximum METS order, which can be calculated, is given by the length of the acquired transient divided by the longer time constant present in the signal.Using the Stirling formula [11], it can be deduced that the integral breadth, for high-order METS signals, is .This is an approximation of the time resolution attainable at a certain METS order.This slow decrease in the kernel width indicates that it would not be worthwhile to compute high orders because the resolution gain would not compensate the noise increase.
Because of the convolution preservation property of the present approach, at any time it is possible to use other deconvolution algorithms such as, for instance, Jansson's deconvolution [12].

III. THEORETICAL NOISE ANALYSIS
The ultimate time resolution in exponential analysis is noise limited.Within the presented framework, deconvolution increases the noise content and hinders the separation of close time constants.A detailed study by Bertero [13] on the eigenfunctions and eigenvalues of the Laplace transform (which were previously obtained by Whirter and Pike [14]) has shown the following: The closest distance between exponential decay rates that can be resolved in exponential analysis when the signal is observed in an infinite domain [0, + ] is given by SNR (11) where the SNR is defined as the transient amplitude divided by the noise standard deviation.The authors have also studied the resolution limits, in particular cases, based on the Cramér-Rao theorem [15].However, a thorough comparison of both approaches has not yet been performed.According to (11) the evolution of the theoretical resolution limit against the SNR may be observed in Fig. 7.As expected, the higher the SNR the better the resolution.For very low noise content (SNR 10 ), this ratio saturates to about 1.22, while for moderate noise content (SNR 10 ), it still gives a 2.45 limit.
A simple approach to the study of the effects of noise in our method is to suppose that the signal contains additive white noise in the transformed domain .Although the real situation is more complex, this simple analysis enlightens the main effects.Our linear deconvolution approach amplifies highfrequency noise, imposing a practical limit to the maximum order that can be computed before the degradation of the SNR prevents further interpretation of the results.Fig. 5 shows the noise gain depending on the METS order and the number of points per decade.In light of the results some conclusions can be drawn.
ii) The noise gain is quite high.
iii) It increases with the METS order, due to the iterative application of a differentiation step.iv) The noise increase is more important as the number of points per decade increases.It can also be observed that the noise gain tends to saturate as the order increases, probably due to the 1/n factor in the differentiation term.This saturation appears at higher orders for higher sampling rates in the transformed domain.
It may appear that given a certain noise floor, Fig. 5 permits estimating which is the highest order attainable and in consequence the time resolution, the only condition being that the SNR remains high enough.However, these arguments disregard the presence of distortion.The METS signals, although of a lowpass nature, are not strictly bandlimited, and, in fact, the high-frequency contents increase with the order.The analytical expression of the spectrum (see Fig. 6) is given by (12) At this point, it is perhaps curious to realize that the application of the METS algorithm, when observed in the frequency domain, corresponds to the application of the Pochhammer symbol property of the Gamma function [11].
To avoid distortion, the sampling rate in the transformed domain has to be high enough so as to permit a correct differentiation of the signal.For instance, for SNR of 10 and ten points per decade, the METS order cannot increase indefinitely, because, at a certain point, distortion impedes further resolution enhancement.

IV. DISCRETE TIME IMPLEMENTATION AND NOISE FILTERING
For the digital implementation of the above algorithm, a very critical step is the transformation from a uniformly sampled transient in time to a uniform sampling in the logarithm of time.Such transformation is carried out together with low-pass filtering.We have to consider that, usually, the experiment needs several time decades to be explored (five to eight decades are common numbers depending on the application).Because the usual hardware for digital signal acquisition only permits uniform time sampling, and to prevent a large amount of data from collapsing our signal processing system, usually signal acquisition is done in several sections with a progressive reduction of the sampling rate.From this nonuniform signal sampling, we interpolate the signal values at times (13) where desired number of points per decade; minimum sampling period in the time domain; sampling frequency in the transformed time domain.For minimum noise plus distortion, the authors have proposed an optimum filter of varying length.This filter is described in detail in [16].A consequence of this procedure is that the noise in the filtered METS signal becomes nonstationary with a variance that decreases with the inverse of time.This is why the noise analysis approach of the previous section can be regarded only as an approximation.
For the calculation of the derivative we use a simple centered two-point difference [9] (14) Of course, more elaborate differentiators can be used [17], but this is well suited to the application because of its minimum length which minimizes border effects at the beginning and final part of the transient.In addition, it also attenuates the high-frequency gain of a perfect differentiator.However, it will introduce distortion earlier because it does not differentiate the whole frequency range.
The evolution of the time resolution of the METS analysis, depending on the signal-to-noise ratio and the points per decade of the signal, may be observed in Fig. 7.This numerical experiment has been carried out with a biexponential signal where both components have the same amplitude.As we have already mentioned, the optimum filter described in [16], induces a nonstationary noise in the METS.To avoid the influence of this effect, synthetic Gaussian white noise of different power has been added before the METS analysis.No filtering, only linear interpolation, has been carried out in the nonlinear time scale transformation.In this case we obtain an upper limit of the time resolution that can be improved by filtering.As expected the ultimate resolution limit, achieved with the METS, follows the same trend given by the theoretical limit but at higher decay rate ratios.Moreover, this degradation is more evident at the lower SNR due to the noise amplification produced by the linear filter.It may be observed that the resolution limit depends on the number of points per decade selected from the initial transient.For the higher noise content, an increase in the number of points per decade produces worse results in agreement with = p =t .Fig. 6.The signal energy becomes concentrated in the lower frequency range where the METS filter has very small gain, while the high-frequency noise is greatly amplified.For the low noise content, the opposite behavior is observed: the higher the sampling frequency, the better the resolution.For very low noise, the maximum resolution achieved for ten points per decade is distortion limited.A higher sampling frequency (for instance, 40 points per decade) allows further improvement in this region.
We have already mentioned that the noise content becomes nonstationary and decreases with time.As a consequence, the time resolution depends on the position of the time constant; more specifically, on the ratio between the time constant and the initial sampling period t0.Results may be observed in Fig. 8.The higher the ratio of the time constant to sampling period, the lower the noise content at the time constant position and the better the resolution.This effect is more evident at lower SNR.

V. ADDITIONAL ITERATIVE DECONVOLUTION
As we have already mentioned, the proposed procedure preserves the signal as a convolution, so we can use alternative deconvolution schemes intending to improve the resolution.Although a great variety of deconvolution methods have been reported in the literature, we will restrict ourselves in this work to the application of the Jansson's method.The Jansson's method is a constrained, nonlinear deconvolution method growing from an evolution of previous linear deconvolution techniques coming from the work of Van Cittert.Its rationale comes from the fact that linear methods sometimes end up in nonphysical results.An example of such a case is the analysis of a thermal one-port system.In this case, only positive amplitudes have physical meaning [18].By forcing the solution to be physically feasible, the overall performance improves.Not only negative peaks can be prevented, but by proper choice of the relaxation factor, the maximum peak amplitude is also limited.In such cases the Jansson's method can be written as METS (15) ( stands for convolution) where the relaxation factor including the above constraints is given by ( 16) " " is a parameter which can be adjusted typically by trial and error, although it is related to the eigenvalues of the impulse response on the Toeplitz matrix form.We have used values between 0.02 and 0.002.First, we present two synthetic examples to illustrate the performances of this technique.The analyzed transients consist of several exponential terms with added Gaussian noise with variance 10 (SNR 10 ).The sampling frequency is 10 kHz, and the total length of the transient is 10 s.The Jansson's method has been applied with = 0.02, and the total number of iterations is 15 000.The obtained results can be seen in Table I.In both examples, peaks separated up to 1.6 octaves (ratio 3) are identified as distinct time constants.However, the obtained peaks are not perfect unity impulses; they present a finite width.To measure the amplitude, we have used the integral amplitude given by the sum of the samples within every peak.Time constants have been estimated from the peak maximum.However, at this stage, neither the positions nor the amplitudes are accurate.This fact has induced in our analysis procedure a final constrained nonlinear least squares fit of the first-order METS signal to estimate the unknown parameters, relying on the deconvolution only for model order determination and for initial values of the parameters.In Table I, see the important increases in the accuracy obtained in this final step.Concerning the evolution of the time resolution against the SNR, we have observed that the Jansson method provides only a small improvement compared to the METS method (see Fig. 7).However, visual inspection of the deconvolution results is easier to interpret.When time resolution is feasible, the Jansson deconvolution provides clearly distinct peaks, while METS provides only a shallow minimum in between.
Here we present the analysis of some experimental results obtained when studying the behavior of thermal microactuators.In thermal systems, the usual way to characterize the system dynamics is the analysis of the thermal response to a step in the applied power.The so-called transient thermal impedance is given by (17) where recorded temperature; reference temperature (usually ambient); amplitude of the power step.This signal consists of a superposition of exponential decays, so that the recovery of the amplitudes and time constants provides the complete information on the thermal transfer function (18) In this particular case, the device under study was a thermopneumatic micropump [15].To study the thermal behavior, a  negative power step (1 W) was applied to the heating resistor, and the temperature of the device was monitored from the same resistance that was previously calibrated as a temperature detector.The transient had total amplitude of 107 C, and the noise floor was about 0.1 C. Fig. 9 shows the results obtained with the METS and the Jansson method.With the METS order eight of this transient signal, one time constant has been identified, but there is a wide peak clearly formed by various time constants that are beyond the practical limits of the METS.Applying the Jansson method, two hidden time constants are now identified.So, as indicated in Fig. 7, the Jansson method (applied to the METS order one) is an improvement of the METS algorithm that avoids visual interpretation of the results because the time constants are clearly identified by narrower peaks.The result of the nonlinear fit is observed in Fig. 10, while the estimated amplitudes and time constants appear at the bottom of Table I.
Although our method shows time-resolutions similar to the Gardner transform [10], [19], [20], it is very difficult to present a consistent comparison.It should be taken into account that there are many factors which affect, in a diverse manner, the performances of the methods.Among them, we can find the noise level, the presence of an offset, the number of exponential components, the ratio of their amplitudes and time constants, the ratio between the sampling time and the time constants, and, finally, the ratio between the total measurement time and the slowest component.A rigorous comparison of the methods requires that all the methods are tested in the same conditions and, moreover, the analyzed transients should cover a wide range of the above parameters.Although some comparisons are present in the literature, in agreement with [8], no one fulfills the above conditions, and a more exhaustive comparison with other methods is beyond the scope of this work.

VI. CONCLUSION
In this paper, we have presented a linear method for the analysis of multicomponent exponential transients.This method, combined with iterative nonlinear deconvolution methods, permits the identification of the model order and provides useful initial values of the signal parameters.We have studied the time resolution of our proposal for different SNR, and it has been compared with the theoretical limit.Obviously, the obtained time resolution is worse than the theoretical limit.However, it is similar or even slightly better than results obtained by the Gardner transform method.The shape preservation characteristic of the proposed algorithm permits the final estimation of the parameters that is carried out by nonlinear fitting of the METS signal.A more complete analysis of the metrological performances of this technique needs further work.

Fig. 7 .
Fig. 7. Evolution of the resolution limit against the SNR: theoretical limit (continuous line), METS with 40 points per decade (circle), METS with 10 points per decade (triangle), and Jansson with 10 points per decade (square).

Fig. 8 .
Fig. 8. Evolution of the resolution against the position of the time constant.

TABLE I ESTIMATION
RESULTS OF TWO SYNTHETIC DATA EXAMPLES AND ONE THERMAL TRANSIENT FROM EXPERIMENTAL DATA