Fitpeaks: Difference between revisions

From Eigenvector Research Documentation Wiki
Jump to navigation Jump to search
imported>Jeremy
(Importing text file)
imported>Jeremy
(Importing text file)
Line 1: Line 1:
===Purpose===
===Purpose===
Peak fitting routine.
Peak fitting routine.
===Synopsis===
===Synopsis===
:[peakdefo,fval,exitflag,out,fit,res] = fitpeaks(peakdef,y,ax,options)
:[peakdefo,fval,exitflag,out,fit,res] = fitpeaks(peakdef,y,ax,options)
===Description===
===Description===
Based on the initial guess in input peakdef, FITPEAKS estimates the peak fit (also the Jacobian and Hessian), and makes a call to LMOPTIMIZEBND to find the best fit of the peaks to the data. (See LMOPTIMIZEBND for additional information.) Results are output to peakdefo.
Based on the initial guess in input peakdef, FITPEAKS estimates the peak fit (also the Jacobian and Hessian), and makes a call to LMOPTIMIZEBND to find the best fit of the peaks to the data. (See LMOPTIMIZEBND for additional information.) Results are output to peakdefo.
Information about individual peaks is stored in standard peak structures (see PEAKSTRUCT). Information on multiple peaks is stored in a multi-record structure. Given a standard peak structure (peakdef) that contains an initial guess of peak locations and widths, FITPEAKS finds new parameters that best fits peaks to the rows of the '' '' data matrix (y). Results are output to a standard peak structure (peakdefo).
Information about individual peaks is stored in standard peak structures (see PEAKSTRUCT). Information on multiple peaks is stored in a multi-record structure. Given a standard peak structure (peakdef) that contains an initial guess of peak locations and widths, FITPEAKS finds new parameters that best fits peaks to the rows of the '' '' data matrix (y). Results are output to a standard peak structure (peakdefo).
Fields of (peakdef) required in the initial guess for each peak are (.fun), (.param), (.lb), (.penlb), (.ub), and (.penub).
Fields of (peakdef) required in the initial guess for each peak are (.fun), (.param), (.lb), (.penlb), (.ub), and (.penub).
====INPUTS====
====INPUTS====
* '''peakdef''' = multi-record standard peak structure with the following fields:
* '''peakdef''' = multi-record standard peak structure with the following fields:
* '''name''': 'Peak', name indicating that this is a standard peak structure.
* '''name''': 'Peak', name indicating that this is a standard peak structure.
* '''id''': ' ', double or character string peak identification.
* '''id''': ' ', double or character string peak identification.
* '''fun''': [ {'Gaussian'} | 'Lorentzian' | 'PVoigt1' | 'PVoigt2' ], defines the peak function (see definitions in the Algorithm section).
* '''fun''': [ {'Gaussian'} | 'Lorentzian' | 'PVoigt1' | 'PVoigt2' ], defines the peak function (see definitions in the Algorithm section).
* '''param''': Parameter list for each peak function. The number of parameters depends on the peak function used:
* '''param''': Parameter list for each peak function. The number of parameters depends on the peak function used:
*  ''''Gaussian'''': [height, location, width],
*  ''''Gaussian'''': [height, location, width],
*  ''''Lorentzian'''': [height, location, width],
*  ''''Lorentzian'''': [height, location, width],
*  ''''PVoigt1'''': [height, location, width, fraction Gaussian],
*  ''''PVoigt1'''': [height, location, width, fraction Gaussian],
*  ''''PVoigt2'''': [height, location, width, fraction Gaussian].
*  ''''PVoigt2'''': [height, location, width, fraction Gaussian].
* '''lb''': [ ], Lower bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
* '''lb''': [ ], Lower bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
* '''penlb''': [ ], Penalty wt for lower bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
* '''penlb''': [ ], Penalty wt for lower bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
* '''ub''': [ ], Upper bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
* '''ub''': [ ], Upper bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
* '''penub''': [ ], Penalty wt for upper bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
* '''penub''': [ ], Penalty wt for upper bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
* '''area''': [ ], Estimated peak area.
* '''area''': [ ], Estimated peak area.
* '''y''' = ''M''x''N'' measured responses with peaks to fit. Each row of (y) is fit to the peaks given in (peakdef).
* '''y''' = ''M''x''N'' measured responses with peaks to fit. Each row of (y) is fit to the peaks given in (peakdef).
====OPTIONAL INPUTS====
====OPTIONAL INPUTS====
* '''ax''' = 1x''N'' x-axis to fit to {default ax=1:N}.
* '''ax''' = 1x''N'' x-axis to fit to {default ax=1:N}.
* '''options''' = discussed below in the Options Section.
* '''options''' = discussed below in the Options Section.
====OUTPUTS====
====OUTPUTS====
* '''peakdefo''' = The input peak structure (peakdef) with parameters changed to correspond to the best fit values.
* '''peakdefo''' = The input peak structure (peakdef) with parameters changed to correspond to the best fit values.
* '''fval''' = Scalar value of the objective function evaluated at termination of FITPEAKS.
* '''fval''' = Scalar value of the objective function evaluated at termination of FITPEAKS.
* '''exitflag''' = Describes the exit condition (see LMOPTIMIZEBND).
* '''exitflag''' = Describes the exit condition (see LMOPTIMIZEBND).
* '''out''' = Structure array with information on the optimization/fitting (see LMOPTIMIZEBND).
* '''out''' = Structure array with information on the optimization/fitting (see LMOPTIMIZEBND).
* '''fit''' = Model fit of the peaks, i.e it is the best fit to (y).
* '''fit''' = Model fit of the peaks, i.e it is the best fit to (y).
* '''res''' = Residuals of fit of the peaks.
* '''res''' = Residuals of fit of the peaks.
===Algorithm===
===Algorithm===
Peaks are fit to the functions defined below based on the definitions in the field (peakdef.fun). The functions can be evaluated using independent functions or a wrapper function PEAKFUNCTION. See PEAKFUNCTION for more help.
Peaks are fit to the functions defined below based on the definitions in the field (peakdef.fun). The functions can be evaluated using independent functions or a wrapper function PEAKFUNCTION. See PEAKFUNCTION for more help.
For peakdef.fun = 'Gaussian' the function is
For peakdef.fun = 'Gaussian' the function is
    
    
where  is the  element of optional input (ax), and  corresponds to the peak parameters in the three-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are  and  .
where  is the  element of optional input (ax), and  corresponds to the peak parameters in the three-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are  and  .
For peakdef.fun = 'Lorentzian' the function is
For peakdef.fun = 'Lorentzian' the function is
   .
   .
