Iterating a curve fit in a varying range

Share icon

Wed, 06/07/2017 - 10:58 am

Hi

I have correlated a Raman spectrum with itself (autocorrelation, see attachement). I now would like to fit a gaussian curve to the central peak of the resulting wave. This should be done multiple times for a varying range of X_values (say from -20 to 20, -19 to 19, -18 to 18, etc. to -2 to 2).

The Gaussian function looks like this: gaus=w[0]*exp(-((x-w[1])/w[2])^2), with w being the coefficient wave, and x the variable. I suppose using the FuncFit operation?

The values of w[2] of all fits should be stored in one result wave.

How can this be done? Ideally there is a certain flexibility for the boundaries and increments chosen, depending on the case.

I am a beginner and this is beyond my programming skills.

Thanks
Remo Auto_corr.pxp
What's your goal?
Fitting with an x-range of -2 to 2 will use three data points.
Do you actually want to perform a deconvolution of your AC into three peaks? Have a look at the multi-fit package (in the Analysis menu).

HJ

June 7, 2017 at 12:35 pm - Permalink

I would like to plot the widths of gaussian curves against the ranges for which they were fitted. ultimately the idea is to fit a curve to these values in turn and extrapolate a width for fitting the range → 0,( "the central peak"). So it is not a multi peak fit, one peak at a time. This is called autocorrelation analysis, developed for vibrational spectroscopic data with strongly overlapping peaks.

June 7, 2017 at 02:20 pm - Permalink

You would have to run curvefit from the menu and use the cursors as input range.
CurveFit gauss Auto_corr[pcsr(A),pcsr(B)] /D and replace "pcsr" by your range (in data points)
If you want to automate it, have a look at the manual regarding loops:
displayhelptopic "loops"

