Fast and Accurate Fitting and Filtering of Noisy Exponentials in Legendre Space

The parameters of experimentally obtained exponentials are usually found by least-squares fitting methods. Essentially, this is done by minimizing the mean squares sum of the differences between the data, most often a function of time, and a parameter-defined model function. Here we delineate a novel method where the noisy data are represented and analyzed in the space of Legendre polynomials. This is advantageous in several respects. First, parameter retrieval in the Legendre domain is typically two orders of magnitude faster than direct fitting in the time domain. Second, data fitting in a low-dimensional Legendre space yields estimates for amplitudes and time constants which are, on the average, more precise compared to least-squares-fitting with equal weights in the time domain. Third, the Legendre analysis of two exponentials gives satisfactory estimates in parameter ranges where least-squares-fitting in the time domain typically fails. Finally, filtering exponentials in the domain of Legendre polynomials leads to marked noise removal without the phase shift characteristic for conventional lowpass filters.


Introduction
Processes with linear first order kinetics, i.e. dy=dt~{ay, are ubiquitous in all fields of science, above all in the life sciences, chemistry, physics and engineering. Noisy exponentials are therefore among the most common experimental outcomes. Typically, noise removal is done by conventional lowpass filtering, and the function's parameters, ie, amplitudes and time constants, are retrieved by nonlinear least squares fitting (NLLSQ) such as the Levenberg -Marquardt algorithm (LMA) [1][2][3]. However, lowpass filters in the Fourier domain distort the signals by introducing phase shifts, and the LMA is too time-consuming in cases where the parameters need to be obtained rapidly or, equivalently, where exponentials need to be fitted in parallel. We therefore set out to find a faster method of filtering and fitting exponentials with the same or higher precision. The goal was to transform noisy signals such as exponentials into a space where signal and noise are mapped essentially onto different subspaces. It turned out that the method described herein is not only faster but in many practical applications also more precise than the LMA. In addition, it allows effective noise removal.

Legendre filter
Consider the outcome of the stopped-flow experiment shown in Fig. 1a as well as its Legendre spectrum (b), obtained by the finite Legendre transform (fLT, see Methods). As Legendre polynomials are orthogonal on the interval [21, 1] (Fig. S1), the first and the last sample of the time function are assumed to occur at -1 and 1, respectively (for scaling and re-scaling of time, time constants and amplitudes, see Methods). The Legendre spectrum (b) is plotted for the first 17 components, each of which indicates the contribution of a specific Legendre polynomial to the signal (eq. 4). While the inverse finite Legendre transform (ifLT) of the entire spectrum would, of course, reconstruct the original noisy function (a), the ifLT of the first, e.g. eight, components of the spectrum results in considerable noise removal (a, red). In analogy to Fourier theory, and for the sake of brevity, we name this filter Legendre lowpass. Its effectiveness is due to the fact that smooth signals such as exponentials are mainly represented by lower Legendre polynomials, i.e., by lower order components of the Legendre spectrum, while the signals' noise is mapped predominantly to higher order components ( Fig. 1 and Fig. S2).
In practical cases, the output of a linear experimental device is affected by the system's response function, r(t). Often the response function is a lowpass so that r(t) is an exponential itself. This can be neglected if its time constant t r is much smaller than the t under investigation (as in the case of Fig. 1a). Generally, however, an experimental outcome y(t) is the convolution of the signal, x(t) (Fig. 2a), with the response function, i.e., y(t)~x(t) Ã r(t). Fig. 2e and f show that the convoluted exponential y(t) (e, noisy trace) can be approximated by Legendre filtering (f, gray, and e, red, cont.). Moreover, in many cases one is interested in the parameters of the original, non-convoluted function x(t) rather than in y(t).
Obtaining these in the time domain (t-domain) from y(t) would require a deconvolution, which is notoriously inaccurate for noisy functions, and often not feasible. In contrast, an approximation of x(t) can readily be obtained in the Legendre domain (L-domain). To this end we calculate the Legendre spectrum of x(t) from r(t) and the Legendre spectrum of y(t) using eq. 8. The inverse transform of its lowpass-filtered spectrum (f, red bars) gives an approximation of x(t) (e, red, dashed; eq. 5).
The marked reduction of noise observed could also be observed when a sum (Fig. 3) or a product (Fig. 4) of noisy exponentials was analysed. The product of exponentials could serve as a model, e.g., for excitatory postsynaptic potentials (EPSPs) of neurons.
Here it is interesting to compare the Legendre lowpass-filtered EPSP (Fig. 4c, dashed, red, calculated from the lower components of the EPSP's Legendre spectrum, b) with the conventionally lowpass-filtered EPSP, which exhibits the characteristic phase shift of Fourier lowpass-filters (c, gray). For obtaining parameters such as time-to-peak and time constants, the Legendre lowpass clearly appears to be more appropriate than the Fourier lowpass.
The noise removal effect of Legendre filters is by no means limited to the above functions as shown for the autocorrelation function (ACF) of an fluorescence correlation spectroscopy (FCS) experiment (Fig. 1c). Its Legendre spectrum as well as the Legendre lowpass-filtered ACF are shown in d and c (red), respectively. Though this cannot replace the fitting of FCS data, it allows a filtered online representation of the ACF.
Taken together, the finite Legendre transform compresses the data from the number of samples in the t-domain to a small number of Legendre components. Second, considerable noise reduction can be achieved by Legendre lowpass-filtering, and third, deconvolution of noisy exponentials can conveniently be carried out in Legendre space.
At this point the question arises whether the parameters A and t of noisy exponentials can be retrieved in the L-domain, and whether this is faster and as least as precise as fitting in the tdomain. Our strategy to analyze these questions is illustrated in Fig. S2. A noisy exponential, which might look similar to part g of the figure, is first transformed (fLT) to its Legendre spectrum. The lower components of the Legendre spectrum are then fitted to the corresponding, amplitude-and time constant-dependent components of a pure exponential (eq. 4). For the comparison of this method with the LMA in the t-domain, we simulated noisy exponentials (g). First, we defined the amplitude and time constant of an exponential (a) and then we added Poisson noise (c) as well as an offset with gaussian noise (e) to it. The right column of the figure clearly shows that the Legendre spectrum of (a) has only low components, while the transforms of (c) and (e) are characterized predominantly by high components. The LMA fit of the resulting noisy exponential in the t-domain (g) now gives estimates A t and t t for the true amplitude A and true time constant t, while the LMA fit of its lower Legendre components (h) gives estimates A L and t L for A and t in the L-domain. . For the experimental details of the system, see [12]. (b) The first 17 components of the Legendre spectrum of x(t). The inverse fLT of the components L o through L 7 of the spectrum (b) gives the red curve in (a). Note that the sharp peak in the noisy trace is virtually not reflected in the filtered curve. (c) Autocorrelation curve (gray) resulting from an experiment where the diffusion constant of tetramethylrhodamine was measured (own data). In this example, fLT and ifLT are performed for nonequidistant samples, and we re-scaled the x-axis to correlation delays.

L-domain
To assess and compare the precision and accuracy of fitting in the t-domain versus the L-domain, we used simulated data as described and applied the LMA with equal weights in both domains, thereby obtaining A t and t t as well as A L and t L . Either way of fitting can be done with the pure signal x(t) or with x(t) being convolved with a system response function r(t). Fig. 5a shows the exponential with the true paramters A and t together with the LMA fit of its noisy variant in the t-domain, yielding A t and t t . On the other hand, Fig. 5b shows the Legendre spectrum of the noisy exponential (gray bars) along with the LMA fit in the L-domain (adjacent red bars), yielding A L and t L . In this example we obtain A t~2 999:034 and t t~0 :6025 on the one hand and A L~2 999:761 and t L~0 :6012 on the other. Both pairs approximate the real values A~3000 and t~0:6, whereby, in this case, A L and t L come closer to A and t. Another realization of the same experiment (same A and t) would, of course, lead to a different pair of 'best fits' and the LMA in the t-domain might come closer to A and t this time. We therefore analyzed a large number of realizations (5000) of the same experiment, fitted them using either method, and plotted the distributions of the normalized fitted valuesÂ A t~At =A andÂ A L~AL =A (Fig. 5c) as well ast t t~tt =t andt t L~tL =t (Fig. 5d). Both amplitudes and time constants are estimated with higher precision in the Ldomain, the effect being more pronounced for the time constants. The precision of fitting A and t in the L-domain or the t-domain can be quantified by the Gaussian sum of variances, resulting in the joint errors e L and e t (see Methods, eq. 14a,b). The probability for the error in the L-domain, e L , to be smaller than e t varied as a function of t and was found to lie in the range between 57:7% and 73:4%.
In the above comparison, we have used the LMA with equal weights, w i~1 , for all samples. Thereby we have implicitly assumed a type of measurement, where no a-priori knowledge is available concerning the measurement and its noise. In such cases and in particular when it is unclear whether the data contain one or two exponentials, and one of them might be suppressed by the importance weighting, all weights should be set to 1. Only in cases where reasonable assumptions on measurement errors can be made, which is, of course, the case for our simulated data, an appropriate weighting should be chosen, since this improves the estimation of the parameters. In fact, the comparison of precision and accuracy of the same data as above but with optimal weighting for Poisson noise, shows virtually no difference between the fitting in L-and t-domain ( Fig. 5d and e). However, gradually increasing the gaussian noise renders this way of importance weighting rapidly sub-optimal (Table 1), and the fit in the L-domain gives the accurate parameters more frequently. Fitting in the L-domain thus appears to be equivalent to optimal weighting in the t-domain. This is useful because the optimal weights in the t-domain are mostly not known. R was chosen such that the curve overlaps with x(t) for large t. (f) Legendre spectrum (gray bars) of convoluted noisy exponential shown in e (continuous curve). The lowpass-filtered inverse transform is shown in e (continuous curve) and approximates the convoluted noisy exponential. In addition, f shows the Legendre spectrum of x(t), obtained through eq. 9. The lowpass-filtered inverse transform of this spectrum is shown as the red dashed curve in e and approximates the original nonconvoluted exponential, from which the noisy convoluted curve was generated. doi:10.1371/journal.pone.0090500.g002 Fitting and Filtering in Legendre Space

Optimal number of Legendre components
In the above, we used the first eight components of a Legendre spectrum for calculating the inverse transform ( Fig. 1a,b) and for the comparison of accuracies in the t-and L-domain (Fig. 5c,d). However, the spectral amplitudes' decays in the latter figure suggest that five or six components might have been sufficient to obtain the same result. This and similar observations led us to analyze the question of how many Legendre components should be taken into account for the LMA in the L-domain. The answer turned out to depend on the ratio of t and the time T r over which the record was taken. This can be seen qualitatively by comparing the Taylor series representation of exponentials with the Legendre representation, and quantitatively by calculating the coefficient of variation of Legendre components (Fig. 6). Specifically, the Taylor series representation of an exponential (eq. 9) shows that the terms of the Taylor series decay with n as 1=n!, but that small time constants t counteract this decay, since (1=t n ) increases with decreasing t. As a consequence, for the same accuracy smaller t 0 s  require higher orders of the Taylor series. Clearly, as there is a unique mapping from the Taylor series representation to the Legendre polynomial representation of exponentials (eq. 12), the number of Legendre components required increases for decreasing time constants. As a quantitative way of obtaining an appropriate cut-off index k c of Legendre components we used the coefficient of variation cv(k)~s L k =vL k w, which describes, for each k, the ratio of the standard deviation of L k to its average, i.e., the ratio of noise to signal. For a noisy exponential with A~10 and t~0:9, the Legendre components including their standard deviations are shown in Fig. 6a (gray). Up to the forth component, the standard deviations are rather small and can hardly be recognized. The corresponding cv(k) is plotted in part b of the figure (black). Assuming cv~1 as threshold for an acceptable signal-to-noise ratio, we obtain k c~3 as cut-off. In the case of a shorter time constant t, e.g. 0:3, the standard deviations can virtually not be recognized in the spectrum (Fig. 6a, red) but the cv(k) clearly show that up to k~6 the signal dominates the noise. It can be taken as a rule of thumb, that for t=T r &0:3 six Legendre components give optimal results.

Computational cost
For comparing the computational cost in the t-domain versus the L-domain, we started out with FLIM data from an experiment, where oxidative stress was measured in cultured hippocampal neurons using the HyPer sensor and two-photon microscopy (Fig. 7). For the 150 data samples recorded the computational cost was 8.78 ms/px (t-domain) versus 0.3 ms/px, amounting to approx. 575 s (t-domain) versus 19,7 s (L-domain) for the whole image. Given an image acquisition time of 30 s, the analysis of exponentials in the L-domain has become feasible in real-time.
As the algorithm does not depend on the specific type of experiment, this result can readily be generalized to higher sample sizes. Expectedly, the computational cost of the LMA in the tdomain is approximately proportional to the number of samples. On the other hand, fitting in the L-domain requires not only (i) the fitting in the L-domain but also (ii) the calculation of the finite Legendre transform (fLT), and the rescaling of the parameters (see Methods). The first step is much more time-consuming than (ii) and (iii), but it is independent of the number of samples in the tdomain. We analyzed this quantitatively and carried out both fitting procedures for increasing numbers of experimental samples.
Under our conditions (Intel i7, 2.3 GHz, Python/Numpy), fitting in the L-domain took less than 350ms=fit and was indeed independent of sample size (Fig. 7, black dots), whereas fitting in the t-domain was approximately 12 to 100 times slower for 1024 to 8192 samples, respectively (Fig. 7, black asterisks). In case the recorded function y(t) is a convolution, y(t)~x(t) Ã r(t), fitting in the t-domain takes approximately four times longer (Fig. 7, red  asterisks). This is because a convolution of the model function with r(t) has to be carried out for each minimization step, as the direct deconvolution of y(t) is often not practicable due to the superimposed noise.
On the other hand, fitting the Legendre spectrum of y(t)~x(t) Ã r(t) is almost as efficient as fitting x(t), since this requires only one additional step, namely a matrix multiplication (eq. 8). As this operation involves the samples size, the computational cost increases accordingly, though rather moderately (Fig. 7, red dots). The same applies to the computation of the fLT itself, which took from 45ms (128 samples) to 210ms (8192 samples) for seven components.
Finally, the computational cost depends also on the initial parameter values used. It turned out that calculating the Legendre spectrum of the experimental outcome and solving the first three lines of the system eq. 12 gives a first guess of the three parameters, A,t, and offset, which can be used as starting values, also for the fitting of double exponentials.

Double exponential decay
The analysis of double exponential decays is relevant in many instances but unfortunatety, at least in many biological recordings, the underlying molecular processes are non-stationary so that the experiment cannot be repeated under identical conditions. A typical situation is a decay curve from experimental data which at first glance seems to be a double exponential and which needs to be analyzed. However, when analyzed in t-domain and L-domain, the same data usually give different results for the parameters as shown in Fig. 8a and b. It was therefore indispensable to analyze which of the two methods is superior in the sense that the resulting parameters come closer to the true values. We therefore generated, on the basis of two known amplitudes and time constants, a large number (1000) of double exponentials, added Poisson and gaussian noise as above, and fitted them in both the t-and Ldomain using the LMA in either domain. In addition to having the resulting x 2 values as a measure to judge the quality of fits, we also used the individual errors e t or e L as a measure of how much the individual fitted parameter values deviated from the (known) true ones at minimized x 2 . This way, the fitting result can be observed separately for amplitudes and time constants.
As four parameters A 1 , A 2 , t 1 and t 2 have to be taken into account now (rather than two as above), the corresponding terms for A 2 and t 2 need to be added to eq. 14a and b, respectively. Fig. 8c shows for the relevant range of t 0 s that the double exponential fit in L-domain is more likely to give a better approximation than the corresponding fit in t-domain, and Fig. 8d shows how much the error in the t-domain differs from that in the L-domain. We calculated, for either domain, the success rate p s of the fit. p s was obtained as the frequency, with which the algorithm results in two different amplitudes and time constants, while the correponding failure rate, 1{p s , describes the cases, in which the algorithm converges to only one time constant or where it does not converge at all. The success rate p s for t 1~0 :2 (dashed vertical line in Fig. 8c) is plotted in e. Clearly, for t 2~0 :6, both methods converge to double exponentials in 100% of all trials. With decreasing t 2 , however, the success rate decreases too, the decrease being more pronounced in the t-domain (e, dashed). In Fig. 8c,d we have marked the 50% and the 95% -thresholds of the success rate by continuous or dashed lines, respectively. In these regions, the probability p Lt of more accurate results in the Ldomain varied between 58:8% and 81:8% (p s~5 0%) and between 66:7% and 81:8% (p s~9 5%). Taken together, for the parameter region where both ways of analysis (i.e., t-domain and L-domain) are successful, the fitting results in the L-domain come, on the average, closer to the true double exponential decay parameters. Furthermore, the minimum ratio of t 2 =t 1 that can be differentiated is smaller for the fit in the L-domain (2:2 vs. 2:6).

Discussion
The method proposed herein consists of two steps, first the deterministic transform of the signal into the L-domain, and second, NLLSQ/LMA fitting of the lower Legendre components. The major advantage of the method is the gain in analysis speed, which is primarily brought about by the fact that the LMA is applied to K Legendre components rather than to the N samples in the t-domain, with KvvN. At the same time, the method yields the accurate parameters with high precision, which is essentially brought about by largely separating signal and noise in the L-domain with subsequent parameter estimation from a truncated Legendre spectrum.
The situation may be compared to a spectrally narrow-banded signal buried in noise, e.g., a very faint radar or radio carrier frequency, where the carrier can be detected in Fourier space or by cross-correlation, because the noise power is small that falls into the signal's spectral range. Likewise, exponentials can be well Black and red curves refer to whether (red) or not (black) a system response function was taken into account. Each point or asterisk is the average duration of 100 computations. As the analysis in the timedomain is fastest when using the Fast Fourier Transform, we choose the sample sizes to be powers of 2, starting with 256. Inset, left, FLIM image of a mouse hippocampus neuron probed for oxidative stress using the hydrogen peroxide sensor HyPer [13] and, right, fluorescence lifetime function for one pixel in the middle of the cell (150 samples, data courtesy K. Kizina and M. Mü ller, CNMPB, Gö ttingen). doi:10.1371/journal.pone.0090500.g007 Fitting and Filtering in Legendre Space PLOS ONE | www.plosone.org detected in the L-domain, because the noise power that falls into the lower components, where the signal is mapped, is low.
Legendre polynomials are important functions in many areas of physics [4]. They are orthogonal on the interval [21,1] which allows expansion of any function that is continuous on this interval into its spectral Legendre components. The operation by which the spectrum is obtained has been called Legendre transform (e.g. [5,6]). Though this usage of the term appears straightforward, it is ambiguous, since in the canonical language of physics and chemistry, Legendre transformations are well-known and widely used to express a function z(t) in terms of its derivative rather than in terms of its independent variable t [7]. We therefore follow Jerri [4] and Méndez-Pérez and Morales [6], and name the operation by which the Legendre spectrum is obtained (eq. 4) finite Legendre transform (fLT), in close analogy to the finite Fourier transform, which is also carried out over a finite interval.
The Legendre spectra of a noisy exponential and its noise-free variant differ in the amplitudes of their components. While the higher components vary considerably with the noise in the time function, the lower components contain little noise, so that backtransforming them into the t-domain results in effective noise removal. In analogy to conventional lowpass filtering in the Fourier domain [8], this filter should be called Legendre lowpass. Evidently, not only exponentials can be lowpass filtered this way. We have, for instance, successfully checked this for the product of two exponentials (Fig. 4), for bleaching time courses in LSM imaging (not shown), and for fluorescence correlation spectroscopy data, where, in the simplest case, the autocorrelation function takes on the form g(t)~(N(1zt=t diff )) {1 , N and t diff being constants (Fig. 1c). However, as this paper is concerned with exponentials, we did not investigate which other classes of functions can advantageously be Legendre lowpass filtered.
Apart from Legendre filtering, the lower Legendre components can also be used to obtain the exponentials' parameters by finding the Legendre spectrum that fits best, in the least squares sense, the fLT of the data. In practice this means that the LMA is applied to the lower Legendre components of an experimental record (or our simulation) using the correspondigly truncated and parameterdependent Legendre spectrum of a noiseless exponential as the model. Prior to the actual fitting, an appropriate record length T r as well as a good guess of the A 0 s and t 0 s involved must be determined. The cases t%T r and t&T r can be excluded for analysis, because in either regime the samples considered carry too little infomation about the exponential. Satisfactory fitting results are obtained in the range 0:1vt=T r v1, and the ratio t=T r^0 :3 turned out to be the best choice. In a next step, we find a good guess for the initial values A and t by comparing the Legendre representation of an exponential with its Taylor series. Both representations involve the powers of time, t 0 , t 1 , …, and the comparison of the respective coefficients leads to a system of linear equations, which allows calculating A and t from the Legendre spectrum of the signal and the constant coefficients of the matrix P of Legendre polynomials. Due to the noise in the signal's Legendre components and the finite size of the system (eq. 12), these values for A and t are, of course, mere approximations, but as such they are optimal starting values for the fitting procedure. (Obviously, the same initial values A and t could also be used for least squares fitting in the t-domain.) Fitting in the L-domain is much faster than in the t-domain, particularly if the exponentials under investigation are convolved with a system response function (Fig. 7). The computaional costs in time and L-domain correspond, in a good approximation, to O(N t ) and O(N L ), with N t and N L being the number of values to be fitted, respectively. This is highly relevant for time-sensitive  . (a,b) The interaction of Ru(II) complexes with DNA (same as in in Fig. 1) shows double exponential decay( [12], data courtesy of F. Secco, Univ Pisa,I). Same data fitted in t-(a) and L-domain (b). (c) Probability for the fitting error in the L-domain, e L , to be smaller than that in the t-domain, e t , represented as a function of t 1 and t 2 . Color code of the probability is shown beneath the plot. For each pair (t 1 ,t 2 ) 1000 trials were computed. (d) Relative difference of fitting errors, (e t {e L )=e t , as a function of t 1 and t 2 . Color code beneath the plot. 1000 trials per pixel. Left of the solid and dashed white lines in a and b, the success rate of the fit in the L-domain is larger that 50% and 95%, respectively. Left of the solid and dashed black lines in c and d, the success rate of the fit in the t-domain is larger that 50% and 95%, respectively. Relative error differences were calculated only for successful trials. In addition to being faster, the fits in the L-domain have the same or a higher probability of giving the accurate paramters A and t. Without importance weighting of the data, fitting in the Ldomain is always better (Fig. 5). Accordingly, fitting in the the Ldomain is also superior to earlier algorithms which either do not converge as well as the LMA or are less precise than the LMA [9].
Regarding the comparisons of precision, a note of caution is necessary. Whenever assumptions on measurement errors can be made, an appropriate weighting should be chosen, since this improves the estimation of the parameters. However, in cases with little or no a-priori knowledge on measurement errors, optimal weights cannot be found, and fitting in the L-domain is more precise.
When analyzing double exponential decays, the probability for obtaining better parameter approximations is always higher in the L-domain. In addition, there is a range of the parameter space where fitting in the L-domain still gives satisfactory estimates, while the analysis in the t-domain fails to converge.
At the beginning of our study we also considered alternative filtering methods such the wavelet transform, Chebychev or Laguerre polynomials. However, the wavelet transform is slower than the fLT and also more appropriate for a different kind of signal, being defined in a given time and frequency window, while the Chebychev and (associated) Laguerre polynomials, albeit similar to the Legendre polynomials, are only orthogonal, when a weighting function applying to the inner product, which makes their use computationally less efficient.
Experimental outcomes often show an offset, which has to be taken into account when using the LMA in the t-domain. In Legendre space, the offset is simply added to the first component, i.e., the first component of the Legendre fit gives the constant amplitude plus the offset. A good guess for the offset can be calculated from eq. 10, where three equations, which have to include the first one, need to be solved to obtain estimates for the offset as well as for A and t (the first equation then results in A plus the offset).
Historically, Legendre appears to play a double role in exponential fitting, since, apart from the fast parameter retrieval described herein, he was the first to publish the idea of least squares fitting in 1806 (though this method is mostly attributed to Gauss, who claimed the first usage of it) [10].

Simulated data and rationale of their usage
To simulate a set of noisy exponentials having the same amplitude and time constant, we first defined the amplitude A and the time constant t for this set. We then calculated the exponential x(t) and added Poisson noise to it. Where indicated, we also added an offset and gaussian noise. We name the (noiseless) exponential x(t) and the trials, which differ only in their noise, x(t). By definition, we thus have a-priori knowledge of A, t, A o , and x(t). The sample size N was chosen to be 1000, except for the comparison of computational costs.
Simulated data and thus the a-priori knowledge of A and t are necessary, since we need to compare the fitting results in both the t-and the L-domain with the real values of the parameters A true and t true . On the one hand, the direct application of the LMA to a noisy exponential x(t) (ie, LMA in the t-domain) yields specific values A t and t t , on the other, transforming the exponentials into the L-domain and applying the LMA to the lower components of the Legendre spectrum gives values A L and t L . Thus, for every generated noisy curve we get four parameters A t , t t , A L and t L . The ratios A i /A true , t i /t true , i meaning t or L, tell by which factor the fitted values deviate from the real parameters. A ratio of 1, for instance, indicates a perfect fitting result. Finally, plotting the probability density functions of the deviations reveals accuracy and precision of the fits in the t-and the L-domain (Fig. 2).

Legendre polynomials and finite Legendre transform (fLT)
Legendre polynomials result from the orthogonalization of the powers 1, t, t 2 , … leading to with P 0 (t)~1. From this we have a recursive definition with P 0 (t)~1 and P 1 (t)~t. The first five Legendre polynomials are given by and plotted in the Fig. S1. Legendre polynomials are an orthogonal set of functions on the interval [21,1], so that any function x(t) defined and continuous on this interval can be transformed into its Legendre spectrum by using the finite Legendre transform (fLT), In our case, t is a discrete variable assuming N values between 21 and 1. The factor (2kz1)=N is a normalization factor.
The time function x(t) can be regained from its Legendre spectrum L x (k) by the inverse fLT (ifLT), Alternatively, one might consider to use the shifted Legendre polynomials [11], which are defined as P Ã k~P k : (2x{1) and orthogonal on the interval [0,1]. While this leads to the same results, it is computationally inconvenient, because the factor (2x{1) needs to be taken into account in a number of operations.

Legendre spectrum of exponentials
Let us consider a class of stochastic processes x(t) characterized by the parameters A and t as well as by non-stationary means g x (t)~A : e t=t , and Poisson-noise fluctuations.
We conveniently represent the respective experimental outcomes, i.e., the realizations x(t) of such a process, in the space of Legendre polynomials P n . While the fLT L x of x(t) is a sequence of random variables L x (k)~fL x (0),:::,L x (k),:::g, the fLT of any particular realization x(t) of x(t) gives a specific Legendre spectrum, L x~f L x (0),:::,L x (k),:::g, i.e., the coefficients of a linear combination of Legendre polynomials.
The first practical step of data analysis is to map the time inverval of the experimental results onto the interval [21,1], i.e., we redefine the time axis so that the exponentials are spanned over [21,1] rather than over ½0,T r , T r being the time of the last sample.
In a second step we calculate the components L k (x) of the the Legendre spectrum (eq. 4).
In cases, where an experimental output y(t) is the convolution of an exponential x(t) and a device response r(t), y(t)~x(t) Ã r(t), we can express y(t) as whereby the L i are the unknown Legendre spectrum components of the noisy exponential to be found. This equation can be rewritten in terms of a matrix P T r consisting of the transposed Legendre polynomials each convoluted with r(t) with L being the vector of (unknown) Legendre coefficients of the pure, non-convoluted exponential x(t) recorded. Using the pseudoinverse of P T 0 r , we obtain The pseudoinverse and the resulting Legendre components L of the experimental outcome y(t) are calculated once only, prior to fitting (A,t) -dependent Legendre components to the experimental spectrum. With increasing sample size the calculation of L increases accordingly (Fig. 7).
From the inverse transform of the vector L of Legendre coefficients, an approximation of the non-convoluted exponential can be obtained.

Direct retrieval of amplitudes and time constants from Legendre components
The finding that the first spectral Legendre components are virtually identical for noisy and pure exponentials with the same parameters suggests a third and analytical way of obtaining A and t. We start out with the Taylor series of an exponential, and represent the same exponential by the superposition of Legendre polynomials Ae {t=t~L (0)P 0 (t)zL(1)P 1 (t)zL(2)P 2 (t)z:::zL(6)P 6 (t)z::: After rearranging with respect to the powers of t, i.e., where the index indicates the cut-off, P stands for the matrix of Legendre polynomial coefficients, L for the truncated Legendre component vector, and T for the vector of Taylor series coefficients. In case the exponential tends to an offset value A o , x(t)?A o for t??, A o adds to the first term of the Taylor series and, consequently, to L(0) Taken together, the analytically derived equation 12 relates, for a given order, the Legendre components of an exponential to its paramters A and t. This is useful because it gives a good guess for A and t from the lower Legendre components.

Fitting errors
As every fitted exponential x(t) depends on A and t, x(A,t)~Ae {ti=t , the resulting error of each value x i can easily be obtained from Gauss' error propagation, i.e., The total errors of a fit in the L-or the t-domain, e L or e t , are thus Scaling If a function is to be fitted in Legendre space, it has to be scaled to the time interval [21,1], where A o is an offset of the exponential. Prior to fitting this does not require any computation as we may assume the function's first and last value to occur at t~{1 and at t~1 rather than at t~0 and t~T r . After fitting, the time scaling factor 2=T r is needed to obtain the real time constants t r . Likewise, the amplitude needs to be rescaled by the factor e 1=t , since, on the interval [21,1], the real amplitude A r is the value at t~{1, while the fit gives the amplitude value A at the scaled time t~0. The exponential in the t-domain, x(t r )~A o zA r e {tr=tr ,0ƒtƒT r can thus be obtained by

Software
All calculations were carried out in Python/numpy/C. The code used in this paper, in particluar filtering routines, fitting routines, the Legendre polynomials, the fLT, ifLT, are available at the public ''python package index'' repository (https://pypi.python.og/pypi), and on our website (https://www.ukmn.gwdg.de/). (a) Exponential of known amplitude and time constant. b) Legendre spectrum of (a). c) Poisson noise corresponding to the amplitude of the exponential in (a) and d) its Legendre transform. e) Stationary gaussian noise and f) its Legendre transform. g) Noisy exponential, superposition of (a), (c), and (e), the latter showing gaussian noise with mean A o~5 and standard deviation s~1:5. The mean simulates the offset of the measurent. h) Legendre transform of (g). For the comparison of methods, the LMA is applied in both domains, and the respective results are compared to the known parameters of the simulated data. A t , t t , A L , and t L are estimates of the true values A true , t true , which, in the following, we call for short A and t. (EPS)