Constraints that should be used are (bounds in peakdef) are  and  .
Constraints that should be used are (bounds in peakdef) are  and  .
For peakdef.fun = 'PVoigt1' the function is
For peakdef.fun = 'PVoigt1' the function is
    
    
where  corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are  and  , while  . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.
where  corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are  and  , while  . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.
For peakdef.fun = 'PVoigt2' the function is
For peakdef.fun = 'PVoigt2' the function is
    
    
where  corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used (bounds in peakdef) are  and  , while  . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.
where  corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used (bounds in peakdef) are  and  , while  . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.
A comparison of the four peaks is given in the figure below, and was generated using the following code:
A comparison of the four peaks is given in the figure below, and was generated using the following code:
:  ax    = 0:0.1:100;
:  ax    = 0:0.1:100;
:  y      = zeros(4,length(ax));
:  y      = zeros(4,length(ax));
:  plot(ax,peakgaussian([2 51 8],ax),'-b', ...
:  plot(ax,peakgaussian([2 51 8],ax),'-b', ...
:      ax,peaklorentzian([2 51 8],ax),'--k', ...
:      ax,peaklorentzian([2 51 8],ax),'--k', ...
:      ax,peakpvoigt1([2 51 8 0.5],ax),':g', ...
:      ax,peakpvoigt1([2 51 8 0.5],ax),':g', ...
:      ax,peakpvoigt2([2 51 8 0.5],ax),'-.r')
:      ax,peakpvoigt2([2 51 8 0.5],ax),'-.r')
:  legend('Gaussian','Lorentzian','PVoigt1','PVoigt2')
:  legend('Gaussian','Lorentzian','PVoigt1','PVoigt2')
   
   
===Options===
===Options===
*'''options'''  = structure array with the following fields:
*'''options'''  = structure array with the following fields:
* '''name''': 'options', name indicating that this is an options structure.
* '''name''': 'options', name indicating that this is an options structure.
* '''display''': [ 'off' | {'on'} ] governs level of display to the command window.
* '''display''': [ 'off' | {'on'} ] governs level of display to the command window.
* '''optimopts''': options structure from LMOPTIMIZEBND. This field is passed to LMOPTIMIZEBND and can be used to control the optimization / fitting.
* '''optimopts''': options structure from LMOPTIMIZEBND. This field is passed to LMOPTIMIZEBND and can be used to control the optimization / fitting.
===Examples===
===Examples===
:%Make a single known peak
:%Make a single known peak
:  ax            = 0:0.1:100;
:  ax            = 0:0.1:100;
:  y            = peakgaussian([2 51 8],ax);
:  y            = peakgaussian([2 51 8],ax);
:%Define first estimate and peak type
:%Define first estimate and peak type
:  peakdef      = peakstruct;
:  peakdef      = peakstruct;
:  peakdef.param = [0.1  43  5]; %coef, position, spread
:  peakdef.param = [0.1  43  5]; %coef, position, spread
:  peakdef.lb    = [0    0  0.0001]; %lower bounds on param
:  peakdef.lb    = [0    0  0.0001]; %lower bounds on param
:  peakdef.penlb = [1e-6 1e-6 1e-6];
:  peakdef.penlb = [1e-6 1e-6 1e-6];
:  peakdef.ub    = [10 99.9 40]; %upper bounds on params
:  peakdef.ub    = [10 99.9 40]; %upper bounds on params
:  peakdef.penub = [1e-6 1e-6 1e-6];
:  peakdef.penub = [1e-6 1e-6 1e-6];
:%Estimate fit and plot
:%Estimate fit and plot
:  yint  = peakfunction(peakdef,ax);
:  yint  = peakfunction(peakdef,ax);
:  [peakdef,fval,exitflag,out] = fitpeaks(peakdef,y,ax);
:  [peakdef,fval,exitflag,out] = fitpeaks(peakdef,y,ax);
:  yfit  = peakfunction(peakdef,ax); figure
:  yfit  = peakfunction(peakdef,ax); figure
:  plot(ax,yint,'m',ax,y,'b',ax,yfit,'r--')
:  plot(ax,yint,'m',ax,y,'b',ax,yfit,'r--')
:  legend('Initial','Actual','Fit')
:  legend('Initial','Actual','Fit')
===See Also===
===See Also===
[[peakfind]], [[lmoptimizebnd]], [[peakfunction]], [[peakgaussian]], [[peaklorentzian]], [[peakpvoigt1]], [[peakpvoigt2]], [[peakstruct]]
[[peakfind]], [[lmoptimizebnd]], [[peakfunction]], [[peakgaussian]], [[peaklorentzian]], [[peakpvoigt1]], [[peakpvoigt2]], [[peakstruct]]

