Suboptimal Filtering and Nonlinear Time Scale Transformation for the Analysis of Multiexponential Decays

Multiexponential decays may contain time-constants differing in several orders of magnitudes. In such cases, uniform sampling results in very long records featuring a high degree of oversampling at the final part of the transient. Here, we analyze a nonlinear time scale transformation to reduce the total number of samples with minimum signal distortion, achieving an important reduction of the computational cost of subsequent analyses. We propose a time-varying filter whose length is optimized for minimum mean square error.


I. INTRODUCTION
T HE ANALYSIS of multiexponential decays appears in many areas of scientific or technical interest: analysis of deep levels in semiconductors [1], thermal responses [2], fluorescence spectroscopy [3], etc.A recent review confirms the impact of this family of signals and reviews the variety of analysis methods [4].In many cases, the transient signal under study may contain time-constants differing in several orders of magnitudes.In some applications such as deep level spectroscopy or thermal analysis, as many as eight time decades are recorded.With linear sampling, the transient would contain 10 points.Such a huge number of points represent an incredible computational burden, and they are the consequence of a high degree of oversampling in the final time decades.Decaying signals are in fact nonstationary, and their bandwidth decreases with time.In qualitative terms, the "instantaneous bandwidth" of a deterministic signal is related to its derivatives.It is commonly accepted that in the regions of maximum curvature the required bandwidth is maximum.However, for exponential decays all derivatives decay with time, and so also the required bandwidth.In consequence, it is not necessary to keep a very high sampling rate to capture the fastest time-constants occurring at the initial part of the transient to analyze the slowest time-constant that appears at the final part.In different application areas, researchers usually prefer a logarithmic representation (Fig. 1) of the time scale which shows up more clearly the presence of fast time constants which appear hidden in a linear time scale plot.In fact, several algorithms, based on the work of Gardner et al. [5], have been proposed which directly operate in this transformed time scale [4], [6].This is because after the time transformation every exponential term presents the same signal shape whatever the time constant is.In other words, an increase (decrease) in the time-constant only represents a right (left) shift.In other words, the transformed signal can be expressed as a convolution [5].The authors have already analyzed in previous works how to take advantage of this invariance [7].
Our purpose in this work is to analyze how to perform this time scale transformation from linear sampling in the time scale to equidistant samples in the logarithm of time for minimum signal distortion in the presence of additive white noise.

II. PROCEDURE
As already mentioned in the previous section, our aim is to carry out a data reduction algorithm that maps the linear time scale to a logarithmic time scale.We will consider first the effect of such transformation in a single exponential (1) A point to discuss concerns the adequate sampling rate in the transformed time domain.In the following we will talk about pseudo-frequency and pseudo-spectrum because the independent variable is no longer time, but the logarithm of time.It can be seen that the pseudo-frequency spectrum of an exponential in the transformed domain is clearly low pass.It can be calculated analytically, and it is related to the Gamma function evaluated along the imaginary axis (2) It may be observed that the pseudo-spectrum amplitude does not depend on the time constant value, because a change in the time constant only represents a shift in position in the transformed domain.A few points per decade are necessary to adequately represent the signal, although the algorithm presented later on permits the calculation of the signal in the transformed domain for a large range of sampling rates (in the transformed domain).Our experience indicates that sampling rates exceeding ten to fifteen points per decade produce high-frequency pseudo-spectral bands where the noise content is more intense that the signal power.However, a rigorous analysis of this problem will be published elsewhere [8].
New samples are equidistant in the logarithm of time (3) where is the initial time, and is the number of points per decade.Note that to represent eight time decades we will need only 80-120 points as compared with 10 .
This "decimation" in time must be accompanied by low-pass filtering to avoid aliasing of the superimposed noise that has, in general, a much wider spectrum.Note that this process resembles, but it is not the same as, the usual sampling rate conversion by a rational factor for which a complete theory is available [9].In addition, the actual time points do not coincide with the original sampling times, so we must necessarily interpolate.Thus, the procedure will consist of: first, filtering for the minimum square error followed by linear interpolation to calculate the signal value at the points .Note that the new time instants do not depend on the filtering method chosen.
Concerning the filtering step, we will propose three different approaches.The first one is an heuristic method, and it is considered as the baseline algorithm.The second is based in an optimum filter introduced by Papoulis [10].The third procedure is a time-varying filter, the length of which is optimized for minimum mean square error.

