Peak Fitter (March 2013)

[Peakfit]  [Peak examples]  [Peakfit demos]  [Interactive Peak Fitter (ipf)]  [ipf demos]  [ipf instructions]  [Hints and Tips]   [Execution time]  [Have a question? Email me]

1. Command line function: peakfit.m (Version 3.6, February, 2013)  

Peakfit.m is a user-defined command window peak fitting function for Matlab or Octave, usable from a remote terminal.  It is written as a self-contained function in a single m-file. (To view of download, click Peakfit.m). It takes data in the form of a 2 x n matrix that has the independent variables (X-values) in row 1 and the dependent variables (Y-values) in row 2, or as a single independent variable vector.  You can also download a ZIP file containing peakfit.m, DemoPeakFit.m, ipf.m, Demoipf.m, and some sample data for testing.  For a discussion of the accuracy and precision of peak parameter measurement using peakfit.m, click here. Version 3.4 and above work in Octave as well as in Matlab. Version 3.6: February, 2013. Addition of fixed-position Gaussian shape (16) and fixed-position Lorentzian shape (17).

Peakfit can be called with several optional additional arguments. Note: In version 2.5 and later, all input arguments (except the signal itself) can be replaced by zeros to use their default values.

peakfit(signal);     
  Performs an iterative least-squares fit of a single Gaussian  peak to the entire data matrix "signal", which has x values in row 1 and Y values in row 2 (e.g. [x y]) or which may be a single signal vector (in which case the data points are plotted against their index numbers on the x axis).
 
peakfit(signal,center,window);       
  Fits a single Gaussian peak to a portion of the matrix "signal". The portion is centered on the x-value "center" and has width "window" (in x units). 

In this and in all following examples, set "center" and "window" both to 0 to fit the entire signal. 

peakfit(signal,center,window,NumPeaks);
  "NumPeaks" = number of peaks in the model (default is 1 if not specified; may be a large as desired in version 3.1 and above). 
  
peakfit(signal,center,window,NumPeaks,peakshape);
   Specifies the peak shape of the model: "peakshape" = 1-15. (1=Gaussian, 2=Lorentzian, 3=logistic, 4=Pearson, 5=exponentially broadened Gaussian; 6=equal-width Gaussians, 7=equal-width Lorentzians, 8=exponentially broadened equal-width Gaussians, 9=exponential pulse, 10=sigmoid, 11=fixed-width Gaussians, 12=fixed-width Lorentzians, 13=Gaussian/Lorentzian blend; 14=bifurcated Gaussian, 15=bifurcated Lorentzian; 16=Fixed-position Gaussians; 17=Fixed-position Lorentzians.
  
peakfit(signal,center,window,NumPeaks,peakshape,extra)
  Specifies the value of 'extra', used in the Pearson, exponentially-broadened Gaussian, Gaussian/Lorentzian blend, and bifurcated Gaussian and Lorentzian shapes to fine-tune the peak shape. 
  
peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials);
  Performs "NumTrials" trial fits and selects the best one (with lowest fitting error). NumTrials can be any positive integer (default is 1).
 
peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start)
  Specifies the first guesses vector "firstguess" for the peak positions and widths, e.g. start=[position1 width1 position2 width2 ...]

peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,autozero)
As above, but sets AUTOZERO mode in the last argument: AUTOZERO=0 does not subtract baseline from data segment;. AUTOZERO=1 (default) subtracts linear baseline from data segment; AUTOZERO=2, subtracts quadratic baseline from data segment.

peakfit(signal,0,0,0,0,0,0,0,2)
Use zeros as placeholders to use the default values of input arguments.  In this case, AUTOZERO is set to 2, but all others are the default values (Version 2.5 and later).

peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,autozero,fixedwidth,0)
11th input argument set to 0 (default is 1) to suppress plotting and command window printing.

Example 8: As above, returns the vector x1 containing 200 interpolated x-values for the model peaks and the matrix y1 containing the y values of each model peak at each xi.  Type plot(xi,yi(1,:)) to plot peak 1 or plot(xi,yi) to plot all peaks.
> [FitResults,LowestError,BestStart,xi,yi]= peakfit(smatrix,.4,.7,2,1,0,10)

Example 9:  Fitting single Gaussian on a linear background, using linear autozero mode (9th input argument = 1)
    >>x=[0:.1:10]';y=10-x+exp(-(x-5).^2);peakfit([x y],5,8,0,0,0,0,0,1)
  autozero=0 does not subtract baseline from data segment.
  autozero=1 (default) subtracts linear baseline from data segment.
  autozero=2, subtracts quadratic baseline from data segment.

Example 10: Fits a group of three peaks near x=2400 in DataMatrix3 with three equal-width exponentially-broadened Gaussians.
>> [FitResults,MeanFitError]= peakfit(DataMatrix3,2400,440,3,8,31,1)
FitResults =
            1       2300.5      0.82546       60.535       53.188
            2       2400.4      0.48312       60.535       31.131
            3       2500.6      0.84799       60.535       54.635
MeanFitError =
      0.19975

