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 TableIO import readTableAsArray,TableFromColumns,writeTable
from math import sqrt

def get_data(filename):
    d = readTableAsArray(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
comments on it below.

Konrad Hinsen is ">"
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
somewhere that admits it.

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
for non-periodic functions (as I didn't, to start with!).

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  
 
 [ 9 post ] 

 Relevant Pages 

1. Quadrature.py - Numerical Integration using Gaussian Quadrature

2. Random Generator using Gaussian distribution!

3. Gaussian Elimination using Prolog ?

4. Numeric Python for Python 1.5

5. Code metrics / line counts using python

6. How to do per-line profiling?

7. Help on script to get numeric difference of alternate lines of a file

8. NumPy v7: syntax error at line 171 of Numeric.py

9. Using single puts line spanning serveral source lines

10. ANNOUNCE: Profile 1.0 (Simple Tcl Profiling)

11. Products Created Using Clarion - Another Profile Exchange

12. Profiling using VxWorks and AdaMULTI.

 

 
Powered by phpBB® Forum Software