[SciPy-User] fitting with convolution?

David Baddeley david_baddeley@yahoo.com...
Fri Jul 19 08:57:32 CDT 2013


Hi Petro,

when I run the code I get very similar, although not identical curves, which are roughly what you'd expect. Your exponential kernel that you're convolving with is not band limited, but FFTs (and hence FFT based convolution) are band-limited to whatever your effective nyquist frequency is (strictly speaking out of band frequencies might be aliased into the pass band - but in your case you can just think of the high frequency components as being discarded). This probably applies generally to any discrete way of doing convolution, when compared to an analytical solution. Put simply, you can never sample fine enough  to reproduce all the detail of the initial spikey bit of your exponential (if you increase your sampling frequency you'll get a better match, but you'll never quite get there). If you are referring to the factor of 14.5, it is a result of fft normalisation, and should more accurately be sqrt(N)/pi.

hope this helps,
David


________________________________
 From: Petro <x.piter@gmail.com>
To: scipy-user@scipy.org 
Sent: Friday, 12 July 2013 2:13 PM
Subject: [SciPy-User] fitting with convolution?
 

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

_______________________________________________
SciPy-User mailing list
SciPy-User@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20130719/33048e09/attachment.html 


More information about the SciPy-User mailing list