傅里叶变换 / 使用numpy/scipy在Python中进行迭代反卷积拟合

4
我希望适配荧光寿命曲线。这些曲线是仪器响应函数(IRF, 假设高斯分布)和(多)指数衰减的卷积结果。

formula

G是高斯函数,F是指数衰减函数。我尝试使用lmfit.minimize在Python中拟合此函数,使用以下函数:

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  

其中,gaus定义了高斯仪器响应函数,exp1是单指数衰减。func1使用scipy.integrate计算两者之间的卷积,它返回的值用于在一组参数给定的情况下计算拟合值和数据之间的差异:

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))

虽然这种方法可以工作,但拟合过程非常缓慢,并且尚未给出任何合理的拟合结果。
有几种方法可以规避卷积过程,例如使用傅里叶变换或迭代去卷积。Origin似乎知道如何做到这一点,但我很难理解这些步骤。
据我所知,迭代去卷积的原理是先将信号与猜测的高斯函数进行去卷积,然后拟合结果,最后调整高斯函数。
傅里叶变换的方法基于实空间中的卷积对应于傅里叶域中的乘法的原理,这将减少计算时间。我猜它会对信号进行傅里叶变换,拟合结果,然后再进行傅里叶反变换。
我想了解如何在Python numpy/scipy中实现这些方法。迭代去卷积似乎最容易实现,所以也许我应该从那里开始。另一方面,根据我的阅读,傅里叶方法应该更快,更可靠。然而,我不知道如何对傅里叶变换后的结果进行拟合。
1个回答

3

一些评论和想法:

1)在您上面的实现中,gtoff和expoff是多余的。同样适用于gamp和expamp。例如使用expoff = 0和expamp = 1。

2)如果您只关心指数衰减与高斯内核的卷积,请使用解析结果,可以使用scipy.special.erf轻松高效地计算:指数衰减exp(-t / tau)(t≥0)与归一化高斯函数1 /(sqrt(2 * pi)* sigma)* exp(-t ^ 2/(2 * sigma ^ 2))的卷积为

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

这显然可以推广到您对exp1和gaus的定义。

3)您可以使用np.convolve有效地卷积向量,其基于您帖子中提到的FFT方法。对于您特定的应用程序,高斯函数应以包裹顺序表示。有关详细信息,请参见例如Numerical Recipes,第13.1章。

4)传统方法基于迭代重新卷积,使用通常恒定的高斯内核。反卷积虽然可能,但是不适用于所有情况。

5)已经开发出一种适用于拟合时间分辨荧光数据的通用软件包,其中实现了您需要的所有内容:trfit


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接