Covolution function fitting problem

Share icon

Wed, 11/09/2016 - 09:55 pm

Hello,

I need to fit experimentally measured data to the convolution of an resolution function (gaussian) with an square function (for anglular resolution test).
I made a procedure with help of example. but it doesn't work well.
so please help me to edit procedure.

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
Function Convolres(pw, yw, xw) : FitFunc    // linear DOS * Fermi, convolved with Gaussian energy resolution
    Wave pw, yw, xw
 
    // pw[0] = offset
    // pw[1] = energy resolution FWHM
    // pw[2] = square function initial position
    // pw[3] = square function finish position
    // pw[4] = square function amplitude
   
 
   
    // Make the resolution function wave W_res.
    Variable x0 = xw[0]
    Variable dx = (xw[inf]-xw[0])/(numpnts(xw)-1)                           // assumes even data spacing, which is necessary for the convolution anyway
    Make/FREE/D/N=(min(max(abs(3*pw[1]/dx), 5), numpnts(yw))) W_res // make the Gaussian resolution wave
    Redimension/N=(numpnts(W_res)+!mod(numpnts(W_res), 2)) W_res    // force W_res to have odd length
    SetScale/P x, -dx*(numpnts(W_res)-1)/2, dx, W_res
    W_res = gauss( x, 0, pw[1]/(2*sqrt(2*ln(2))) ) +pw[0]
    Variable a = sum(W_res)
    W_res /= a
 
    // To eliminate edge effects due to the convolution, add points to yw, make the spectrum,
    // do the convolution, and then delete the extra points.
    Redimension/N=(numpnts(yw)+2*numpnts(W_res)) xw, yw
    xw = (x0-numpnts(W_res)*dx) + p*dx

    // square function
    variable i
    yw=0+pw[0]
    for(i=pw[2];i<pw[3];i+=1)
    yw[i]=pw[4]
    endfor
   
    Convolve/A W_res, yw
    DeletePoints 0, numpnts(W_res), xw, yw
    DeletePoints numpnts(yw)-numpnts(W_res), numpnts(W_res), xw, yw
End


Thanks!

Irene
Welcome to the forum. There are a few things missing before people can help you:
- You should state more clearly what the function should do and what exactly is failing ("doesn't work well" is rather ambiguous).
- Please put your function in <igor> tags and also don't forget to post the full code including any helping functions (your current post seems cut off).
- Even better would be if you could provide a minimal experiment file together with an explanation what you want to do and at which step stuff fails.
- It also never hurts to state your Igor-Version.

Hope that brings us a step closer to the solution.

November 10, 2016 at 12:23 am - Permalink

First I have a question. Does yw have a fixed scale? If so, I cannot see why you need xw. If not, I cannot see how the convolution of W_res (*) yw will end with the proper result. When the scaling of yw is not identically constant, W_res will be inconsistently scaled each time the position of the Gaussian is moved.

I will assume that yw has a constant scale for what follows.

Now I have some suggestions.

* I think you can use this or an equivalent to create the W_res wave

variable npts = numpnts(yw)
npts = npts - mod(npts,2) // forces value to odd either same as yw or one less than yw
Make/FREE/N=(npts) W_res
SetScale/P x, -round(npts/2), 1, W_res
W_res = gauss(x, 0, pw[1])


I am rather certain the gauss(...) function in Igor returns a normalized Gauss. So the steps to integrate and divide (normalize the curve) are unnecessary.

* I think you can use this to create the square function ...

yw = 0
yw[pw[2],pw[3]] = pw[4]


* Your function has no return value. So, it will do nothing.

* Once you have the changes made, test your function. Dump the results in to a test wave ...

duplicate/O yw testw
testw = Convolve(pw,yw)



--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

November 13, 2016 at 08:23 am - Permalink

It has to have xw in order to match the required function format for an all-at-once fitting function.
An all-at-once fitting function doesn't use the return value; instead it alters the yw input wave.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

November 17, 2016 at 05:10 pm - Permalink

johnweeks wrote:
It has to have xw in order to match the required function format for an all-at-once fitting function.
An all-at-once fitting function doesn't use the return value; instead it alters the yw input wave.


Oh. Yes of course. Thanks!

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH

November 18, 2016 at 06:54 am - Permalink