Multipeak Fit: Fit your experimental data with ... your experimental data?!
Igor's Multipeak Fit is convenient to fit analytical peak functions to data. But there might be times when you want to include some other data (wave) into your fit. For example, to find a match between different measurements or to include reference data into the fit. This code snippet provides the framework to fit experimental data including Gaussian broadening with Multipeak Fit (MPF).
The data is prepared in numbered folders (such as 'root:myData_1:') with numbered wave names (such as 'myPeak0', 'myPeak1' ...). Since MPF can only interact via numerical values, these numbers are later used to select the folder and desired data. After loading MPF, first navigate to Analysis => Multipeak Fit => Prepare Data Peak Function. A small prompt appears where the base folder and wave names are set ("root:myData_" and "myPeak" in above example).
Then start a fit set. Now, the DataPeak fit function is available from within MPF. Here, the last two parameters control the folder and wave number of the data you have just set up. These two parameters need to be on Hold at all times, or the fit will fail. Dial in the desired value to select the data to fit (e.g., folder 1 and wave 2), then run the fit. The first three parameters control the relative x-shift, Gaussian broadening and height (scale) of the data. If you want to switch off broadening, just set and hold this parameter at zero.
You can provide data in two variants: as-is or normalized. If you provide data as-is, then the first parameter can be seen as a 'relative shift' and the third parameter is a 'relative scale'. But you might want to normalize the data before using it in the fit: This would mean to normalize the data height to a maximum of 1 and/or adjust the x-scaling to be centered at zero (so that the 'peak' in your data is placed at x=0). For normalized data, the first parameter can be seen as a 'peak position' and the third parameter is 'height'.
Some notes and caveats:
- You might want to at least remove the baseline from your input data or this baseline will appear in (and be scaled by) the fit.
- If the data's x-range is smaller than needed for the fit, endpoints will be repeated to extend the range.
- The output Location, Amplitude, Area, and FWHM values are numerically determined and may be junk depending on how much your data resembles 'a peak'.
#pragma rtGlobals = 3
#pragma IgorVersion = 8.00
#pragma moduleName = MPFDataFit
#pragma DefaultTab = {3,20,4}
static StrConstant kMPF_WorkingDir = "root:Packages:MultiPeakFit2:"
Menu "Analysis"
Submenu MPFDataFit#CreateDataFitMenu() // add menu entry after loading MPF
"Prepare Data Peak Function", /Q, MPF_PrepareDataFit()
End
End
static Function/S CreateDataFitMenu()
return SelectString(DataFolderExists(kMPF_WorkingDir),"-","Multipeak Fit")
End
// +++++++++
Function MPF_PrepareDataFit() // setup for data locations
String setBaseFolder = "", setBaseName = ""
SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder") // load previous settings
SVAR/Z baseName = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
if (SVAR_Exists(baseFolder))
setBaseFolder = baseFolder
endif
if (SVAR_Exists(baseName))
setBaseName = baseName
endif
Prompt setBaseFolder,"Set base folder name without trailing number (e.g., \"root:myFit_\"):"
Prompt setBaseName,"Set base data name without trailing number (e.g., \"peak_\"):"
DoPrompt "Setup MPF data fit", setBaseFolder, setBaseName
if (V_Flag)
return -1
endif
if (!strlen(setBaseFolder) || !strlen(setBaseName))
Abort "You need to enter both a folder path and base data name."
endif
int i, j, sourceFound = 0
if (!DataFolderExists(kMPF_WorkingDir)) // if MPF was not started yet => create basic folder structure
DFREF saveDF = GetDataFolderDFR()
SetDataFolder root:
for (i = 1; i < ItemsInList(kMPF_WorkingDir,":"); i++)
NewDataFolder/O/S $StringFromList(i,kMPF_WorkingDir,":")
endfor
Variable/G currentSetNumber = 0
SetDataFolder saveDF
endif
setBaseFolder = RemoveEnding(setBaseFolder , ":") // make sure no colon was entered
String/G $(kMPF_WorkingDir+"MPF_DataFit_baseFolder") = setBaseFolder
String/G $(kMPF_WorkingDir+"MPF_DataFit_baseName") = setBaseName
for (i = 0; i < 10 && !sourceFound; i++) // try to look for the data in the first 10 folder / waves
for (j = 0; j < 10 && !sourceFound; j++)
Wave/Z source = $(setBaseFolder+num2str(i)+":"+setBaseName+num2str(j))
sourceFound = WaveExists(source)
endfor
endfor
if (!sourceFound) // still nothing found? => have a simple check to remind the user
DoAlert/T="Checking folder and data" 0, "Your first folder and/or data does not seem to exist yet. Make sure to prepare the data in the format: "+setBaseFolder+"[NUM]:"+setBaseName+"[NUM]"
endif
return 0
End
//########
static Function/Wave getSource(Wave cw) // get source data from saved base names
SVAR/Z baseFolder = $(kMPF_WorkingDir+"MPF_DataFit_baseFolder")
SVAR/Z baseName = $(kMPF_WorkingDir+"MPF_DataFit_baseName")
if (SVAR_Exists(baseFolder) && SVAR_Exists(baseName))
Wave/Z source = $(baseFolder+num2str(cw[3])+":"+baseName+num2str(cw[4]))
if (WaveExists(source) && WaveType(source,1) == 1 && WaveDims(source) == 1)
return source
endif
endif
return $""
End
//########
static Function MPF_setHoldStateForPeak(Wave cw, String holdStr) // injects hold state for the current peak
DFREF DFpath = GetWavesDataFolderDFR(cw)
int peakNum
sscanf NameOfWave(cw), "Peak %d Coefs", peakNum
if (!DataFolderRefStatus(DFpath) || numType(peakNum) != 0 || peakNum < 0)
return -1
endif
SVAR gName = DFPath:GraphName
String pName = gName+"#MultiPeak2Panel" // the MPF panel
if (WinType(pName) == 7)
Wave/T HoldStrings = DFPath:HoldStrings
Execute/Z "WMMultiPeakFit#MPF2_RefreshHoldStrings(\""+pName+"\")" // first, grab the latest hold settings from the panel
HoldStrings[peakNum+1] = holdStr // inject hold string
Execute/Z "WMMultiPeakFit#SetHoldCheckboxesFromWave(WMMultiPeakFit#GetSetNumberFromWinName(\""+pName+"\"))" // rewrite to panel
endif
return 0
End
//####### MPF peak format functions #####
Function/S DataPeak_PeakFuncInfo(Variable which) // pass names to MPF
Switch (which)
case 0:
return "LocShift;Broaden;Height;Folder;DataNo;"
case 1:
return "MPF_DataPeakFunc"
case 2:
return "MPF_GaussToDataPeakGuess"
case 3:
return "MPF_DataPeakParams"
case 4:
return "Location;Height;Area;FWHM;"
EndSwitch
End
//+++++++
Function MPF_GaussToDataPeakGuess(Wave cw) // initializes parameters, and set hold state
// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
Redimension/N=5 cw
int i,j, sourceFound = 0
for (i = 0; i < 10 && !sourceFound; i++)
for (j = 0; j < 10 && !sourceFound; j++)
cw[3] = i
cw[4] = j
Wave/Z source = getSource(cw) // grabs the first valid source wave
sourceFound = WaveExists(source)
endfor
endfor
int normInt = 1, normPos = 1
if (sourceFound)
Variable wMax = WaveMax(source)
Variable left = leftx(source) < rightx(source) ? leftx(source) : rightx(source)
Variable right = leftx(source) < rightx(source) ? rightx(source) : leftx(source)
normInt = wMax > 0.9 && wMax < 1.1 ? 0 : 1
normPos = left < 0 && right > 0 ? 0 : 1
endif
cw[0] = normPos ? 0 : cw[0] // if data seems centered around zero, peak position can stay
cw[1] = 0.5 // initial Gaussian broadening
cw[2] = normInt ? 1 : cw[2] // if data is normalized, peak height can stay
// this needs to be executed later or the hold state will be reset
Execute/P/Q "MPFDataFit#MPF_setHoldStateForPeak("+GetWavesDataFolder(cw,2)+", \"00011\")"
return 0
End
//+++++++
Function MPF_DataPeakFunc(Wave cw, Wave yw, Wave xw) : FitFunc // the fit function
Wave/Z source = getSource(cw)
if (!WaveExists(source))
yw = NaN
if (Exists("MPF2_DisplayStatusMessage")) // display an error in the MPF graph
Execute/Z "MPF2_DisplayStatusMessage(\"Base folder setup wrong or data not found.\","+GetWavesDataFolder(cw,1)+")"
endif
return -1
endif
Duplicate/free source, yWave
SetScale/P x, DimOffset(source,0)+cw[0], DimDelta(source,0), ywave // prepare shifted convolution wave
PeakConvolveGaussHelper(yWave, cw[1]) // Gauss convolution
Variable minX = pnt2x(yWave,0), maxX = pnt2x(yWave,numpnts(yWave)-1) // make sure to not run out-of-bounds (end-point data is repeated)
yw = xw[p] >= minX ? ( xw[p] <= maxX ? cw[2]*yWave(xw[p]) : cw[2]*yWave(maxX) ) : cw[2]*yWave(minX)
return 0
End
//+++++++
Function MPF_DataPeakParams(Wave cw, Wave sw, Wave outWave)
Variable Amplitude = NaN, Location = NaN, FWHM = NaN, PeakArea = NaN // set to NaN for now
FindPeakStatsHelper(cw, Amplitude, Location, FWHM, PeakArea) // approximately calculates all results
// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
// output values are: "Location;Height;Area;FWHM;"
outWave[0][0] = Location // position
outWave[0][1] = NaN
outWave[1][0] = Amplitude // height
outWave[1][1] = NaN
outWave[2][0] = PeakArea // area
outWave[2][1] = NaN
outWave[3][0] = FWHM // FWHM
outWave[3][1] = NaN
return 0
End
//#######
Static Function PeakConvolveGaussHelper(Wave in, Variable width) // convolve with a prepared Gaussian
if (width == 0 || numtype(width) != 0)
return 0
endif
Variable pnt, dx = DimDelta(in,0), size = DimSize(in,0)
width = abs(width/dx)/sqrt(2) // width normalized by x scaling and convert to expected Gauss 'width'
pnt = round(max(abs(10*width),11)) // create a wave 10 times the Gauss width (5 sigma on each side)
pnt = min(pnt, 2*numpnts(in)+1) // doesn't have to be larger than two times the full input data
pnt = pnt+!mod(pnt, 2) // force odd size
Make/FREE/D/N=(pnt) GaussWave
GaussWave = Gauss(x, (pnt-1)/2, width)
Variable A = sum(GaussWave) // make sure the Gauss area limited even for a few points
if (A > 1)
GaussWave /= A
endif
Make/FREE/D/N=(2*pnt+size) temp = in[0] // enlarge the wave to get the full convolution range
temp[pnt,size+pnt-1] = in[p-pnt]
temp[size+pnt,size+2*pnt-1] = in[size-1]
MatrixOP/Free temp = replaceNaNs(temp,0) // makes sure there are no NaN values in-between
Convolve/A GaussWave temp
in[0,size-1] = temp[p+pnt]
End
//#########
// recalculate the (almost) full peak and extract all result values numerically
Static Function FindPeakStatsHelper(Wave cw, Variable &Amplitude, Variable &Location, Variable &FWHM, Variable &PeakArea)
Wave/Z sw = getSource(cw)
if (!WaveExists(sw))
return -1
endif
// cw values are: "PosShift;GaussWidth;Height;Folder;DataNo;"
Variable left = pnt2x(sw,0)+cw[0], right = pnt2x(sw,numpnts(sw)-1)+cw[0]
Variable pntSize = DimSize(sw,0)*10 // create a big wave (fine increments)
pntSize = pntSize > 100000 ? 100000 : pntSize // make sure its not too big
Make/FREE/D/N=(pntSize+1) yw, xw
Make/FREE/D/N=2 Levels
SetScale/I x, left, right, xw; xw = x
//+++++ find values +++++
MPF_DataPeakFunc(cw, yw, xw) // recalculate the peak with the coef set
PeakArea = AreaXY(xw, yw, left, right) // the area is calculated for the full source range
WaveStats/Q yw
Amplitude = cw[2] > 0 ? V_max : V_min
Location = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
FindLevels/Q/D=Levels yw, Amplitude/2 // find the half-maximum positions
if (V_LevelsFound == 2)
left = xw[Levels[0]]
right = xw[Levels[1]]
FWHM = abs(left - right)
endif
if (numtype(FWHM) != 0)
return 0
endif
//++++++ recalculate with higher resolution ++++++
// choose a range just a bit bigger than the FWHM
SetScale/I x, left - 2*FWHM/(pntsize+1), right + 2*FWHM/(pntsize+1), xw; xw = x
MPF_DataPeakFunc(cw, yw, xw)
WaveStats/Q yw
Amplitude = cw[2] > 0 ? V_max : V_min
Location = cw[2] > 0 ? xw[V_maxloc] : xw[V_minloc]
FindLevels/Q/D=Levels yw, Amplitude/2
if (V_LevelsFound == 2)
left = xw[Levels[0]]
right = xw[Levels[1]]
FWHM = abs(left - right)
endif
return 0
End
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More