Revision as of 15:25, 3 September 2008

Purpose

Peak fitting routine.

Synopsis

[peakdefo,fval,exitflag,out,fit,res] = fitpeaks(peakdef,y,ax,options)

Description

Based on the initial guess in input peakdef, FITPEAKS estimates the peak fit (also the Jacobian and Hessian), and makes a call to LMOPTIMIZEBND to find the best fit of the peaks to the data. (See LMOPTIMIZEBND for additional information.) Results are output to peakdefo.

Information about individual peaks is stored in standard peak structures (see PEAKSTRUCT). Information on multiple peaks is stored in a multi-record structure. Given a standard peak structure (peakdef) that contains an initial guess of peak locations and widths, FITPEAKS finds new parameters that best fits peaks to the rows of the data matrix (y). Results are output to a standard peak structure (peakdefo).

Fields of (peakdef) required in the initial guess for each peak are (.fun), (.param), (.lb), (.penlb), (.ub), and (.penub).

INPUTS

  • peakdef = multi-record standard peak structure with the following fields:
  • name: 'Peak', name indicating that this is a standard peak structure.
  • id: ' ', double or character string peak identification.
  • fun: [ {'Gaussian'} | 'Lorentzian' | 'PVoigt1' | 'PVoigt2' ], defines the peak function (see definitions in the Algorithm section).
  • param: Parameter list for each peak function. The number of parameters depends on the peak function used:
  • 'Gaussian': [height, location, width],
  • 'Lorentzian': [height, location, width],
  • 'PVoigt1': [height, location, width, fraction Gaussian],
  • 'PVoigt2': [height, location, width, fraction Gaussian].
  • lb: [ ], Lower bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
  • penlb: [ ], Penalty wt for lower bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
  • ub: [ ], Upper bounds on the function parameters. This is a row vector with the same number of elements as peakdef.param.
  • penub: [ ], Penalty wt for upper bounds, >=0. This is a row vector with the same number of elements as peakdef.param. If set to 0 this constraint is not employed.
  • area: [ ], Estimated peak area.
  • y = MxN measured responses with peaks to fit. Each row of (y) is fit to the peaks given in (peakdef).