However, it seems that the result for a narrow range is quite similar for the result from the deconvoluted peaks. It sounds to me like this method is some kind of approximation (eliminating overlapping sidebands -- I haven't read details, I admit -- just my feeling).
Have a look at the attached image. The black bold line is the center peak from a multi peak fit (middle graph; AC data as '+') and it looks like the original signal in the range of approximately (-7 ...7).
HJ
AC.png

June 7, 2017 at 03:46 pm - Permalink

I might ask, as a follow up to @HJDrescher, what reference source are you using for the algorithm to do what you propose? We could benefit by seeing an example of why and how someone else does what you want to do (presumably using some other software application).

The option proposed by @HJDrescher to peak fit and thereby resolve (and apologies but ... speaking strictly and pedantically, deconvolution is an improper term here) the components to the broad peak may also miss something fundamental that is determined by the approach you want to take. We'd need to see a reference otherwise.

Finally, I would be surprised when the approach that you propose does not lead to a "step-function-like" plot of half-width versus fit range. At the largest fit range, the peak width will be consistently broad. At some range, the width with snap-down to the narrow peak. At some final range, the width will essentially snap-down to the step size between points. Again, a reference to see where what you want to do has been applied (successfully) somewhere else would help.

ps ... Yes, what you want to do can be programmed to be a routine that you can run automatically. Further recommendations may follow later.

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

June 8, 2017 at 05:50 am - Permalink

Thanks for all your comments.

I have attached an article describing this approach. The main recipe is on pp508.

I have invested a few hours and tackled this problem myself. The whole thing is a patch work of snippets and adaptions thereof from stuff I found.

It seems to work in principle, although I do get "out of range" errors at the end. The Aim is to be able to select any subrange of the input wave and perform the analysis. The Value of interest eventually is y of x=0 of the polynomial fit to the various gaussian widths.


#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.
Macro AutoCorrelation(scrwave)
    string scrwave,subrange,prefix,bline,k0_coef,k0_coef_sigma,Dcorr,Dcorr_sigma,autocorr
   
        Prompt scrwave,"Spectrum to analyse",popup Wavelist("*",";","") //select source wave, subrange, and lable in pop-up menu
   
      rangeinput()
     
        subrange=lable+"_range_"+scrwave
        bline="bl_"+scrwave
       
        autocorr=lable+"_autocorr_"+scrwave
       
        k0_coef="k0_"+lable+"_"+scrwave
        k0_coef_sigma="k0_sigma_"+lable+"_"+scrwave
       
        Dcorr="Dcorr_"+scrwave
        Dcorr_sigma="DcorrSigma_"+scrwave
       
            make/N=0   $k0_coef
        make/N=0   $k0_coef_sigma
        make/N=0   $Dcorr
        make/N=0   $Dcorr_sigma
       
        duplicate/O/R=(gX1,gX2) $scrwave,$subrange    //Duplicate a subrange of the source wave, specified by x-value
        duplicate/O $subrange,$bline
   
        Variable x1,x2,y1,y2
        x1=leftx($subrange)
        x2=rightx($subrange)
        y1=$subrange[0]
        y2=$subrange[numpnts($subrange)-1]
        $bline= ((y1-y2)/(x1-x2))*x+x1*y2
        $subrange=$subrange-$bline
       
           duplicate/O $subrange,$autocorr
           Correlate/AUTO $subrange, $autocorr             
   
        GaussIteration(max_deltaw,increment,$autocorr,$k0_coef,$k0_coef_sigma)
   
       CurveFit poly 3, $k0_coef[start_xtrp,max_deltaw] /F={0.95, 4}  //fit the widths of the Gaussian curve fits ("$k0_coef") by poly3
       
       Variable Dc
        Dc = W_coef[0]
    Redimension /N=(numpnts($Dcorr)+1) $Dcorr
       $Dcorr [numpnts($Dcorr)] = Dc
   
        Variable Dc_sigma
        Dc_sigma = W_ParamConfidenceInterval[0]
      Redimension /N=(numpnts($Dcorr_sigma)+1) $Dcorr_sigma
        $Dcorr_sigma [numpnts($Dcorr_sigma)] = Dc_sigma
End

Function GaussIteration(plusminusX,incrementX,ywave,output,output_sigma)
 
   variable plusminusX,incrementX
   wave ywave, output, output_sigma
    variable i

    for(i= 3; i<=plusminusX; i+=incrementX)
        CurveFit/Q gauss ywave((-i),i) /D  
        Wave W_coef, W_sigma

        Variable k0
        k0 = W_coef[3]
        Redimension /N=(numpnts(output)+1) output
        output [numpnts(output)-1] = k0
        SetScale/P x 3,incrementX,"", output
       
       Variable k0_sigma
       k0_sigma = W_Sigma[3]
       Redimension /N=(numpnts(output_sigma)+1) output_sigma
       output_sigma [numpnts(output_sigma)-1] = k0_sigma
       SetScale/P x 3,incrementX,"", output_sigma
    endfor
End
<pre><code class="language-igor">

EDIT: here is a tidier version of my macro
Salje-Autocorrelation analysis of infrared spe.pdf Spectrum.pxp

June 9, 2017 at 04:28 am - Permalink

That is certainly an interesting approach to analysis. Thank you for the reference. I have to wonder whether it can be translated to other spectroscopies and how it might supplement such statistical analysis methods as PCA/PCR.

I cannot comment directly now on improvements to what you have but I do want to follow up once my schedule permits.

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

June 8, 2017 at 11:37 am - Permalink

Here is an example experiment that only calculates and stores the correlation results (width and standard uncertainty) when it has what you term as the auto_corr wave.

Eventually, I'd be interested to expand this to a method to do the entire analysis starting from a spectrum input. That's going to take a bit more though.

--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAH
Auto Correlation Analysis.pxp

June 9, 2017 at 12:46 pm - Permalink

Great job

when it does the gaussian fit to a subrange of the auto_corr wave, it fits a value for y0. In the orignal procedure proposed in literature however the gaussian fit is done always with y0 = 0.
I have tried to implement this behaviour in my own attempt by defining K0=0 and add the flag /H="1000", which corresponds to y0 for a gaussian curve, this seems to work, however there is an ugly step in the resulting fit curves. How would you do this?

here is my humble try
Function GaussIteration(plusminusX,incrementX,ywave, output,output_sigma,graphname)
 
   variable plusminusX,incrementX
   wave ywave, output, output_sigma
   string graphname
    variable i
    variable y0=0
       
    K0 = y0
   
    Display /N=graphname
   
    for(i= 15; i<=plusminusX; i+=incrementX)
       String Gauss_fit=graphname+"_"+num2str(i)
       Duplicate/O  ywave, $Gauss_fit
       
   
        CurveFit/Q/H="1000" gauss ywave((-i),i)/D=$Gauss_fit

        Wave W_coef, W_sigma
        Variable width
        width = W_coef[3]
        Redimension /N=(numpnts(output)+1) output
        output [numpnts(output)-1] = width
        SetScale/P x 15,incrementX,"", output
       
       Variable width_sigma
       width_sigma = W_Sigma[3]
       Redimension /N=(numpnts(output_sigma)+1) output_sigma
       output_sigma [numpnts(output_sigma)-1] = width_sigma
       SetScale/P x 15,incrementX,"", output_sigma
       
       appendtograph /W=graphname $Gauss_fit,ywave
    endfor
End
<pre><code class="language-igor">
Ugly_steps.png

June 11, 2017 at 09:19 am - Permalink

My next recommendations are as follows:

* Write the exponential function as a self-made fitting function where y0 is absent. Fit to the home-made function instead of the in-built gauss function. This avoids the need to use a HOLD wave or coefficient term.

* Write the procedure so that it can be added by an #include statement and works on a wave in the front-most graph. This will make the package portable. You will be able to generate auto_corr waves from spectra, store them in specific folders, #include the correlation analysis package, run it on the different auto_corr waves in different graphs, and obtain the output values of correlation value versus step size.

On the front end ... loading and correlation analysis of the spectra are left to you or to a much later project. On the back end, to save analysis results as independent experiment sets, I recommend the SnapIt Package found here http://www.igorexchange.com/project/SnapIt.

So, for the next step, here's what you might do. Write a GaussFixedPosition function to replace the built-in Gauss function. Modify the line in the code I provided to use FitFunc with the modified GaussFixedPosition function instead of CurveFit with the built-in Gauss function. I'll wrap the functions provided in a UI panel that works on a wave in the frontmost graph (much like the SnapIt package works on waves on the frontmost graph).

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

June 12, 2017 at 09:50 am - Permalink