Gaussian line profile using numeric python
Author 
Message 
Geoff Lo #1 / 9

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((vvo)^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 


Konrad Hinse #2 / 9

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((vvo)^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.: +332.38.25.56.24 Rue Charles Sadron  Fax: +332.38.63.15.17 45071 Orleans Cedex 2  Deutsch/Esperanto/English/ France  Nederlands/Francais 

Mon, 15 Sep 2003 19:45:02 GMT 


Tim Hochber #3 / 9

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 


John J. Le #4 / 9

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((vvo)^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 


Konrad Hinse #5 / 9

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 zeropadding in case of nonperiodic 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.: +332.38.25.56.24 Rue Charles Sadron  Fax: +332.38.63.15.17 45071 Orleans Cedex 2  Deutsch/Esperanto/English/ France  Nederlands/Francais 

Tue, 16 Sep 2003 18:23:39 GMT 


John J. Le #6 / 9

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 zeropadding in case of nonperiodic 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 nonperiodic functions (as I didn't, to start with!). John

Wed, 17 Sep 2003 22:04:51 GMT 


Robert Ame #7 / 9

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 zeropadding in case of nonperiodic 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 


John J. Le #8 / 9

Gaussian line profile using numeric python
Quote:
[...] > >Except for zeropadding in case of nonperiodic 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 nonperiodic 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 


Robert Ame #9 / 9

Gaussian line profile using numeric python
Quote:
>> [Convolution using FFT] >He was just pointing out that it doesn't make sense to take a >nonperiodic 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 


