Normalized Derivative Peak Functions

The built-in fit functions "gauss" and "lor" have parameters optimized for speed rather than readability or physical meaning. Also, they cannot be used to fit "derivative" line shapes occurring, e.g., in electron spin resonance or any other lock-in based spectroscopy. The "Multipeak Fitting" package allows to create fit reports with physical parameters such as the full width at half maximum (FWHM), but it involves a lot of clicking.

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

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More