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 TextEncoding = "UTF-8"
#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