## 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: 297
Joined: 2015-01-20
Location: Germany

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: 297
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
Posts: 1236
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

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
Posts: 1236
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
Posts: 1236
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
Posts: 1236
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
Posts: 1236
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: 50
Joined: 2007-08-28