User Defined Peak Shape in Multipeak fitting 2.0

Share icon

Hi everyone,

I was hoping to add a user defined peak shape to my multipeakfitting 2.0 toolbox, but I am struggling to understand how to do thins. I have looked into the Peakfunctions2.ipf macro that has the peak functions, but I am unsure how to go about adding my own function. If someone could walk me through doing this, or may already have this peak shape written into a macro that they are willing to share I would appreciate it.

The function I want to add is a Breit–Wigner Fano peak shape. This is an asymmetric peak that is pretty common when looking at Raman spectroscopy of amorphous carbons. The function is as follows, sorry for the horrible formatting of the equation here.

I(ω) = (A/(1+(2(ω+ω)/T)^2)) + B[1+2((ω- ω2)/qT2)]2/[1+2((ω- ω2)/T2)]2.

Any help or guidance on how to properly get this added to my fitting options would be appreciated.

Thanks
Yeah, Multipeak Fit 2 makes it possible to add your own peak shapes, not easy to do so.

I started looking at the problem and I think your equation has a typo in it. Should this: 2(ω+ω) be something more like 2(ω-ω2)? I guess ω2 is what is commonly rendered as x0.

And can you give me some info on what controls width and amplitude?

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

June 10, 2015 at 11:26 am - Permalink
Yes, it is not extremely easy, but if you understand how it works it is rather straight-forward. You basically need to write four functions:
1) a function for naming you input/output parameters, as well as telling MPF the names of the other three functions.
2) a function calculating the shape of you peak using an x & y-wave input and a parameter set
3) a function calculating useful initial guesses from the standard asymmetric gaussian (the shape you draw when adding a peak) handed from MPF, or at least plug in the ones like width, height at the right places for your peak to work
4) a function calculating useful output parameters from you fit parameters (i.e., FWHM from width etc.). If you feel lazy just hand over the fit parameters and do the rest later.
Peeking at the at the already implemented one's (PeakFunctions2.ipf) might help a bit to get a grasp how stuff works. Also, you might want to read here a bit: http://www.igorexchange.com/node/5750
I attach my definition of the Doniac-Sunjic peak-shape for reference (this one is already rather advanced, as it involves Gauss-convolution to simulate experimental broadening, so please don't feel demotivated by this).


//________________________________________________________________________________________________//
//	MultiPeak-Fit function:	Doniac-Sunjic shape
//________________________________________________________________________________________________//

Function/S DS_PeakFuncInfo(InfoDesired)
	Variable InfoDesired
 
	String info=""
 	Switch (InfoDesired)
		case PeakFuncInfo_ParamNames:
			info = "Location;Width;Height;Asymmetry;GaussW;"
			break;
		case PeakFuncInfo_PeakFName:
			info = "DSPeak"
			break;
		case PeakFuncInfo_GaussConvFName:
			info = "GausstoDSGuess"
			break;
		case PeakFuncInfo_ParameterFunc:
			info = "DSParams"
			break;
		case PeakFuncInfo_DerivedParamNames:
			info = "Location;Height;Area;FWHM;Asymmetry;GaussWidth;"
			break;
		default:
			break;
	endSwitch
 	return info
End
//################################################################################################
Function DSPeak(w,yw,xw)
	Wave w
	Wave yw, xw
 
 	Variable Resolution	= 10								 // the sample resolution (higher = slower)
	Variable Xdelta		= abs(xw[0] - xw[NumPnts(xw)-1]) / (Resolution * NumPnts(xw))		 // calculate an high-res xdelta for scaling
	Variable yPoints	= Resolution * NumPnts(yw)			 // make the wave X times bigger than the original y wave
	Make/FREE/D/N=(2*yPoints+1) yWave						 // create a wave two times as big to have space before and after the peak
	SetScale/P x xw[0]-yPoints/2*Xdelta, Xdelta, yWave			 // reserve 1/2 size of space before the x wave starts (1/2 after)
 
	yWave = w[2]*cos(Pi*w[3]/2+(1-w[3])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[3])/2)
	
	PeakConvolveGauss(yWave, w[4])
	
	yw = yWave(xw[p])		 // write the calculated values into the resultwave
End
//################################################################################################
Function GaussToDSGuess(w)
	Wave w
	
	Variable x0	= w[0]	 // extract the initial guesses
	Variable width	= w[1]
	Variable height	= w[2]
	Variable lwidth	= w[3]
	Variable rwidth	= w[4]
	
	Redimension/N=5 w		 // redimension to apropriate values
	// output values are: "Location;Width;Height;Asym;GaussW;"
	w[1] = 0.8*width	 // 80% into lorentzian width
	w[2] = 0.8*height	 // make the peak a bit smaller
	w[3] = 5*(rwidth-lwidth)/width	 // asymmetry depending on width difference
	w[4] = 0.3*width	 // 30% into gaussian width
	return 0
end
//################################################################################################
Function DSParams(cw, sw, outWave)
	Wave cw, sw, outWave

	Variable Amplitude, Location, FWHM, PeakArea
	Peak_FindStats(cw, "DSPeak", Amplitude, Location, PeakArea, FWHM) // approximately calculates all the results from the parameters

	Redimension/N=(6,2) outWave
	// output values are: "Location;Height;Area;FWHM;Asym;GaussWidth"	
	outWave[0][0] = Location				 // position
	outWave[0][1] = sqrt(sw[0][0])
	outWave[1][0] = Amplitude				 // height
	outWave[1][1] = NaN					 // no idea how to calculate the error
	outWave[2][0] = PeakArea				 // area (probably wrong)
	outWave[2][1] = NaN
	outWave[3][0] = FWHM				 // FWHM
	outWave[3][1] = NaN
	outWave[4][0] = cw[3]				 // asymmetry value
	outWave[4][1] = sqrt(sw[3][3])
	outWave[5][0] = cw[4]				 // gauss broadening width
	outWave[5][1] = sqrt(sw[3][3])
End
//###################################### Value finder helper function ###################################
Function Peak_FindStats(cw, PeakFunc, Amplitude, Location, PeakArea, FWHM) // recalculate the FULL peak and extract all result values (approximately)
	Wave cw
	String PeakFunc
	Variable &Amplitude, &Location, &PeakArea, &FWHM
	
	Funcref PCISimplePeak PeakFunction = $PeakFunc		 // reference the passed peak function
	
	Variable x0	= cw[0]
	Variable width	= cw[1] + cw[4]
	Variable PosH	= cw[2] > 0				 // Look if height is positive
	
	Variable Size = 3000						 // create a big wave (fine increments)
	Make/FREE/N=(Size+1) yw, xw
	xw = (p/Size-0.5) * 30 * width + x0			 // create a broad x scale (to make sure the full peak fits in)
	PeakFunction(cw, yw, xw)					 // recalculate the peak with the coef set
	WaveStats/Q yw							 // find the statistics

	if (PosH)
		Amplitude	= V_max
		Location	= xw(V_maxloc)
	else
		Amplitude	= V_min
		Location	= xw(V_minloc)
	endif
	
	FindLevel/Q/R=(0,Size/2) yw, Amplitude/2
	Variable Left	= V_flag == 1 ? NaN : xw(V_LevelX)
	FindLevel/Q/R=(Size/2,Size) yw, Amplitude/2
	Variable Right	= V_flag == 1 ? NaN : xw(V_LevelX)	
	
	FWHM		= abs(Left -Right)
	PeakArea	= AreaXY(xw, yw)
	return 0
End

June 10, 2015 at 07:03 pm - Permalink