A. Heuristic Method
The heuristic method is as follows: new time points are generated in the logarithmic time scale.Typically the sampling rate in this transformed domain is chosen to be between 10 and 20 points per decade.In the computation of each new sample we use all original samples ( samples) within the interval ( ) where .Moreover, the filter length increases linearly with time, and consequently the cutoff frequency decreases linearly (4) where " " is the sample index in the linear domain.To calculate the new sample we use polynomial linear or quadratic interpolators, obtained after least squares fitting of the points within the interval to a straight line or parabola, respectively.This polynomial filtering can be considered as a generalization of the Savitzky-Golay SG-FIR filters [11], which provide linear phase, and maximum noise rejection ratio subjected to the condition of frequency flatness at DC.The main difference is: since the new sample does not coincide with the original, the filter coefficients may vary slightly depending on the parameter which represents the fractional position of the new sample within the original time sampling period .In this way, we can write the filter operation as (5) with and, , and the filter coefficients that are linear or parabolic functions of if we are considering first-or second-order Savitzky-Golay filters, respectively.Their closed forms are not written here due to space constraints.
The first consequence of the above procedure is that the noise signal becomes nonstationary and its variance decreases with time.Very important noise reduction is achieved in the final time-decades because of the very low cutoff frequency of the interpolating filter, which is a consequence of the huge number of points within the interval.From this point of view the process is very advantageous because it provides simultaneously significant data reduction and reduces the noise content.
However, at this point we must consider the signal distortion that can arise after this time-varying filtering and "decimation."As real exponentials are not bandlimited signals, any filtering may induce some degree of signal distortion and, in fact, this is the case.Figs. 2 and 3 show the mean square error which appears when filtering a single exponential decay with the first-order modified S-G filter according to the above procedure.In Fig. 2, we observe the distortion depending on the time constant value: maximum distortion appears in the vicinity of the time constant value.In Fig. 3, we may observe that this distortion is more apparent for low noise levels.

B. Papoulis Filter
Up to now, we have seen the limitations imposed by the first procedure.To introduce our second proposal, we wonder if the filter length evolution along the transient is optimal.Papoulis [10] has obtained the optimal filter length when filtering a deterministic signal embedded in white noise, and he has found that it depends on the noise content and on the second derivative of the noise signal (6) where original signal; filtered signal; optimal filter length; linear sampling period; noise variance; second derivative of the underlying noise-free transient.The practical application of this method in our case has a problem: the second derivative has to be estimated.It is well known that differentiation of noisy signals is a risky signal processing operation [12]- [14].A variety of methods have been proposed in the literature.Among them we have obtained the best results by double differentiation of a signal model obtained by fitting an exponential series [3] to our original data (ESM).Such series contain three exponentials per decade at fixed logarithmically distributed times.Only the amplitudes are estimated by linear least squares.The obtained fitting is very good, and the residuals are white.Fig. 4 shows the goodness of the fitting, and Fig. 5 compares the real and estimated amplitudes for a biexponential decay.Despite the very good time response, we still suffer from unacceptable errors in the estimation of the second derivative at the initial part of the transient (Fig. 6).These errors are undoubtedly related to the inadequacy of the model used.Poor estimation of the second derivative will lead to erroneous filter lengths that will provide either distortion or insufficient noise rejection.Because of that  we have devised a method which avoids the calculation of the instantaneous second derivative of the transient.
Even when the second derivative is known, ( 6) is only valid when (see the original demonstration for further details [10]), where is the sampling period and is one of the time constants of the analyzed transient.