OPTIONAL INPUTS

  • ax = 1xN x-axis to fit to {default ax=1:N}.
  • options = discussed below in the Options Section.

OUTPUTS

  • peakdefo = The input peak structure (peakdef) with parameters changed to correspond to the best fit values.
  • fval = Scalar value of the objective function evaluated at termination of FITPEAKS.
  • exitflag = Describes the exit condition (see LMOPTIMIZEBND).
  • out = Structure array with information on the optimization/fitting (see LMOPTIMIZEBND).
  • fit = Model fit of the peaks, i.e it is the best fit to (y).
  • res = Residuals of fit of the peaks.

Algorithm

Peaks are fit to the functions defined below based on the definitions in the field (peakdef.fun). The functions can be evaluated using independent functions or a wrapper function PEAKFUNCTION. See PEAKFUNCTION for more help.

For peakdef.fun = 'Gaussian' the function is


where is the element of optional input (ax), and corresponds to the peak parameters in the three-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are and .

For peakdef.fun = 'Lorentzian' the function is

  .

Constraints that should be used are (bounds in peakdef) are and .

For peakdef.fun = 'PVoigt1' the function is


where corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used are (bounds in peakdef) are and , while . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.

For peakdef.fun = 'PVoigt2' the function is


where corresponds to the peak parameters in the four-element vector (peakdef.param). Constraints that should be used (bounds in peakdef) are and , while . The Pseudo-Voigt peak shape is an estimate of the Gaussian and Lorentzian peak shapes convolved.

A comparison of the four peaks is given in the figure below, and was generated using the following code:

ax = 0:0.1:100;
y = zeros(4,length(ax));
plot(ax,peakgaussian([2 51 8],ax),'-b', ...
ax,peaklorentzian([2 51 8],ax),'--k', ...
ax,peakpvoigt1([2 51 8 0.5],ax),':g', ...
ax,peakpvoigt2([2 51 8 0.5],ax),'-.r')
legend('Gaussian','Lorentzian','PVoigt1','PVoigt2')


Options

  • options = structure array with the following fields:
  • name: 'options', name indicating that this is an options structure.
  • display: [ 'off' | {'on'} ] governs level of display to the command window.
  • optimopts: options structure from LMOPTIMIZEBND. This field is passed to LMOPTIMIZEBND and can be used to control the optimization / fitting.

Examples

%Make a single known peak
ax = 0:0.1:100;
y = peakgaussian([2 51 8],ax);
%Define first estimate and peak type
peakdef = peakstruct;
peakdef.param = [0.1 43 5]; %coef, position, spread
peakdef.lb = [0 0 0.0001]; %lower bounds on param
peakdef.penlb = [1e-6 1e-6 1e-6];
peakdef.ub = [10 99.9 40]; %upper bounds on params
peakdef.penub = [1e-6 1e-6 1e-6];
%Estimate fit and plot
yint = peakfunction(peakdef,ax);
[peakdef,fval,exitflag,out] = fitpeaks(peakdef,y,ax);
yfit = peakfunction(peakdef,ax); figure
plot(ax,yint,'m',ax,y,'b',ax,yfit,'r--')
legend('Initial','Actual','Fit')

See Also

peakfind, lmoptimizebnd, peakfunction, peakgaussian, peaklorentzian, peakpvoigt1, peakpvoigt2, peakstruct