Normalized Derivative Peak Functions
harneit
As a "pedestrian" alternative for those preferring command-line operations, here's a collection of properly normalized peak fitting functions often used in spectroscopy: Gaussian, Lorentzian, and pseudo-Voigt shapes; both for normal and derivative peaks.
There are four versions for each shape: two versions have the visible amplitude of the peak as the intensity parameter, which makes providing good guesses easier. The other two have the integrated area as the intensity parameter (useful for estimating scattering strength, analyte concentrations, etc.). Each of these variants are available either for absorption or for derivative line shapes. For the absorption shapes, the width parameter is the FWHM while for the derivative peaks, it is the peak-to-peak width.
For more detail, see the documentation in the code. It also explains how to use these functions with
FuncFit's cool feature of "fitting sums of functions". The functions provided here have no constant term, so it is not necessary to hold those terms in order to avoid the "linear-dependency" problem. Rather, one could add another fitting function (e.g., lin) to simultaneously fit the background.I checked the formulae used in these functions thoroughly but if you detect errors, please let me know.
#pragma rtGlobals=1 // Use modern global access method.
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Normalized peak functions (Gauss, Lorentz, pseudo-Voigt) in absorption and dispersion mode.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// by W. Harneit (Berlin), November 2009
// August 2010: Updated comments, included note on multi-peak fitting
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Each function takes three parameters (in wave w): intensity, position, and width.
//
// For the ...Amp functions, the intensity parameter means amplitude (from zero).
// The area is not normalized in this case.
// For the ...Area functions, intensity means integrated area (or double integral for Diff...).
// The intensity parameter does not correspond to the peak amplitude in this case.
//
// For the absorption functions (without Diff...), width means full width at half maximum (FWHM).
// The Diff... functions are the derivatives (dispersion mode). Width is peak-to-peak width (p.p. width).
// p.p. width (dispersion) and FWHM (absorption) are linearly related, depending on the peak form.
//
// The Pseudo-Voigt functions are a simple linear interpolation between Gauss and Lorentz form, with
// _equal_ width. They take a fourth parameter 0 <= w[3] <= 1, which is the weight of the Lorentzian
// (the weight of the Gaussian is 1 - w[3]). Strange things will happen if w[3] > 1.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Examples:
//
// GaussAmp is the absorptive Gaussian, with amplitude w[0], position w[1], fwhm w[2].
//
// DiffLorArea is the dispersive Lorentzian, with double integral w[0], position w[1], p.p. width w[2].
//
// To do multi-peak fitting, follow the steps below.
// (see also DisplayHelpTopic "Fitting Sums of Fit Functions")
//
// We assume that a two-line absorption spectrum is present in spec_INT. The position parameters
// are the most sensitive parameters, followed by the width. Intensities (here: areas) may be way off.
// Tip: You can use "Decommentize" (Edit menu) on the following lines.
//
// // Step 1: provide guesses for the two peaks
// Make/D/O low={4e-10,0.33707,2e-5}, high={4e-10,0.33818,2e-5, 0.77}
// // Step 2: provide hold strings (if any)
// String hstr = "0001" // to hold the shape parameter p=0.77 for the "high" peak
// // Step 3: do the fit
// FuncFit/M=0/L=4096/W=0/X {{LorArea,low}, {PseudoVoigtArea,high,hold=hstr}}, spec_INT /D
// // /M=0 => no cov. matrix; /L= numpnts; /W=0 => don't wait; /X => full X range; /D => auto dest.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// code note: Numerical constants are given in the comments with 5 digits precision.
// If speed is an issue, one could use global constants "constant kln2 = 0.69315" etc.
//
function GaussAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2) // ln2 = 0.69315
return w[0] * exp(-ln2*xi*xi)
end
function GaussArea(w, x)
wave w
variable x
return GaussAmp(w, x) / ( sqrt( pi / ln(2) ) / 2 * w[2] ) // divisor = 1.0645 * w[2]
end
function LorAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2]
return w[0] / (1 + xi*xi)
end
function LorArea(w, x)
wave w
variable x
return LorAmp(w, x) / ( pi / 2 * w[2] ) // divisor = 1.5708 * w[2]
end
function DiffGaussAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2) // ln2 = 0.69315
return -xi * w[0] * exp(-0.5*(xi*xi-1))
end
function DiffGaussArea(w, x)
wave w
variable x
return DiffGaussAmp(w, x) / ( sqrt( pi * e / 8 ) * w[2] * w[2] ) // divisor = 1.0332 * w[2]^2
end
function DiffLorAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2]
return -xi * w[0] * ( 4 / (3 + xi*xi) )^2
end
function DiffLorArea(w, x)
wave w
variable x
return DiffLorAmp(w, x) / ( sqrt( 4 / 3 ) * w[2] * w[2] ) // divisor = 1.1547 * w[2]^2
end
function PseudoVoigtAmp(w, x)
wave w
variable x
return w[3] * LorAmp(w, x) + (1 - w[3]) * GaussAmp(w, x)
end
function PseudoVoigtArea(w, x)
wave w
variable x
return w[3] * LorArea(w, x) + (1 - w[3]) * GaussArea(w, x)
end
function DiffPseudoVoigtAmp(w, x)
wave w
variable x
return w[3] * DiffLorAmp(w, x) + (1 - w[3]) * DiffGaussAmp(w, x)
end
function DiffPseudoVoigtArea(w, x)
wave w
variable x
return w[3] * DiffLorArea(w, x) + (1 - w[3]) * DiffGaussArea(w, x)
end
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Normalized peak functions (Gauss, Lorentz, pseudo-Voigt) in absorption and dispersion mode.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// by W. Harneit (Berlin), November 2009
// August 2010: Updated comments, included note on multi-peak fitting
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Each function takes three parameters (in wave w): intensity, position, and width.
//
// For the ...Amp functions, the intensity parameter means amplitude (from zero).
// The area is not normalized in this case.
// For the ...Area functions, intensity means integrated area (or double integral for Diff...).
// The intensity parameter does not correspond to the peak amplitude in this case.
//
// For the absorption functions (without Diff...), width means full width at half maximum (FWHM).
// The Diff... functions are the derivatives (dispersion mode). Width is peak-to-peak width (p.p. width).
// p.p. width (dispersion) and FWHM (absorption) are linearly related, depending on the peak form.
//
// The Pseudo-Voigt functions are a simple linear interpolation between Gauss and Lorentz form, with
// _equal_ width. They take a fourth parameter 0 <= w[3] <= 1, which is the weight of the Lorentzian
// (the weight of the Gaussian is 1 - w[3]). Strange things will happen if w[3] > 1.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Examples:
//
// GaussAmp is the absorptive Gaussian, with amplitude w[0], position w[1], fwhm w[2].
//
// DiffLorArea is the dispersive Lorentzian, with double integral w[0], position w[1], p.p. width w[2].
//
// To do multi-peak fitting, follow the steps below.
// (see also DisplayHelpTopic "Fitting Sums of Fit Functions")
//
// We assume that a two-line absorption spectrum is present in spec_INT. The position parameters
// are the most sensitive parameters, followed by the width. Intensities (here: areas) may be way off.
// Tip: You can use "Decommentize" (Edit menu) on the following lines.
//
// // Step 1: provide guesses for the two peaks
// Make/D/O low={4e-10,0.33707,2e-5}, high={4e-10,0.33818,2e-5, 0.77}
// // Step 2: provide hold strings (if any)
// String hstr = "0001" // to hold the shape parameter p=0.77 for the "high" peak
// // Step 3: do the fit
// FuncFit/M=0/L=4096/W=0/X {{LorArea,low}, {PseudoVoigtArea,high,hold=hstr}}, spec_INT /D
// // /M=0 => no cov. matrix; /L= numpnts; /W=0 => don't wait; /X => full X range; /D => auto dest.
//
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// code note: Numerical constants are given in the comments with 5 digits precision.
// If speed is an issue, one could use global constants "constant kln2 = 0.69315" etc.
//
function GaussAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2) // ln2 = 0.69315
return w[0] * exp(-ln2*xi*xi)
end
function GaussArea(w, x)
wave w
variable x
return GaussAmp(w, x) / ( sqrt( pi / ln(2) ) / 2 * w[2] ) // divisor = 1.0645 * w[2]
end
function LorAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2]
return w[0] / (1 + xi*xi)
end
function LorArea(w, x)
wave w
variable x
return LorAmp(w, x) / ( pi / 2 * w[2] ) // divisor = 1.5708 * w[2]
end
function DiffGaussAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2], ln2 = ln(2) // ln2 = 0.69315
return -xi * w[0] * exp(-0.5*(xi*xi-1))
end
function DiffGaussArea(w, x)
wave w
variable x
return DiffGaussAmp(w, x) / ( sqrt( pi * e / 8 ) * w[2] * w[2] ) // divisor = 1.0332 * w[2]^2
end
function DiffLorAmp(w, x)
wave w
variable x
variable xi = 2 * (x - w[1]) / w[2]
return -xi * w[0] * ( 4 / (3 + xi*xi) )^2
end
function DiffLorArea(w, x)
wave w
variable x
return DiffLorAmp(w, x) / ( sqrt( 4 / 3 ) * w[2] * w[2] ) // divisor = 1.1547 * w[2]^2
end
function PseudoVoigtAmp(w, x)
wave w
variable x
return w[3] * LorAmp(w, x) + (1 - w[3]) * GaussAmp(w, x)
end
function PseudoVoigtArea(w, x)
wave w
variable x
return w[3] * LorArea(w, x) + (1 - w[3]) * GaussArea(w, x)
end
function DiffPseudoVoigtAmp(w, x)
wave w
variable x
return w[3] * DiffLorAmp(w, x) + (1 - w[3]) * DiffGaussAmp(w, x)
end
function DiffPseudoVoigtArea(w, x)
wave w
variable x
return w[3] * DiffLorArea(w, x) + (1 - w[3]) * DiffGaussArea(w, x)
end
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More