Fitting a convolution of gaussian and exponential

Share icon

Im new to Igor, and need to fit some data with a convolution of a gaussian and multiexponential function. The code i've written returns a funcfiterror: "the fitting function returned NaN for at least one X value". Here's a copy of the code, your help is much appreciated!! Thanks

Function convfunc(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0]+(x>=0)*(1-exp(-xw/pw[4]))*(pw[1]*exp(-xw/pw[5])+pw[2]*exp(-xw/pw[6]))
// pw[0] = c0
// pw[1] = c1
// pw[2] = c2
// pw[3] = tg
// pw[4] = tau
// pw[5] = tau1
// pw[6] = tau2

Variable resolutionFactor = 10
Variable dT = min(pw[3]/resolutionFactor, pw[4]/resolutionFactor)
Variable nExpWavePnts = round((resolutionFactor*pw[6])/dT)*2 + 1
Make/D/O/N=(nExpWavePnts) expwave
Variable nYPnts = max(resolutionFactor*numpnts(yw), nExpWavePnts)
Make/D/O/N=(nYPnts) yWave
setscale/P x -dT*(nExpWavePnts/2),dT,expwave
setscale/P x xw[0],dT, yWave
expwave = pw[0]+(x>=0)*(1-exp(-xw/pw[4]))*(pw[1]*exp(-xw/pw[5])+pw[2]*exp(-xw/pw[6]))

ywave = exp((-xw/pw[3])^2)
Variable sumywave
sumywave = sum(ywave, -inf,inf)
ywave /= sumywave
convolve/A expwave, ywave
yw = yWave(xw[p])
End testfit.pxp (212.3 KB)
Thanks for the reply John.

I noticed a mistake in the ywave gaussian formula (minus sign was inside parenthesis) but i fixed that. Now there is a different error though...

new error: "while executing Make, the following error occurred: the number of pints in a wave must be between 0 and 2147 million"

Here is what pops up in the command line:

•FuncFit/H="1000000"/NTHR=0 convfunc wave2 wave0 /X=wave1 /D /C=T_Constraints
Fit converged properly
--Curve fit with constraints--
No constraints active or violated
Duplicate/O fit_wave0,WMCF_TempAutoXWave
convfunc(wave2,fit_wave0,WMCF_TempAutoXWave)
KillWaves/Z WMCF_TempAutoXWave
wave2={0,1.2081,0.99999,90.964,70.751,358.18,28498}
V_chisq= 32.6;V_npnts= 154;V_numNaNs= 0;V_numINFs= 0;
V_startRow= 0;V_endRow= 153;
W_sigma={0,3.42e+04,7.29e+03,1.26e+06,3.24e+06,5.94e+06,7.36e+08}
Coefficient values ± one standard deviation
K0 =0 ± 0
K1 =1.2081 ± 3.42e+04
K2 =0.99999 ± 7.29e+03
K3 =90.964 ± 1.26e+06
K4 =70.751 ± 3.24e+06
K5 =358.18 ± 5.94e+06
K6 =28498 ± 7.36e+08


*wave2 is the coefficient wave*

May 14, 2011 at 02:27 pm - Permalink

Quote:
new error: "while executing Make, the following error occurred: the number of pints in a wave must be between 0 and 2147 million"


The most likely explanation (in fact the only explanation that I can think of) is that your code is asking Igor to make a very big wave:
Variable nExpWavePnts = round((resolutionFactor*pw[6])/dT)*2 + 1
Make/D/O/N=(nExpWavePnts) expwave


Add a print statement before Make:
Print nExpWavePnts, resolutionFactor, pw[6], dT


Or better yet, use the debugger:
DisplayHelpTopic "The Debugger"

May 14, 2011 at 03:37 pm - Permalink

Thanks, the debugger helped me to fix the size problem. Now it will work without any errors, but the fit has huge uncertainty and when graphing it just looks like a straight line...?

*EDIT*

Now i can get it to fit decently with proper constraints and initial guesses. There are 2 problems now:

1) The fitting is very strange as it approaches the last data points. I'm guessing this is due to the convolution.

2) The fitted curve is offset for some reason.

Any suggestions as to fixing these. Thank you!
Graph0.pxp (17.32 KB)

May 15, 2011 at 01:02 pm - Permalink

Aaron_ser wrote:
Thanks, the debugger helped me to fix the size problem. Now it will work without any errors, but the fit has huge uncertainty and when graphing it just looks like a straight line...?

*EDIT*

Now i can get it to fit decently with proper constraints and initial guesses. There are 2 problems now:

1) The fitting is very strange as it approaches the last data points. I'm guessing this is due to the convolution.

2) The fitted curve is offset for some reason.

Any suggestions as to fixing these. Thank you!


