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 [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
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 <[hidden email]> To: [hidden email] 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 [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
David Baddeley <david_baddeley <at> yahoo.com.au> writes:
> > > 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 I this case the problem is much simpler, the analytical solution is not circular, while the FFT convolution is. With the right zero-padding it is possible to get identical results. But my testing did show this is slower than just using the analytical solution. greetings Till _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |