Gaussian line profile using numeric python
Author Message
Gaussian line profile using numeric python

Hi.

I have an extensive set of xy data that I want to convolute with a gaussian
profile.  Below is my algorithm,  which doesn't work.  Can anyone see where
I'm being dumb?

import Numeric
from MLab import sum
from Numeric import arange,zeros,transpose,shape,exp,power,pi
from math import sqrt

def get_data(filename):
x = d[:,0]
y = d[:,1]
del d
return x,y

def main(inpfile,outfile='res.out',sfreq=350,efreq=5050,stepf=2,fwhm=2.5):
fin,iin = get_data(inpfile)
xpt = arange(sfreq,efreq,stepf,'d')
ypt = calcspec(fin,iin,xpt,fwhm)
del fin,iin
put_data(outfile,xpt,ypt)

def put_data(filename,x,y):
tmp = TableFromColumns([x,y])
tmp['filename'] = filename
writeTable(tmp)
del tmp

def calcspec(fin,iin,x,fw):
norm,expterm = 0.0,0.0
inten = zeros((len(x)),'d')
norm = float(1/(sqrt(2.0*pi)*fw))
expterm = float(1/(2*(fw**2)))
for i in range(len(x)):
print "Starting step ",i+1,"of ",len(x)
print "sum fin: ",sum(fin)," sum iin: ",sum(iin)
print "x[",i,"] = ",x[i],"; SUM(x[",i,"] - fin) = ",sum((x[i]-fin[:]),0)
#inten[i] = sum((norm*iin[:]*exp(-1.0*(power((expterm*(xpts[i]-fin[:])),2)))),0)
inten[i] = norm*sum(iin[:]*exp(-1.0*(power(((x[i]-fin[:])*expterm),2))))
#for j in range(len(fin)):
#inten[i] = inten[i] + float(norm*iin[j]*exp(-1.0*((expterm*(xpts[i]-fin[j]))**2)))
return inten

P.S.  I'm no math wizard but this is my interpretation of the gaussian lineshape
function
G(v) = 1/(sqrt(2*PI)*FWHM)*exp(-(v-vo)^2/(2*FWHM^2))

Cheers
Geoff

http://www.*-*-*.com/

######
Fortune favours the brave....
but keeps an axe aside for the stupid!
######

Mon, 15 Sep 2003 19:03:58 GMT
Gaussian line profile using numeric python

Quote:

> I have an extensive set of xy data that I want to convolute with a gaussian
> profile.  Below is my algorithm,  which doesn't work.  Can anyone see where
> I'm being dumb?

To be honest, I am too lazy to look through this. However,

1) The function Numeric.convolve is certainly a lot more efficient
than your code, and should make it simpler.

2) If you work with large data sets, you ought to use FFTs to
compute the convolution, that's O(N*log(N)) in the size of the
data set, instead of O(N**2) for the straightforward method.
Any good book on FFTs should explain how this works in detail.

Quote:
> P.S.  I'm no math wizard but this is my interpretation of the gaussian lineshape
> function
> G(v) = 1/(sqrt(2*PI)*FWHM)*exp(-(v-vo)^2/(2*FWHM^2))

Assuming that FWHM stands for full width at half minimum, no, it's not
correct. What you call FWHM is the variance sigma of the Gaussian, the
FWHM is around 2.35*sigma (quick calculation, no guarantee).
--
-------------------------------------------------------------------------------

Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

Mon, 15 Sep 2003 19:45:02 GMT
Gaussian line profile using numeric python

Like Konrad, I'm too lazy to completely figure this out, but I do have some

Geoff Low is ">>":

Quote:

> > I have an extensive set of xy data that I want to convolute with a
gaussian
> > profile.  Below is my algorithm,  which doesn't work.  Can anyone see
where
> > I'm being dumb?

> To be honest, I am too lazy to look through this. However,

Ditto. However some thoughts / questions. Is your data evenly spaced in x?
If it is you might want to pass in deltaX instead of x to make that clear.
If it's not, you're going to need a considerably more sophisticated
algorithm. If it is, go with Numeric.convolve as Konrad suggests.

Quote:
> 2) If you work with large data sets, you ought to use FFTs to
>    compute the convolution, that's O(N*log(N)) in the size of the
>    data set, instead of O(N**2) for the straightforward method.
>    Any good book on FFTs should explain how this works in detail.

One thing to watch for is that convolution using FFT makes an implicit
assumption that your data is periodic. If it's not, and you want to try this
route, look up zero padding  in that FFT book that I'm _sure_ you'll consult
before you try this <0.5 wink>.

-tim

Mon, 15 Sep 2003 23:45:12 GMT
Gaussian line profile using numeric python

Quote:

> > I have an extensive set of xy data that I want to convolute with a gaussian
> > profile.  Below is my algorithm,  which doesn't work.  Can anyone see where
> > I'm being dumb?

What goes wrong?  Couldn't see your original message.

Quote:
> To be honest, I am too lazy to look through this. However,

Me too.  :)  I do have some working routines that will do gaussian
convolution, using the standard Numeric FFT module -- very simple.  Or at
least would be if I didn't have a lot of stuff for using the FFTW library
in there, which doesn't actually work (the docs for the rfftw part of the
FFTW python wrapper seem to be wrong:  anybody got this working?).

Quote:
> 1) The function Numeric.convolve is certainly a lot more efficient
>    than your code, and should make it simpler.

I seem to recall that it doesn't cover all the possibilities, though.  I
don't remember exactly why, but since I wrote my own functions to do this
there must have been a reason!

Quote:
> 2) If you work with large data sets, you ought to use FFTs to
>    compute the convolution, that's O(N*log(N)) in the size of the
>    data set, instead of O(N**2) for the straightforward method.
>    Any good book on FFTs should explain how this works in detail.

but you probably don't need to know in detail, it boils down to

inv_FFT(FFT(data)*FFT(kernel))

But if you have small data sets it might be slower, of course.

Quote:
> > P.S.  I'm no math wizard but this is my interpretation of the gaussian lineshape
> > function
> > G(v) = 1/(sqrt(2*PI)*FWHM)*exp(-(v-vo)^2/(2*FWHM^2))

> Assuming that FWHM stands for full width at half minimum, no, it's not
> correct. What you call FWHM is the variance sigma of the Gaussian, the
> FWHM is around 2.35*sigma (quick calculation, no guarantee).

That number is about right, from memory.

John

Tue, 16 Sep 2003 02:01:44 GMT
Gaussian line profile using numeric python

Quote:
> > 2) If you work with large data sets, you ought to use FFTs to
> >    compute the convolution, that's O(N*log(N)) in the size of the
> >    data set, instead of O(N**2) for the straightforward method.
> >    Any good book on FFTs should explain how this works in detail.

> but you probably don't need to know in detail, it boils down to

> inv_FFT(FFT(data)*FFT(kernel))

Except for zero-padding in case of non-periodic signals. That's why
I always recommend to look at a detailed description before doing
this; some of my colleagues keep telling me that FFT techniques "don't
work", and I suspect this is the reason.
--
-------------------------------------------------------------------------------

Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24
Rue Charles Sadron                       | Fax:  +33-2.38.63.15.17
45071 Orleans Cedex 2                    | Deutsch/Esperanto/English/
France                                   | Nederlands/Francais
-------------------------------------------------------------------------------

Tue, 16 Sep 2003 18:23:39 GMT
Gaussian line profile using numeric python

Quote:

[...]
> > but you probably don't need to know in detail, it boils down to

> > inv_FFT(FFT(data)*FFT(kernel))

> Except for zero-padding in case of non-periodic signals. That's why
> I always recommend to look at a detailed description before doing
> this; some of my colleagues keep telling me that FFT techniques "don't
> work", and I suspect this is the reason.

:)

Yes.  And that's (IIRC) what Numeric's convolve() doesn't do quite as
generally as it might -- but there is (IIRC again) an XXX in the manual

Good old 'Numerical Recipes' has a reasonable description of the
accountancy for the various kinds of convolutions, as long as you've
grasped that you do inevitably lose some of your function to end effects

John

Wed, 17 Sep 2003 22:04:51 GMT
Gaussian line profile using numeric python

Quote:

>> > 2) If you work with large data sets, you ought to use FFTs to
>> >    compute the convolution, that's O(N*log(N)) in the size of the
>> >    data set, instead of O(N**2) for the straightforward method.
>> >    Any good book on FFTs should explain how this works in detail.

>> but you probably don't need to know in detail, it boils down to

>> inv_FFT(FFT(data)*FFT(kernel))

>Except for zero-padding in case of non-periodic signals. That's why
>I always recommend to look at a detailed description before doing
>this; some of my colleagues keep telling me that FFT techniques "don't
>work", and I suspect this is the reason.

Well, why dont you/they try a different kind of integral transform,
that might give a better result. Perhaps something like wavelets? Of
course, a fast algorithm must be available to do the transform and its
inverse but if there is it might be worth to do some experimenting.

Robert Amesz

Thu, 18 Sep 2003 06:29:20 GMT
Gaussian line profile using numeric python

Quote:

[...]
> >Except for zero-padding in case of non-periodic signals. That's why
> >I always recommend to look at a detailed description before doing
> >this; some of my colleagues keep telling me that FFT techniques "don't
> >work", and I suspect this is the reason.

> Well, why dont you/they try a different kind of integral transform,
> that might give a better result. Perhaps something like wavelets? Of
> course, a fast algorithm must be available to do the transform and its
> inverse but if there is it might be worth to do some experimenting.

[...]

He was just pointing out that it doesn't make sense to take a non-periodic
function, smear it out, and then expect there to be no end effects.  With
FFT convolution, for example, you end up smearing from 'the other end' of
your function, whether you like it or not.  You have to reckon on
calculating more of your function that you want to end up with, so you can
throw the wrongly smeared regions out.

Err, hope that's clear...  As Konrad says, this is explained properly in
many Numerical Analysis books.

John

Fri, 19 Sep 2003 22:27:55 GMT
Gaussian line profile using numeric python

Quote:

>> [Convolution using FFT]

>He was just pointing out that it doesn't make sense to take a
>non-periodic function, smear it out, and then expect there to be no
>end effects.  With FFT convolution, for example, you end up smearing
>from 'the other end' of your function, whether you like it or not.
>You have to reckon on calculating more of your function that you
>want to end up with, so you can throw the wrongly smeared regions
>out.

Oh, I know, but the smear from other types of fast integral transforms
might not be as objectionable (i.e. broad) as that of the FFT, so you
might not need as much empty guard space to prevent it from folding
back. It was just a thought, untried and untested.

Robert Amesz

Sat, 20 Sep 2003 10:37:27 GMT

 Page 1 of 1 [ 9 post ]

Relevant Pages