[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