1) You are probably seeing convolution end effects. To remove this problem you can a) Calculate the fit function over a wider range than the data requires and trim before you leave the fit function. b) Do a linear extrapolation. c) (if you are using gencurvefit) specify your own cost (chi2) function and tell it to ignore the last points. I rank them in the order a>c>b in order of correctness.

2) Although I've not opened the experiment I'm guessing it's offset because the convolution assumes unit scaling. You have to multiply (or divide) by the spacing between successive x points. (Can't remember which, even though I do it all the time). If you check the offset factor you'll soon find out.

May 15, 2011 at 04:53 pm - Permalink

Aaron_ser wrote:
1) The fitting is very strange as it approaches the last data points. I'm guessing this is due to the convolution.

Quote:
expwave = pw[0]+(x>=0)*(1-exp(-xw/pw[4]))*(pw[1]*exp(-xw/pw[5])+pw[2]*exp(-xw/pw[6]))

If pw[0] is non-zero, then your expwave is non-zero all the way to the ends. That's bad. As the non-zero values fall off the ends of the convolution, the summed products decrease. If the intent is to have a fitted vertical offset, add it at the end of the function rather than in the exponential component:
expwave = (x>=0)*(1-exp(-xw/pw[4]))*(pw[1]*exp(-xw/pw[5])+pw[2]*exp(-xw/pw[6])) // don't add vertical offset here

ywave = exp((-xw/pw[3])^2)
Variable sumywave
sumywave = sum(ywave, -inf,inf)
ywave /= sumywave
convolve/A expwave, ywave
yw = yWave(xw[p]) + pw[0]       // add it here instead


Quote:
2) The fitted curve is offset for some reason.

It's really tricky to get that right in a fit function like this one. When I was developing the example in the documentation, I graphed all the intermediate waves so that I could see what was going on during the computations. If you can't figure it out, send me a copy of your experiment file with your latest code and the data set you show in the graph.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

May 17, 2011 at 02:25 pm - Permalink

Thanks for all the help. I got it working, with a re-dimension to eliminate the end effects and offsetting the gaussian function to fix the offset. i have to say that igor is a lot cooler now that i realize how much you can do with it

May 22, 2011 at 01:52 pm - Permalink

I have one other point to mention. Because of the way that IGOR does the convolution, you _may_ need to horizontally offset(rotate) the convolved data, depending on whether:
a) there is an odd/even number of points (I think an even number of points produces a 1/2 point offset).
b) the dest/src wave are located at x = 0. If they aren't you may need to use the rotate operation.

You can test these effects out by doing convolutions on gaussians for each of the cases. The net effect of them is to shift the convolved data lower/higher in x. If you are using a horizontal offset and you don't care about its value, then you don't need to worry. If the horizontal offset means something you should check this out.

THe offset I was talking about before was the vertical scaling. Try the following code. The result of this will be that dest will be centred at -0.5, not 0 (as you might expect it to be). Change the number of points to 201 and there won't be a problem.

//this codes produces a half point offset.
make/n=200/o/d src,dest
display src, dest
SetScale/I x -100,100,"", src,dest
src=gauss(x, 0, 10)
dest=gauss(x, 0, 10)
convolve/a src, dest
CurveFit/M=2/W=0/TBOX=(0x300) gauss, dest/D
//to correct that offset
dest[] = dest[p + 0.5]

//now lets investigate the effect of offsetting src + dest
//you'll see that the convolved wave is offset to +50.  Thus you have to rotate back to the left.
•SetScale/I x -150,50,"", src,dest
•src=gauss(x, 0, 10)
dest=gauss(x, 0, 10)
convolve/a src,dest
//to correct this offset you have to use rotate

May 22, 2011 at 05:13 pm - Permalink

Aaron_ser wrote:
Thanks for all the help. I got it working, with a re-dimension to eliminate the end effects and offsetting the gaussian function to fix the offset. i have to say that igor is a lot cooler now that i realize how much you can do with it


Hello Aaron,

I am a new igor user and try to use the program to fit the pump-probe data. I think it is very similar with your problem. I am wondering could you please upload your new program. I am trying to learn something form you.

Many thanks,
Best regards
Henry

January 24, 2018 at 11:54 am - Permalink

Aaron_ser wrote:
Thanks for all the help. I got it working, with a re-dimension to eliminate the end effects and offsetting the gaussian function to fix the offset. i have to say that igor is a lot cooler now that i realize how much you can do with it


FWIW, a trick I used to work with FFT smoothing is to baseline correct and mirror the spectrum first. This creates a truer representation of one full cycle of a wave. Here is a brief schematic of the idea.



This helps to remove the Gibb's oscillations that appear in the IFFT operation. Also, the algebra to remove the baseline and mirror+flip the spectrum is entirely reversible on the IFFT result.

ps -- I may have inadvertently meant to post this as a response to questions about the use of FFT to do convolution but will leave it here in any case.

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH
fftprep.png (45.79 KB)

January 25, 2018 at 07:38 am - Permalink