Fourier transform / iterative deconvolution fitting in python using numpy/scipy

Question

I would like to fit a fluorescent lifetime curve. These curves are given by the convolution of the instrument response function (IRF, assumed Gaussian) and a (multi) exponential decay:

$\int_{0}^{t} G(t - T)F(T) dT$

Where G is the gaussian and F exponential decay. I have tried fitting this in python using lmfit.minimize, with the following functions:

def gaus(t, gamp, gwidth, gtoff):
    i = gamp*math.exp(-(t-gtoff)**2/gwidth)
    return i


def exp1(t, expamp, exptime, exptoff):
    i = expamp*math.exp((t-exptoff)/(exptime))
    return i

def func1(t, gamp, gwidth, gtoff, expamp1, exptime1, exptoff1):

    i = [(scipy.integrate.quad(lambda Tpr: gaus((ti - Tpr), gamp, gwidth, gtoff)*exp1(ti, expamp1, exptime1, exptoff1), 0, ti))[0]
            for ti in t]
    return i

Where gaus defines the Gaussian instrument response function, exp1 is single exponential decay. func1 calculates the convolution between the two with scipy.integrate, the values it returns are used to calculate the difference between the fit and the data given a set of parameters:

def fitfunc(params, x, data):
    ..getting parameters part here..

    values = func1(t, gamp, gwidth, toff, expamp1, exptime1, toff)
    return values - data

result = lmfit.minimize(fitfunc, fit_params, args = (t, test))

Although this sort of works, the fitting process is very slow and has not yet given my any reasonable fits.

There are several approaches to circumvent the convolution process by using Fourier transforms or iterative deconvolution. Origin seems to know how to do it, but I am having trouble understanding the procedures.

As far as I can tell, iterative deconvolution works by deconvolving the signal with a guessed Gaussian function and then fitting the result, and then adjusting the Gaussian.

The Fourier transform approach would be based on the principle that convolution in real space corresponds to a multiplication in the Fourier domain, which would reduce computation time. I'm guessing that it would be Fourier transforming the signal, fitting that, and then Fourier transforming back.

I'd like some input on how to implement these methods in python numpy/scipy. The iterative deconvolution seems to be the easiest to do, so maybe I should start with that. On the other head, from what I've read, the Fourier method should be faster and more reliable. However, I have no idea how to perform fitting on the Fourier transformed result.

hth hth · Accepted Answer · 2013-05-10T20:38:35

Some comments and thoughts:

1) In your above implementation gtoff and expoff are redundant. The same applies for gamp and expamp. Use, e.g., expoff=0 and expamp=1.

2) If you are only after the convolution of multi-exponential decays with Gaussian kernels, use the analytic result, which is easily and efficiently evaluated using scipy.special.erf: The convolution of exp(-t/tau) (t >= 0) with the normalized Gaussian 1/(sqrt(2*pi)*sigma)*exp(-t^2/(2*sigma^2)) is

1/2 * exp( - (2*t*tau - s2)/(2*tau^2) ) * (1 + erf( (2*t*tau - s2)/(sqrt(2)*tau*sigma) )).

This obviously generalizes to your definitions of exp1 and gaus.

3) You can efficiently convolve vectors using np.convolve, which is based on the FFT approach mentioned in your post. For your particular application, the Gaussian ought to be represented in wrapped-around order. See e.g. Numerical Recipes, Chapter 13.1 for details (otherwise your time origin will shift).

4) The conventional approach is based on the iterative re-convolution using a typically constant Gaussian kernel. The deconvolution, though possible, is ill-posed.

5) A versatile package for fitting time resolved fluorescence data has been devised that implements all what you require: trfit

Fourier transform / iterative deconvolution fitting in python using numpy/scipy

1 Answers