Example 11:   Example of an unstable fit to a signal consisting of two Gaussian peaks of equal height (1.0). The peaks are too highly overlapped for a stable fit, even though the fit error is small and the residuals are unstructured. Each time you re-generate this signal, it gives a different fit, with the peaks heights varying about 15% from signal to signal.  
>> x=[0:.1:10]'; y=exp(-(x-5.5).^2)+exp(-(x-4.5).^2)+.01*randn(size(x)); [FitResults,MeanFitError]=peakfit([x y],5,19,2,1)
FitResults =
            1       4.4059      0.80119       1.6347       1.3941
            2       5.3931       1.1606       1.7697       2.1864
MeanFitError =
      0.598
  Much more stable results can be obtained using the equal-width Gaussian model (peakfit([x y],5,19,2,6)), but that is justified only if the experiment is legitimately expected to yield peaks of equal width. See http://terpconnect.umd.edu/~toh/spectrum/CurveFittingC.html#Peak_width_constraints.

Example 12: Accuracy of autozero, for a single noiseless Gaussian on a curved background, specifying center (5.0) and window (5.5), but using placeholders (zeros) to accept the default values for NumPeaks, peakshape, extra, NumTrials, and start (Version 2.5 and above). The "true" values of peak parameters are 5.00, 1.00, 1.66, and 1.77.  Comparing autozero OFF (AUTOZERO=0), linear autozero (AUTOZERO=1), and quadratic autozero (AUTOZERO=2): 
>> x=[0:.1:10]';y=1./(1+x.^2)+exp(-(x-5).^2);peakfit([x y],5,5.5,0,0,0,0,0,0)
ans =
            1       4.9725       1.0179       1.8037       1.9538
>> x=[0:.1:10]';y=1./(1+x.^2)+exp(-(x-5).^2);peakfit([x y],5,5.5,0,0,0,0,0,1)
ans =
            1       5.0024      0.96197       1.5875       1.6256
>> x=[0:.1:10]';y=1./(1+x.^2)+exp(-(x-5).^2);peakfit([x y],5,5.5,0,0,0,0,0,2)
ans =
            1       5.0121       1.0007       1.6526       1.7602 
In this particular case, the quadratic autozero is more accurate for peak height, width, and area.

Example 13: Same as example 4, but with fixed-width Gaussian (shape 11), width=1.666. (added in version 2.6)
>> x=[0:.1:10];y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.1*randn(size(x));
>> [FitResults,MeanFitError]=peakfit([x' y'],0,0,2,11,0,0,0,0,1.666)
FitResults =
            1       3.9943      0.49537        1.666      0.87849
            2       5.9924      0.98612        1.666       1.7488

Example 14: Peak area measurements. Same as the example in Figure 15 on Integration and Peak Area Measurement.  All four peaks have the same theoretical peak area (1.772). The four peaks can be fit together in one fitting operation using a 4-peak Gaussian model, with only rough estimates of the first-guess positions and widths.  The peak areas thus measured are much more accurate than the perpendicular drop method:.
>> x=[0:.01:18];
>> y=exp(-(x-4).^2)+exp(-(x-9).^2)+exp(-(x-12).^2)+exp(-(x-13.7).^2);
>> peakfit([x;y],0,0,4,1,0,1,[4 2 9 2 12 2 14 2],0,0)

         Peak#    Position       Height       Width         Area
            1            4            1       1.6651       1.7725
            2            9            1       1.6651       1.7725
            3           12            1       1.6651       1.7724
            4         13.7            1       1.6651       1.7725

This works well even in the presence of substantial amounts of random noise:
>> x=[0:.01:18]; y=exp(-(x-4).^2)+exp(-(x-9).^2)+exp(-(x-12).^2)+exp(-(x-13.7).^2)+.1.*randn(size(x));
>> peakfit([x;y],0,0,4,1,0,1,[4 2 9 2 12 2 14 2],0,0)
ans =
            1       4.0086      0.98555       1.6693       1.7513
            2       9.0223       1.0007        1.669       1.7779
            3       11.997       1.0035       1.6556       1.7685
            4       13.701       1.0002       1.6505       1.7573