C. Suboptimal Method
In this method, we carry out first a varying time filtering in the linear domain.Transformation to the logarithmic scale is performed after filtering by a simple linear interpolator.
The error can be theoretically evaluated with (7) where is the noise content, and is the noise-free multiexponential decay.We suppose now that a multiexponential model of the signal is available, which can be obtained as it has been previously outlined (8) where filter impulse response; amplitudes; .Then, the error power can be estimated for a general case.Note that we consider a centered filter for zero group delay (9) where (10) Expression (9) shows that the error has two contributions: distortion and filtered noise.
Also, when an estimation of the noise and the exponentials contained in the signal is available, (9) can be used to find the optimum length of the filter in each point to achieve the minimum error power in the filtered signal.We have minimized the above expression in the case of a moving average filter (length , filter coefficients ) (11) An algorithm based on golden section search and parabolic interpolation [15] is used to find the local minima corresponding to the optimum filter length .The computational burden of this iterative minimization is highly reduced if the last value found is provided as the initial value in the next time step.
This kind of filter provides the maximum noise rejection ratio.Moreover for the large filter lengths obtained at the end of the transient, it is not practical to try to optimize the filter coefficients, because huge matrices should be inverted.Note that the optimum length depends on the sample index " ," so we obtain again a time-varying filter as expected due to the signal characteristics.The proposed filtering method is optimal,   Finally, we would like to remark that due to the final linear interpolation step, it is not necessary to filter the complete transient but only those samples needed for the final interpolation.This reduces drastically the number of floating point operations needed to complete the algorithm.

D. Discussion
Results in Figs.7 and 8 show the evolution of the filter length and the mean square error with time.The heuristic, the optimal and the suboptimal methods are compared.The estimated signal model is obtained as it has been previously outlined (Exponential Series Method) for the second derivative estimation.It can be seen that, despite the signal model parameters being clearly far from reality, the predicted filter length with the estimated model (suboptimal) is almost identical to the exact result (optimal).This very good result is obtained with the same signal model that produced problems in the second derivative estimation.Moreover, the results clearly outperform the mean square errors obtained with the heuristic method: the noise rejection is larger and distortion peaks are avoided.The large increase in the error at the final part of the transient is a border effect.

III. CONCLUSIONS
It is not practical to work with linear sampled multiexponential signals containing time-constants differing in several orders of magnitudes.A procedure has been proposed to achieve a dramatic data reduction and noise reduction based on a moving average filter, the length of which is optimized for the minimum mean square error.However, a previous estimation of a signal model and noise level is needed.This estimation is based on the linear least squares of an exponential series with three terms per decade.Despite the fact that the signal model parameters could be quite far from reality, very good time agreement is obtained and the filter length prediction is almost optimal, providing maximum noise rejection with minimum distortion.

Suboptimal
Filtering and Nonlinear Time Scale Transformation for the Analysis of Multiexponential Decays Jordi Palacín, Santiago Marco, and Josep Samitier Abstract-Multiexponential decays may contain time-constants differing in several orders of magnitudes.In such cases, uniform sampling results in very long records featuring a high degree of oversampling at the final part of the transient.Here, we analyze a nonlinear time scale transformation to reduce the total number of samples with minimum signal distortion, achieving an important reduction of the computational cost of subsequent analyses.We propose a time-varying filter whose length is optimized for minimum mean square error.Index Terms-Data reduction, Gardner transform, multiexponential decay, noise, nonuniform sampling, transients.
Authorized licensed use limited to: Universitat de Barcelona.Downloaded on February 16, 2009 at 07:11 from IEEE Xplore.Restrictions apply.
Authorized licensed use limited to: Universitat de Barcelona.Downloaded on February 16, 2009 at 07:11 from IEEE Xplore.Restrictions apply.