[SciPy-User] fitting with convolution?
Petro
x.piter at gmail.com
Fri Jul 12 14:13:47 EDT 2013
Hi all,
I try to fir a time-resolved dataset with multiple exponents convoluted
with a Gaussian instrument response function (IRF).
I had a look how it is done in Origin
http://wiki.originlab.com/~originla/howto/index.php?title=Tutorial:Fitting_With_Convolution
There fft_fft_convolution calculates the circular convolution of an
exponent with IRF.
I have found a similar function for python here:
http://stackoverflow.com/questions/6855169/convolution-computations-in-numpy-scipy
This convolution also can be calculated analytically as, for example, in
this package:
http://www.photonfactory.auckland.ac.nz/uoa/home/photon-factory/pytra
def convolutedexp(tau,mu,fwhm,x):
d = (fwhm/(2*sqrt(2*log(2))))
return 0.5*exp(-x/tau)*exp((mu+(d**2.)/(2.*tau))/tau)*
(1.+erf((x-(mu+(d**2.)/tau))/(sqrt(2.)*d)))
def gaussian(mu,fwhm,x):
d = (fwhm/(2.*sqrt(2.*log(2.))))
return exp(-((x-mu)**2.)/(2.*d**2.))
My problem is if I compare analytical and circular convolution they do
not match:
_____source_________
import numpy
from scipy.special import erf
def cconv(a, b):
'''
Computes the circular convolution of the (real-valued) vectors a and b.
'''
return fft.ifft(fft.fft(a) * fft.fft(b)).real
def convolutedexp(tau,mu,fwhm,x):
d = (fwhm/(2*sqrt(2*log(2))))
return
0.5*exp(-x/tau)*exp((mu+(d**2.)/(2.*tau))/tau)*(1.+erf((x-(mu+(d**2.)/tau))/(sqrt(2.)*d)))
def gaussian(mu,fwhm,x):
d = (fwhm/(2.*sqrt(2.*log(2.))))
return exp(-((x-mu)**2.)/(2.*d**2.))
t = array(linspace(-10.0,1000.0,2040.0))[:-1]
mu = 0
fwhm = 4.0
tau = 20.0
uf = gaussian(mu,fwhm,t)
vf = exp(-t/tau)
figure(figsize=[12,12])
plot(t,uf)
#plot(t,vf)
uvf1 = cconv(uf,vf)
plot(tuv,uvf1/14.5)
uvf2 = convolutedexp(tau,mu,fwhm,t)
plot(t,uvf2)
xlim([-10,20])
____source_end___
My feeling is that I miss something about convolution?
Can anybody give me a hint?
Thanks.
Petro
More information about the SciPy-User
mailing list