Iterating a curve fit in a varying range

remolorio
Posts: 10
Joined: 2016-09-15
Location: United Kingdom

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

AttachmentSize
Auto_corr.pxp4.96 KB

HJDrescher
Posts: 332
Joined: 2015-01-20
Location: Germany

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


remolorio
Posts: 10
Joined: 2016-09-15
Location: United Kingdom

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.


HJDrescher
Posts: 332
Joined: 2015-01-20
Location: Germany

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

AttachmentSize
AC.png26.85 KB

jjweimer
jjweimer's picture
Posts: 1288
Joined: 2007-08-14
Location: United States

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


[ last edited June 8, 2017 - 04:50 ]
remolorio
Posts: 10
Joined: 2016-09-15
Location: United Kingdom

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
<igor>
 
EDIT: here is a tidier version of my macro

AttachmentSize
Spectrum.pxp23.12 KB
Salje-Autocorrelation analysis of infrared spe.pdf601.76 KB

[ last edited June 9, 2017 - 03:28 ]
jjweimer
jjweimer's picture
Posts: 1288
Joined: 2007-08-14
Location: United States

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


jjweimer
jjweimer's picture
Posts: 1288
Joined: 2007-08-14
Location: United States

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

AttachmentSize
Auto Correlation Analysis.pxp14.13 KB

remolorio
Posts: 10
Joined: 2016-09-15
Location: United Kingdom

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
<igor>

AttachmentSize
Ugly_steps.png195.86 KB

[ last edited June 11, 2017 - 08:19 ]
jjweimer
jjweimer's picture
Posts: 1288
Joined: 2007-08-14
Location: United States

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


jjweimer
jjweimer's picture
Posts: 1288
Joined: 2007-08-14
Location: United States

Here is a new version. The panel works on traces on the frontmost graph. It uses FuncFit to a Gauss function centered at zero.

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

AttachmentSize
Auto Correlation Analysis v1.1.pxp16.06 KB
autocorranalysis.png287.04 KB

tony
Posts: 104
Joined: 2007-08-28
Location: Canada

I have played a bit with this sort of analysis.
Here's my procedure file. It hasn't been tested much, but worked when I used it for some Raman and IR spectra.

AttachmentSize
AutocorrelationAnalysis.ipf5.65 KB

Back to top