Sometimes experimental peaks are effected by exponential broadening, which does not by itself change the true peak areas, but does shift peak positions and increases peak width, overlap, and asymmetry, making it harder to separate the peaks.  Peakfit.m (and ipf.m) have an exponentially-broadened Gaussian peak shape (shape #5) that works well in those cases:, recovering the original peak positions, heights, and widths:
>> y1=expbroaden(y',-50);
>> peakfit([x;y1'],0,0,4,5,50,1,[4 2 9 2 12 2 14 2],0,0)
ans=
            1            4            1       1.6651       1.7725
            2            9            1       1.6651       1.7725
            3           12            1       1.6651       1.7725
            4         13.7      0.99999       1.6651       1.7716

Sometimes experimental peaks are effected by exponential broadening, which does not by itself change the true peak areas, but does shift peak positions and increases peak width, overlap, and asymmetry, making it harder to separate the peaks.  Peakfit.m (and ipf.m) have an exponentially-broadened Gaussian peak shape (shape #5) that works well in those cases:, recovering the original peak positions, heights, and widths:
>> y1=expbroaden(y',-50);
>> peakfit([x;y1'],0,0,4,5,50,1,[4 2 9 2 12 2 14 2],0,0)
ans=
            1            4            1       1.6651       1.7725
            2            9            1       1.6651       1.7725
            3           12            1       1.6651       1.7725
            4         13.7      0.99999       1.6651       1.7716

Example 15: Prints out a table of parameter error estimates; Version 3 only. See DemoPeakfitBootstrap for a self-contained demo of this function.
>> x=0:.05:9; y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.01*randn(1,length(x));
>> [FitResults,LowestError,BestStart,xi,yi,BootstrapErrors]= peakfit([x;y],0,0,2,6,0,1,0,0,0);

Peak #1         Position    Height       Width       Area
Bootstrap Mean: 2.9987      0.49717     1.6657      0.88151
Bootstrap STD:  0.0039508   0.0018756   0.0026267   0.0032657
Bootstrap IQR:  0.0054789   0.0027461   0.0032485   0.0044656
Percent RSD:    0.13175     0.37726     0.15769     0.37047
Percent IQR:    0.18271     0.55234     0.19502     0.50658
 
Peak #2         Position    Height       Width       Area
Bootstrap Mean: 4.9997     0.99466       1.6657      1.7636
Bootstrap STD:  0.001561   0.0014858     0.00262     0.0025372
Bootstrap IQR:  0.002143   0.0023511     0.00324     0.0035296
Percent RSD:    0.031241   0.14938       0.15769     0.14387
Percent IQR:    0.042875   0.23637       0.19502     0.20014

Demonstration script for peakfit.m

DemoPeakFit.m is a demonstration script for peakfit.m. It generates an overlapping Gaussian peak signal, adds normally-distributed noise, fits it with the peakfit.m function (in line 78), repeats this many times (NumRepeats in line 20), then compares the peak parameters (position, height, width, and area) of the measurements to their actual values and computes accuracy and precision (percent relative standard deviation). You can change any of the initial values in lines 13-30.  Here is a typical result for a two-peak signal with Gaussian peaks:

Peak shape: 1 (Gaussian)
Number of peaks: 2
Noise amplitude: 0.05
Number of data points in signal: 91
Number of repeat signals: 100
Number of trial fits per sample signal: 10
Actual parameters:
Position      Height       Width       Area
4.00           1           1.666      1.7735
6.00           0.5         1.666      0.88676
Average percent fitting error: 4.4847
Average measured parameters:
Position      Height       Width       Area
4.0019        1.0008       1.6681      1.777
6.0014        0.50191      1.66        0.88616
Percent RSDs of measured parameters:
Position      Height       Width       Area
0.71843       1.7824      4.0692      4.3941
0.96778       4.2114      8.0177      8.0543
Average %RSD of measured parameters for all peaks:
0.8431       2.9969      6.0435      6.2242
Percent errors of measured parameters:
Position      Height       Width       Area
0.048404     0.07906     0.12684     0.19799
0.023986     0.38235    -0.36158    -0.067655
Average Percent Parameter Error for all peaks:
0.036195      0.2307     0.24421     0.13282

In the above results, you can see that the accuracy and precision (%RSD) of the peak position measurements are typically the best, followed by peak height, and then the peak width and peak area, which are usually the worst.  

Which to use: iPeak or Peakfit?  Download these Matlab demos that compare iPeak.m with Peakfit.m for signals with a few peaks and signals with many peaks and that shows how to adjust iPeak to detect broad or narrow peaks. DemoPeakfitBootstrap demonstrates the ability of peakfit version 3 to computer estimates of the errors in the measured peak parameters. These are self-contained demos that include all required Matlab functions. Just place them in your path and click Run or type their name at the command prompt.  Or you can download all these demos together in idemos.zip.

2. Interactive keypress-operated version: ipf.m (Version 9.2, March, 2013)   

ipf.m is a version of the peak fitter for x,y data that uses keyboard commands and the mouse cursor. It is written as a self-contained Matlab function, in a single m-file. Version 9 has more flexible syntax and allows you to specify the initial focus by adding “center” and “window” values as additional input arguments, where 'center' is the desired x-value in the center of the upper window and “window” is the desired  width of that window.  For example:
   1 input argument:
    ipf(y) or ipf([x;y]) or ipf([x;y]');
   2 input arguments:
    ipf(x,y) or ipf([x;y],center) or ipf([x;y]',center);
   3 input arguments:
    ipf(x,y,center) or ipf(y,center,window) or
    ipf([x;y],center,window) or ipf([x;y]',center,window);
   4 input arguments:
    ipf(x,y,center,window);

0.91,0.18) focuses on second peak

In this example, pan and zoom controls are used to isolate a segment of a chromatogram that contains three weak peaks near 5.8 minutes The lower plot shows the whole chromatogram and the upper plot shows the segment.  Only the peaks in that segment are subjected to the fit. Pan and zoom are adjusted so that  the signal returns to the local baseline at the ends of the segment.

Pressing T turns on the autozero mode, causing the program to compute and subtract a linear baseline between those points.  Pressing 3, E, F performs a 3-peak exponentially-broadened Gaussian fit. The A and Z keys are then used to adjust the time constant ("Extra") to obtain the best fit (only used for the exponentially broadened and Pearson shapes).

Operating instructions for ipf.m (version 9.2)

1. At the command line, type ipf(x,y), (x = independent variable, y = dependent variable) or  ipf(datamatrix) where "datamatrix" is a matrix that has x values in row or column 1 and y values in row or column 2. Or if you have only one signal vector y, type ipf(y).  You may optionally add to additional numerical arguments:  ipf(x,y,center,window); where 'center' is the desired x-value in the center of the upper window and “window” is the desired  width of that window.
2. Use the four cursor arrow keys on the keyboard to pan and zoom the signal to isolate the peak or group of peaks that you want to fit in the upper window. (Use the < and > and ? and " keys for coarse pan and zoom and the square bracket keys [ and ] to nudge one point left and right). The curve fitting operation applies only to the segment of the signal shown in the top plot. The bottom plot shows the entire signal.  Try not to get any undesired peaks in the upper window or the program may try to fit them.
3. Press the number keys (1– 9) to choose the number of model peaks, that is, the minimum number of peaks that you think will suffice to fit this segment of the signal. For more than 9 peaks, press 0, type the number, and press Enter.
4. Select the desired model peak shape by pressing the '-' key and selecting the desired shape by number from the list that is displayed. Alternatively, you can select peak shapes directly in one keystroke by pressing the following keys: Gaussian (lower case g), Equal-width Gaussians (lower case h), Fixed-width Gaussians (Shift-G), fixed-position Gaussians (Shift-P), exponentially-broadened Gaussian (e), exponentially-broadened equal-width Gaussians (j); bifurcated Gaussian (Shift-H), Lorentzian (lower case L), Fixed-width Lorentzians (upper case L), fixed-position Lorentzians (Shift [ ), Equal-width Lorentzians (lower case ;), bifurcated Lorentzian (:), Logistic (lower case o), Pearson (lower case p), exponential pulse (u); sigmoid (s); and Gaussian/Lorentzian blend (`) models.  If the peak widths of each group of peaks is expected to be the same or nearly so, select the "equal-width" shapes. If the peak widths or peak positions are known from previous experiments, select the "fixed-width" or "fixed position" shapes. These more constrained fits are faster, easier, and much more stable than regular all-variable fits, especially if the number of model peaks is greater than 3 (because there are fewer variable parameters for the program to adjust - rather than a independent value for each peak). 
5. A set of vertical dashed lines are shown on the plot, one for each model peak. Try to fine-tune the Pan and Zoom keys so that the signal goes to the baseline at both ends of the upper plot and so that the peaks (or bumps) in the signal roughly line up with the vertical dashed lines, as show in the figures on the left. This does not have to be exact.
6. Press F to initiate the curve-fitting calculation. Each time you press F, another fit of the selected model to the data is performed with slightly different starting positions, so that you can judge the stability of the fit. Keep your eye on the residuals plot and on the "Error %" display. Do this several times, trying for the lowest error and the most unstructured random residuals plot. At any time, you can refine the signal region to be fit (step 2), the baseline position (step 9), change the number or peaks (step 3), or peak shape (step 4). 
7. The model parameters of the last fit are displayed in the upper window. For example, for a 3-peak fit:
   Peak# Position  Height   Width    Area
    1    5.33329   14.8274  0.262253 4.13361
    2    5.80253   26.825   0.326065 9.31117
    3    6.27707   22.1461  0.249248 5.87425
   The column are, left to right: the peak number, peak position, peak height, peak width, and the peak area. (Note: for exponential pulse (U) and sigmoid (S) shapes, Position and Width are replaced by Tau1 and Tau2). Press R to print this table out in the command window. Peaks are numbered from left to right. (The area of each component peak within the upper window is computed using the trapezoidal method and displayed after the width). Pressing Q prints out a report of settings and results in the command window, like so:

   Peak Shape = Gaussian
   Number of peaks = 3
   Fitted range = 5 - 6.64
   Percent Error = 7.4514    Elapsed time = 0.19741 Sec.
   Peak#  Position      Height       Width        Area
   1      5.33329      14.8274     0.262253      4.13361
   2      5.80253       26.825     0.326065      9.31117
   3      6.27707      22.1461     0.249248      5.87425

8. To select the Autozero mode, press the T key repeatedly; it cycles thorough Autozero OFF, Linear Autozero, and Quadratic Autozero.  When autozero is linear, a straight-line baseline connecting the two ends of the signal segment in the upper panel will be automatically subtracted. When autozero is quadratic, a parabolic baseline connecting the two ends of the signal segment in the upper panel will be automatically subtracted.  Use the quadratic autozero if the baseline is curved, as in these examples:

   Autozero OFF

   Linear autozero

   Quadratic autozero

9. If you prefer to set the baseline manually, press the B key, then click on the baseline to the LEFT the peak(s), then click on the baseline to the RIGHT the peak(s). The new baseline will be subtracted and the fit re-calculated. (The new baseline remains in effect until you use the pan or zoom controls). Alternatively, you may use the multi-point background correction for the entire signal: press the Backspace key, type in the desired number of background points and press the Enter key, then click on the baseline starting at the left of the lowest x-value and ending to the right of the highest x-value. Press the \ key to restore the previous background to start over.
10. If you want to manually specify the first-guess peak positions, press C, then click on your estimates of the peak positions in the upper graph, once for each peak. A fit is automatically performed after the last click. Peaks are numbered in the order clicked.  (The custom start positions remain in effect until you change the number of peaks or use the pan or zoom controls). It's possible to fit a peak whose maximum is outside the 
11. The A and Z keys control the "extra" parameter that is used only if you are using the variable-shape peak models: Pearson, exponentially-broadened Gaussian (ExpGaussian), bifurcated Gaussian, bifurcated Lorentzian, or Gaussian/Lorentzian blend. For the Pearson shape, an "extra" value of 1.0 gives a Lorentzian shape, a value of 2.0 gives a shape roughly half-way between a Lorentzian and a Gaussian, and a larger values give a nearly Gaussian shape. For the exponentially broadened Gaussian shapes, "extra" controls the exponential "time constant" (expressed as the number of points). For the Gaussian/Lorentzian blend and the bifurcated Gaussian and Lorentzian shapes, "extra" controls the peak asymmetry (a values of 50 gives a symmetrical peak). You can enter an initial value of "extra" by pressing Shift-X , typing in a value, and pressing the Enter key. Then you can adjust this value using the A and Z keys (hold down the Shift key to fine tune). Seek to minimize the Error % or set it to a previously-determined value. Note: if fitting multiple overlapping peaks with an "extra" parameter, it's easier to fit a single peak first, to get a rough value for the "extra" parameter, then just fine-tune that parameter for the multi-peak fit if necessary.
    
12. For difficult fits, it may help to press X to perform ten silent trial fits with slightly different first guesses and take the one with the  lowest fitting error. (You can change the number of trials,"NumTrials", in or near line 182 - the default value is 10). Remember: equal-width fits,  fixed-width fits, and fixed position shapes (added in version 9.1) are both faster, easier, and much more stable than regular variable fits, so use equal-width fits whenever the peak widths are expected to be equal or nearly so, or fixed-width (or fixed position) fits when the peak widths (or positions) are known from previous experiments. 
    
13. Press Y to display the entire signal full screen without cursors, with the last fit displayed in green. In version 5.6 and later, the residual is displayed in red, on the same y-axis scale as the entire signal.
14. Press M to switch the Y-axis of the full-signal plot in the bottom window between log and linear scales (added in version 7). This does not effect the fit, only the full-signal plot.
15. Press the D key to print out a table of model data in the command window: x, y1, y2, ..., where x is the column of x values of the fitted region and the y's are the y-values of each component, one for each peak in the model. You can then Copy and Paste this table into a spreadsheet or data plotting program of your choice.
16. Press W to print out the ipf.m function in the command window with the current values of  'center' and 'window' as input arguments. Useful when you want to return to that specific data segment later. Also prints out peakfit.m with all input arguments, including the last best-fit values of the first guess vector. You can copy and paste the peakfit.m function into your own code. 
17. Version 8 and later of ipf.m, and version 3 and later of peakfit.m, are able to estimate the expected variability of the peak position, height, width, and area from the signal, by using the bootstrap sampling method. This involves extracting 100 bootstrap samples from the signal, fitting each of those samples with the model, then computing the percent relative standard deviation (RSD) and the interquartile range (IQR) of the parameters of each peak. Basically this method calculates weighted fits of a single data set, using a different set of different weights for each sample. This process is computationally intensive can take several minutes to complete, especially if the number of peaks in the model and/or the number of points in the signal are high.  
    
    To activate this process in ipf.m, press the V key. It first asks you to type in the number of "best-of-x" trial fits per bootstrap sample (the default is 1, but you may need higher number here if the fits are occasionally unstable; try 5 or 10 here if the initial results give NaNs or wildly improbable numbers). (To activate this process in peakfit.m, you must use version 3.1 or later and include all six output arguments, e.g. [FitResults, LowestError, BestStart, xi, yi, BootstrapErrors]=peakfit...).  The results are displayed as a table in the command window. For example, for a three-peak fit (to the same three peaks used by the Demoipf demonstration script described in the next section) and using 10 as the number or trials:
    
    Number of fit trials per bootstrap sample (0 to cancel): 10
    Computing bootstrap sampling statistics....May take several minutes.
         
          Peak #1         Position   Height      Width      Area
          Bootstrap Mean: 800.5387   2.969539    31.0374   98.10405
          Bootstrap STD:  0.20336    0.02848      0.5061    1.2732
          Bootstrap IQR:  0.21933    0.027387     0.5218    1.1555
          Percent RSD:    0.025402   0.95908      1.6309    1.2978
          Percent IQR:    0.027398   0.92226      1.6812    1.1778
         
          Peak #2         Position   Height      Width      Area
          Bootstrap Mean: 850.0491   1.961368    36.4809     75.98684
          Bootstrap STD:  6.4458     0.14622      3.0657      4.7596
          Bootstrap IQR:  0.52358    0.058845     1.4303      3.9161
          Percent RSD:    0.75828    7.4549       8.4035      6.2637
          Percent IQR:    0.061594   3.0002       3.9205      5.1537
         
          Peak #3         Position   Height      Width      Area
          Bootstrap Mean: 896.23     1.029522   54.7959    59.87246
          Bootstrap STD:  6.9831     0.1342      4.1827     4.6316
          Bootstrap IQR:  1.9564     0.035681    3.1709     4.1454
          Percent RSD:    0.77916    13.035      7.6332     7.73578
          Percent IQR:    0.21829     3.4658     5.7868     6.9238
          Elapsed time is 98.394381 seconds.
          Min/Max Fitting Error of the 100 bootstrap samples: 3.0971/2.5747

Demoipf.m

Execution time

a) The execution time can vary over a factor of 4 or 5 or more between different computers, (e.g. a small laptop with 1.6 GHz, dual core Athlon CPU with 4 Gbytes RAM, vs a desktop with a 3.4 GHz i7 CPU with 16 Gbytes RAM).  Run the Matlab "bench.m" benchmark test to see how your computer stacks up.

b) The execution time increases with the square of the number of peaks in the model. (See PeakfitTimeTest.m).

c) The execution time varies greatly (sometimes by a factor of 100 or more) with the peak shape, with the exponentially-broadened Gaussian being the slowest and the fixed-width Gaussian being the fastest.  See PeakfitTimeTest2.m and PeakfitTimeTest2a.m  The equal-width and fixed-width shape variations are always faster.  

d) The execution time increases directly with NumTrials in peakfit.m. The "Best of 10 trials" function (X key in ipf.m) takes about 10 times longer than a single fit.  

3. Interactive version with mouse-controlled sliders (2006 version, for Matlab 6.5): InteractivePeakFitter.m

The opening screen view shows the entire signal below and the zoomed-in segment above. Pan and zoom the signal to select the desired group of peaks to be fit, adjusting so that the signal goes to the baseline at both ends of the upper plot and so that the peaks in the signal roughly line up with the vertical dashed lines. After selecting a speak shape and number of of peaks, press the F key to perform a fit.  A table of fit results are shown in the upper window, with the component  peaks shown in green. The baseline can be subtracted automatically or manually. The bottom graph shows the residuals and standard deviation of the error.  Pressing F again repeats the fit with slightly different start positions.  Adjusting the zoom pan and zoom settings, refining the background subtraction, and using custom start positions may help reduce the fitting error.   This figure shows a peak fit with an exponentially-broadened Gaussian shape.  The time constant of the exponential is controlled by the "Extra" slider. Other available peak shapes include Gaussian ('g'), Lorentzian ('l'), Logistic ('o'), and Pearson ('p'). Other shapes can be added if you know the mathematical form.

Operating Instructions for InteractivePeakFitter.m (2006 version, for Matlab 6.5)

1. Load your data set into the vectors x and y (x = independent variable, y = dependent variable) and then execute InteractivePeakFitter.m.
2. Use the Pan and Zoom sliders on the left and right to isolate the peak or group of peaks that you want to fit. You can fine-tune the Pan and Zoom by using the cursor arrow keys on the keyboard. The curve fitting operation applies only to the segment of the signal shown in the top plot. (The bottom plot shows the entire signal).  Try not to get any undesired peaks in the upper window or the program will try to fit them.
3. Use the # Peaks slider (or press the 1,2,3,4 or 5 keys) to choose the number of model peaks, that is, the minimum number of peaks that you think might suffice to fit this segment of the signal.
4. Use the Shape slider or press the G,L,O,P, or E keys to select the desired model peak shape. In this version there are five choices: Gaussian ('g'), Lorentzian ('l'), Logistic ('o'), Pearson ('p'), and Exponentially-broadened Gaussian ('e').
5. A set of vertical dashed lines are shown on the plot, one for each model peak. Try to fine-tune the Pan and Zoom using the sliders or the cursor arrow keys so that the signal goes to the baseline at both ends of the upper plot and so that the peaks (or bumps) in the signal roughly line up with the vertical dashed lines, as show in the figures on the left. This does not have to be exact.
6. Click on the Re-fit slider or press 'f' to initiate the curve-fitting calculation. Each time you click on this slider, another fit of the selected model to the data is performed. Keep your eye on the residuals plot and on the "Error %" display. Do this several times, trying for the lowest error and the most unstructured random residuals plot. At any time, you can refine the baseline position (step 8), change the number or peaks (step 3), or peak shape (step 4).
7. The model parameters of the last fit are displayed in the upper window. For example, for a 3-peak fit:
   

   1 5.33329 14.8274 0.262253 4.13361
   2 5.80253 26.825 0.326065 9.31117
   3 6.27707 22.1461 0.249248 5.87425

   The column are, left to right: the peak number, peak position, peak height, peak width, and peak area. (The area of each component peak within the upper window is computed using the trapezoidal method and displayed after the width). Peaks are numbered from left to right. To print this out in the Matlab command window, type "FitResults" in the Matlab command window and press Enter. FitResults(:,2) lists the peak positions, FitResults(:,3) lists the peak heights, FitResults(:,4) lists the peak widths, and FitResults(:,5) lists the peak areas falling within the upper window.
8. If you wish to set the baseline manually, press the 'B' key, then click on the baseline to the LEFT the peak(s), then click on the baseline to the RIGHT the peak(s). The new baseline will be subtracted and the fit re-calculated. (The new baseline remains in effect until you use the pan or zoom controls to change the fitting region; then you have to re-do the baseline for that new region).
9. If you want to manually specify the first-guess peak positions, click the Custom button or press 'c', then click on your estimates of the peak positions in the upper graph, once for each peak. A fit is automatically performed after the last click. Peaks are numbered in the order clicked.  (The custom start positions remain in effect until you change the number of peaks or use the pan or zoom controls).
10. The Extra slider (or the A and Z keys, in the keyboard-operated version  ipf.m) is used only if you are using the Pearson or the exponentially-broadened Gaussian (ExpGaussian) shapes. It controls the extra shape factor of the Pearson (a value of 1.0 gives a Lorentzian shape, a value of 2.0 gives a shape roughly half-way between a Lorentzian and a Gaussian, and a larger values give a nearly Gaussian shape) and it also controls the time constant of the exponentially broadened Gaussian. You can either adjust this to minimize the Error % or set it to a previously-determined value. You can change this slider's range to suit your proposes in lines 47-48 of InteractivePeakFitter.m.

DemoInteractivePeakFitter.m

>> ColumnLabels =
 Peak# Positon Height Width
ActualParameters =
 1 100 1 30
 2 130 2 40
 3 200 3 50
 4 400 1 30
 5 600 4 20
 6 630 1 20
 7 800 3 30
 8 850 2 40
 9 900 1 50
 10 1200 2 100

ZIP file with 2006 version, for Matlab 6.5 only.

Notes for all versions

Note 1: When a fit is performed, the lower half of the figure window displays the residuals (differences between the model and the data).  The "Error" or "Percent Error" reported by this program is the RMS difference between the best-fit model and the raw data over the fitted segment.

Note 2: When the number of peaks is 2 or greater, the green lines displayed in the upper plot after a fit is performed are the individual model peaks; the red line is the total model, i.e. the sum of the components, which is a least-squares best-fit to the blue data points.  

Note 3: In the slider version, if the range of the sliders is inappropriate for your signal, you can adjust the slider ranges in lines 45-52 of InteractivePeakFitter.m.

Note 4: The Extra slider (or the A and Z keys, in the keypress-operated version ipf.m) is used to manually control the global variable "extra", which could be used for a variety of purposes.  If you add your own peak shape (see "How to add a new peak shape", below) , you can use the extra variable for any purpose, is you need another adjustable parameter.

Note 5: To select the Autozero mode, press the T key repeatedly; it cycles thorough Autozero OFF, Linear Autozero, and Quadratic Autozero.  When autozero is linear, a straight-line baseline connecting the two ends of the signal segment in the upper panel will be automatically subtracted. When autozero is quadratic, a parabolic baseline connecting the two ends of the signal segment in the upper panel will be automatically subtracted.  Try to adjust the pan and zoom to include some of the baseline at the beginning and end of the segment in the upper window, allowing the automatic baseline subtract gets a good reading of the baseline  If that's not possible, you can still refine the baseline subtraction manually if you wish, by using the the BG button or press 'B' as described in step 8, above. To select the autozero mode, press T.  

Note 6: The peak areas are calculated only for the portion of the component peaks falling within the upper window.  If some of the component peaks are close enough to the edges of the window that they tail off outside the window, their peak areas will be inaccurate.

Note 7:  The "Error" or "Percent Error" reported by this program is simply the RMS difference between the best-fit model and the raw data over the fitted segment. It is not the error in the parameters derived from the fit (peak position, height, and width).  The accuracy of those parameters (the difference between the true and measured values of peak parameters) is much harder to estimate, because of course you don't know the true values of the peak parameters. However, the precision of measurement can be estimated if you can obtain several repeat samples of the same signal (with independent random noise), fit each sample using the same procedure, and calculate the standard deviation of the parameters from all those fits.  If the random errors arising from random noise in the signal or random variations is the accuracy of background subtraction are the dominant sources of error, then an estimate of the precision is a good indication of the accuracy of parameter measurement.  More on accuracy and precision of peak parameter measurement.

Hints and Tips

1. It's best not to smooth your data prior to curve fitting.  Smoothing can distort the signal shape and the noise distribution, making it harder to evaluate the fit by visual inspection of the residuals plot.  Smoothing your data beforehand makes it impossible to achieve the goal of obtaining a random unstructured residuals plot and it increases the chance that you will "fit the noise" rather that the actual signal.
2. The most important factor in non-linear iterative curve fitting is selecting the underlying model peak function, for example, Gaussian, Equal-width Guardians, Lorentzian, etc. (Some examples of this). It's worth spending some time finding and verifying a suitable function for your data. (See the instructions below for adding a new peak shape to the program). In particular, if the peak widths of each group of peaks is expected to be the same or nearly so, select the "equal-width" shapes;equal-width fits (available for the Gaussian and Lorentzian shapes in this version) are faster, easier, and much more stable than regular variable-width fits. 
3. You should always use the minimum number of peaks that adequately fits your data. (Some examples of this). Using too many peaks will result in an unstable fit - the green lines in the upper plot, representing the individual component peaks, will bounce around wildly for each repeated fit, without significantly reducing the Error %.  
4. This program uses an iterative non-linear search function (modified Simplex) to determine the peak positions and widths that best match the data. This requires first guesses for the peak positions and widths. (The peak heights don't require first guesses, because they are linear parameters; the program determines them by linear regression). The default first guesses for the peak positions are made by the computer on the basis of the pan and zoom slider settings and are indicated by the magenta vertical dashed lines. The first guesses for the peak widths are computed from the Zoom setting, so the best results will be obtained if you zoom in so that the particular group of peaks is isolated and spread out as suggested by the peak position markers (vertical dashed lines). 
5. If the peak components are very unevenly spaced, you might be better off entering the first-guess peak positions yourself by pressing the C key (clicking on the Custom slider) and then clicking on the top graph where you think the peaks might be. None of this has to be exact - they're just first guesses, but if they are too far off it can throw the search algorithm off.
6. Each time you perform a repeat fit (e.g. pressing the F key), the program adds small random deviations to the first guesses, in order to determine whether an improved fit might be obtained with slightly different first guesses. This is useful for determining the robustness or stability of the fit. If the error and the values of the peak parameters vary slightly in a tight region, this means that you have a robust fit (for that number of peaks). If the error and the values of the peak parameters bounces around wildly, it means the fit is not robust (try changing the number of peaks, peak shape, and the pan and zoom settings), or it may simply be that the data are not good enough to be fit with that model. The variability in the peak parameters from fit to fit is an estimate of the uncertainty caused by the curve fitting procedure (but not of the uncertainty caused by the noise in the data, because this is only for one sample of the data and noise).  Pressing F and R alternately will record the results of each successive fit in the command window.

How to add a new peak shape

1. To add a new peak shape to the five initial shapes, you must have a function for that peak shape, similar to gaussian.m, that takes as arguments a vector of x-values, a peak-center value, and a peak width, and returns a vector of that function evaluated at each value of x. Let's call it NewPeakShape. Add NewPeakShape to the PeakFitter folder (or if you are using the keypress version ipf.m, there is already a "NewPeakShape" function at the end of the code; just replace "g=sin(a.*x+b)" with your own function of x, a, and b, and assign your calculated function to the variable g).
2. Create a fitting function for this peak shape, similar to fitgaussian.m.

   function err = fitNewPeakShape(lambda,t,y)
   % Fitting function for NewPeakShape.
   global c
   A = zeros(length(t),round(length(lambda)/2));
   for j = 1:length(lambda)/2,
    A(:,j) = NewPeakShape(t,lambda(2*j-1),lambda(2*j))';
   end
   c = A\y';
   z = A*c;
   err = norm(z-y');

   Save this function as "fitNewPeakShape" into the PeakFitter folder (or if you are using the single-function keypress version ipf.m, there is already a "fitNewPeakShape" rear the end of the code; just modify its name to suit).
3. Change FitAndPlot.m to add a 6th case to the switch statement that starts at at line 20:

    case 6
    FitParameters=FMINSEARCH('fitNewPeakShape',start,options,x,y);
    ShapeString='NewPeakShape';

   and at to the one that starts at line 52:

    case 6
    A(m,:)=NewPeakShape(x,FitParameters(2*m-1),FitParameters(2*m));

   Note: If your function requires three adjustable arguments rather than two, you can assign one to the Extra slider, as is done here for the exponentially-broadened Gaussian and Pearson-type shape. In that case, add "extra" as the last argument of the function calls above, as is done in case 4 and 5 of the two switch statements in FitAndPlot.m. Change this slider's range to suit your proposes in lines 47-48 of InteractivePeakFitter.m.
4. Change ipfshape.m to add a 6th case to the switch statement:

    case 6
    ShapeString='NewPeakShape';

5. Change ipfpeaks.m to add a 6th case to the switch statement:

    case 6
    ShapeString='NewPeakShape';

6. Assign a key to the new peak shape and change KeyPressTest.m to add a case to the switch statement in KeyPressTest to interpret this keypress, similar to the other shapes.
   
7. Change line 39 in InteractivePeakFitter.m to read MaxShapes=6.

Number of unique visits since May 17, 2008: