[SciPy-User] fitting with convolution?

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] fitting with convolution?

x.piter
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
Reply | Threaded
Open this post in threaded view
|

Re: fitting with convolution?

David Baddeley
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
Reply | Threaded
Open this post in threaded view
|

Re: fitting with convolution?

Till Stensitzki
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