Code covered by the BSD License

### Highlights fromPeak Fitter

4.3

4.3 | 10 ratings Rate this file 98 Downloads (last 30 days) File Size: 10.9 KB File ID: #23611

# Peak Fitter

09 Apr 2009 (Updated 21 Feb 2013)

Command-line peak fitter for time-series signals. Version 3.6: February, 2013.

File Information
Description

[FitResults,MeanFitError]=PEAKFIT(signal,center,window...)
A command-line peak fitting program for time-series signals,
written as a self-contained Matlab function in a single m-file.
Uses an non-linear optimization algorithm to decompose a complex,
overlapping-peak signal into its component parts. The objective
is to determine whether your signal can be represented as the sum of
fundamental underlying peaks shapes. Accepts signals of any length,
including those with non-integer and non-uniform x-values. Fits
Lorentzian, equal-width Lorentzians, Pearson, Logistic, exponential
pulse, abd sigmoid shapes (expandable to other shapes). This is a command
line version, usable from a remote terminal. It is capable of making
multiple trial fits with sightly different starting values and taking
the one with the lowest mean fit error.

T. C. O'Haver (toh@umd.edu). Version 3.6: February, 2013. Added fixed-position
Gaussian shape (16) and fixed-position Lorentzian shape (17).

Example 1:
>> x=[0:.1:10];y=exp(-(x-5).^2);peakfit([x' y'])
Fits exp(-x)^2 with a single Gaussian peak model.
ans =
1 5 1 1.665 1.7725
Peak number Peak position Height Width Peak area

Example 2:
x=[0:.1:10];y=exp(-(x-5).^2)+.1*randn(1,length(x));peakfit([x' y'])
Like Example 1, except that random noise is added to the y data.
ans =
1 5.0279 0.9272 1.7948 1.7716

Example 3:
x=[0:.1:10];y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.1*randn(1,length(x));
peakfit([x' y'],5,19,2,1,0,1)
Fits a noisy two-peak signal with a double Gaussian model (NumPeaks=2).
ans =
1 3.0001 0.49489 1.642 0.86504
2 4.9927 1.0016 1.6597 1.7696

Example 4:
>> x=[0:.005:1];y=humps(x);peakfit([x' y'],.3,.7,1,4,3);
Fits a portion of the humps function, 0.7 units wide and centered on
x=0.3, with a single (NumPeaks=1) Pearson function (peakshape=4)
with extra=3 (controls shape of Pearson function).

Example 5:
>> x=[0:.005:1];y=(humps(x)+humps(x-.13)).^3;smatrix=[x' y'];
>> [FitResults,MeanFitError]=peakfit(smatrix,.4,.7,2,1,0,10)
Creates a data matrix 'smatrix', fits a portion to a two-peak Gaussian
model, takes the best of 10 trials. Returns FitResults and MeanFitError.
FitResults =
1 0.4128 3.1114e+008 0.10448 3.4605e+007
2 0.3161 2.8671e+008 0.098862 3.0174e+007
MeanFitError =
0.68048

Example 6:
>> peakfit([x' y'],.4,.7,2,1,0,10,[.3 .1 .5 .1]);
As above, but specifies the first-guess position and width of the two
peaks, in the order [position1 width1 position2 width2]

Example 7:
[FitResults,LowestError,BestStart,xi,yi]=peakfit(smatrix,.4,.7,2,1,0,10)

As above, returns the vector x1 containing 100 interploated 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.

For more details, see
http://terpconnect.umd.edu/~toh/spectrum/CurveFittingC.html and
http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

Acknowledgements

Interactive Peak Fitter (Version 9.2) inspired this file.

This file inspired Demo Functions For Peak Detection And Fitting..

MATLAB release MATLAB 7.8 (R2009a)
Tags for This File
Everyone's Tags
Tags I've Applied
25 May 2013

Thanks, Avigdor,
The #9 exppulse function itself (starting on line 1152 in version 3.6) is a unit-height function, as are all the shape functions in peakfit.m. The amplitude (peak height) is supplied by the fitting function ("fitexppulse") beginning on line 1140; it's called PEAKHEIGHT. That's because the amplitude is a linear variable that is calculated by linear regression (in line 1148) rather than by iteration using fminsearch.

24 May 2013

Great script!
I was wondering if you could give an example of what is needed to modify an existing function. Specifically I am looking to add an amplitude variable to the #9 exppulse function.
Keep up the good work!

29 Apr 2013

Hi Tom,

What do you think some basic 2D shape, such as Gaussian, elongated Gaussian (Gabor), high peak and (approximately) pole with half sphere on the top?

I will send you an example data later. The data belong to my old project.

Ben

26 Apr 2013

Ben, not so far, but it's an interesting possibility. Can you give me an example of a 2D data set with peaks and the kind of information your would like to obtain from a 2D peak fitter program?

26 Apr 2013

Do you have a 2D version of this peak fitter?

21 Feb 2013

Version 3.6: February, 2013. Addition of fixed-position Gaussian (shape 16) and fixed-position Lorentzian (shape 17). See example 19 on http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm.

18 Feb 2013

So I guess It's about that that I add such a capability. I'll work on it and try to have something together in a few days.

18 Feb 2013

Thanks for the quick reply. I noticed in the meantime that this issue was also raised by Kresten (20 Feb 2012). I would like to work with only the position fixed (the non-trivial solution!). It would indeed be great to have a new peak type (gaussian) for which separate fixed parameters can be set for each peak.

18 Feb 2013

Thanks, Rems. Normally the easy way to add a new peak type to either peakfit.m or its interactive cousin ipf.m is to cannibalize one of the other shapes that you are not using and that has the same variable structure. Unfortunately none of my current peak types have separate fixed parameters for each peak, so it's a non-trivial job. Do you want the peak widths to be individually fit by the program or do you want all of the widths to be constrained to be equal? If all the widths were preset and fixed, and only the peak heights needed to be fit, then the problem is MUCH simpler and you use multilinear regression to obtain an instant solution (http://terpconnect.umd.edu/~toh/spectrum/CurveFittingB.html).

17 Feb 2013

Hi,
Best fitter available for Matlab, thanks for this wonderful work.
I am facing the following problem: I need to fit a signal with 4 gaussians but I want to fix the centre of each gaussian to a specific value. In others words, only the width and height of the gaussian need to be fitted. Of course I can choose the initial peak position and width as 'first guess' in argument 'start' but that does not guarantee that the final position will the be the same. Could you advise me a way to work with fixed position values? Many thanks

20 Sep 2012

Version 3.3: September, 2012. Added 11th input argument ('plots') to turn plotting OFF (plots=0) or on (plots=1); added Gaussian/Lorentzian blend (shape 13) bifurcated Gaussian (shape 14), and bifurcated Lorentzian (shape 15), whose shapes are controlled by 'extra'. See
http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

19 Sep 2012

That's a good suggestion, Himmst. I'll add an graphing ON/OFF switch to the next version. Thanks.

07 Sep 2012

First of all, thank you so much for sharing this valuable piece of code. An ON/OFF option for plotting would be nice. How can I switch off the plot? I am only interested in the fitted data and want to use it in my own plot. Your fit routine (in my case only interested in Gaussian) and my plotting would run in until-key-pressed loop for live-update of the datafit. As of now the built-in plot function always wants to take over my figure window. I tried to go through the hundreds of lines of code and delete everything that does plot, subplot or disp. But I still get an axes canvas. Any suggestions?

12 Jul 2012

These are very good script. Thank you very much. The scripts are very useful to help to solve my problem, and my headache

22 Jun 2012

Version 2.6: June, 2012. Added fixed-width Gaussian and Lorentzian peak shapes (shape numbers 11 and 12).

12 Jun 2012

Thank you for the valuable information!

12 Jun 2012

12 Jun 2012

I don't have repeated measurements,my application is as follows, I have only one measured spectrum of a part of mouse brain, and that spectrum is suspected to be an addition of several components. I know the spectrum of those components, and as you adviced me, I am trying with multilinear regression, to find the amplitude of each component so that if I add them I get the best fit to the main spectrum. In this way I should have an estimate of concentration or to know if this component is present or not.
do you advice me to continue with this approach, and maybe take several measurements to enhance the statistics?

12 Jun 2012

It means that the amplitudes of those components are so low, compared to the standard deviation of repeated measurements, that some of the measurements fall below zero. Remember, if the true result is really zero, then half of the individual measurements, on average, must be negative.

12 Jun 2012

but what does it mean when when a spectrum for a mixture of multiple fluorophores is decomposed into several individual spectra, and some of these spectra have negative amplitudes.

12 Jun 2012

Negative coefficients do have physical meaning, from a statistical point of view. They are an important part of the ensemble of results whose mean is zero (or very low) and whose standard deviation is relatively high. If the true average result is zero, the individual measurements will be negative about half the time. If you ignore or discard the negatives ones, you will bias the mean, which is never good.

12 Jun 2012

Thank you very much Tom ,multilinear regression worked fine,
but how to avoid getting negative coefficients? how to add constraints on the multilinear regression to get only coefficients with physical meaning, not just coefficients which results in a good fit.

08 Jun 2012

Version 2.5, June 2012, allows zeros as placeholders for unspecified input arguments.
e.g. peakfit([x y],0,0,0,0,0,0,0,2) to specify only the autozero mode. This makes it possible to control the autozero setting while leaving all the other input arguments at their default values.

08 Jun 2012

Using the Autozero function (default) tends to lead to zeros and Nan's being returned in some occasions. Commenting out the autozeroeing line fixed the problem but it would be nice to have it available.

07 Jun 2012

Hussein, for that purpose you can use the much simpler technique of multilinear regression. For example, see http://terpconnect.umd.edu/~toh/spectrum/CurveFittingB.html

07 Jun 2012

very useful program, but how can we predefine the position of the fitted peaks. this could be helpful for spectroscopy application where you can fit predetermined peaks or functions in a main plot that contain the contribution of all the sub peaks. in other words when you have a main plot, which has to be fitted using multiple peaks with known characteristics.

17 May 2012

Version 2.4: May, 2012, has modified
exponential broadening to using normal rather than circular convolution by zero-padding the FFTs. Thanks to Kalamaya for encouraging me to do this.

31 Mar 2012

Yes, you are correct. An error on my part. I'll attempt to add the zero padding. Thanks.

23 Mar 2012

@Tom O'Haver, Nice detailed function, I have yet to study it but look forward to it, thanks for sharing. One thing though, in your function 'ExpBroaden' you list:

% ExpBroaden(y,t) convolutes y by an exponential decay of time constant t
% by multiplying Fourier transforms and inverse transforming the result.
fy=fft(y);
a=exp(-[1:length(y)]./t);
fa=fft(a);
fy1=fy.*fa';
yb=real(ifft(fy1))./sum(a);

As it currently is coded, this will return the circular convolution, and not the 'usual' convolution. For a normal convolution, the result must be of length(a) + length(y) - 1, in which case the FFT's need to be zero padded for correct result. Otherwise, result is circular convolution. Was this your intent?

07 Mar 2012

Bogdan, if you could send some of your problematic data sets to me via email (address on any of my web pages), it would help. Thanks.

05 Mar 2012

Hi, Tom
I'm finding this code very useful, although I was wondering if you could look into a couple more things. Some of these instabilities may be cause by the fact that I have a truncated peak that I zero padded (see above):

1. In some datasets, I get NaN's and 0's if I try to fit for example 4 peaks, but 3 works ok. In the latter case, I tend to get peaks positioned at something like -4e004, whereas my data range is 0 to 300. I havenâ€™t been able to find a trend as to when this occurs.

2. When I run any fit other than Gaussian, I get NaN's and 0's for all peak parameters. [Results1,MeanFitError1,BestStart,xi,yi1]=peakfit([x1 y1],100,200,2,1); -- Is ok, but changing the "1" at the end causes issues.

3. If I let it choose Gaussian as a default, ( i.e. [Results1,MeanFitError1,BestStart,xi,yi1]=peakfit([x1 y1],100,200,2); ), it takes about 150 seconds to process my 300 data points. If I explicitly ask for a Gaussian shape (add a "1" as the last input), it takes seconds.
Thanks for the work

20 Feb 2012

That option does not exist in my program yet. If you want to fix BOTH the position and the width, then the problem is no longer an non-linear one, and you can obtain the best fit by multi-linear regression in a single statement. But the option to fix position OR width, leaving only one non-linear variable to be iterated, would be computationally easy, if I can figure out a way to do the UI. Thanks.

20 Feb 2012

Great program!
Can I fix a value like the gauss center value, width etc.

07 Jan 2012

Thanks for the careful evaluation and feedback, Bogdan. I'll look into the "close to the edge" problem and see if I can either bypass the problem or add zero padding as an option.

06 Jan 2012

Does a pretty good job. May require zero padding to work if one of your peaks is close to the edge of data. In my case it would only fit one peak no matter how many I specified until I zero padded around my data.
I briefly tried using window+center for the same effect but didn't get the same result.

04 Nov 2011

I read your last comment first, and check my linspace function and it does not have a default value of 100 points between values if non is specified. So I just put it in as 100 and it works! Not sure why my linspace is different. But I will also download the most current version as well.
Thanks!

04 Nov 2011

AM, I just checked it, and these examples are working just fine on my system (Matlab 7.8 R2009a). This error is thrown by the built-in "linspace" function, not by my peakfit.m function. I'm at a loss to understand what is going on on your system. Could there be that much difference between R2009a and R2009b?

04 Nov 2011

AM, I think you are using the wrong function. In the current version of my peakfit.m (version 2.2), line 23 is a comment line. There is no line reading "x = x2fx((0:n-2)');" in any of my code. Download a fresh copy of the current version and try again.

03 Nov 2011

With Matlab R2009b I am getting an error even when I try example 1.

??? Input argument "n" is undefined.
Error in ==> linspace at 23
x = x2fx((0:n-2)');
Error in ==> peakfit at 392
xxx=linspace(min(xx),max(xx));

I tried this on Matlab 2008b and got the same errors. Not sure what I am missing.
Thanks for any help!
Would love to try this

01 Nov 2011

nice program. thanks for sharing

12 Sep 2011

Good job!

31 Aug 2011

Load the .xls file into the Matlab workspace using File > Import... command. Then write a loop that goes through the matrix one row (or column) at a time and passes it to peakfit.m

31 Aug 2011

how can use this for a matrix file? ( .xls file)
thanks

10 Aug 2011

Version 1.8 takes AUTOZERO setting as 9th input argument (O=OFF, 1=ON), and it has improved data input argument flexibility.

25 Jul 2011

The latest version peakfit 1.6 removes the offending 'TypicalX', seemingly with no ill effects. There are also some other small changes to the format of the peak table on the figure. See http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

24 Jul 2011

Sorry about that! It works fine in R2009a. I'll check it out and see if I can come up with something that works in both versions.

24 Jul 2011

This seems to use a parameter no longer supported in Matlab R2010B.. I get this erorr message:
??? Error using ==> optimset at 204
Unrecognized parameter name 'TypicalX'. Please see the optimset reference page in the documentation for a list of acceptable option parameters. Link to
reference page.

26 Apr 2011

As suggested by MB46, the latest version (1.4) allows specifying that the widths of all gaussians be the same (select peak shape 6), and also the same for lorentzian peak fits (peak shape 7).

19 Apr 2011

No, no toolboxes required.

19 Apr 2011

Is it require any toolbox to run the program?

16 Apr 2011

I've added a demo script for this function. It's described on http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

14 Apr 2011

Joe K, I've just uploaded an updated version that includes the option to extract the data points (x,y) for the newly split peaks. See http://terpconnect.umd.edu/~toh/spectrum/peakfit.m

19 Nov 2010

This program is fantastic. However, I've run into a slight problem. When fitting multiple gaussians, I want to be able to put in a condition requiring the width of all the gaussians to be the same. Is there an easy way to do this? And indeed to restrict some of the other fitting parameters?

10 May 2010

Wonderful program. Thank you very much for posting it Tom. My question is, how do I extract the data points (x,y) for the newly split peaks ?

22 Mar 2010

Try turning your 'signal" matrix into its transpose, signal'

21 Mar 2010

Hi can someone help:
(I am a complete Matlab novice), I have my matrix 'signal' arranged in two columns, x ranges non uniformly from 0 to 1 and y starts at zero has a small peak on the start tail of the main peak (around10% size of main), I use peakfit(signal) and it throws up these errors:

??? Error using ==> optimset>checkfield at 310
Invalid value for OPTIONS parameter TypicalX: must be a matrix.

Error in ==> optimset at 248
checkfield(Names{j,:},arg,optimtbx);

Error in ==> peakfit at 155
options = optimset('TolX',.001,'TypicalX',center,'Display','off' );

and I have no idea why?

thanks for any help!

14 Apr 2009

Bug fixes. No major features added.

27 Apr 2009

Expanded description

13 Apr 2011

Added ability to access the x-y data for the individual fitted peaks, via the optional output parameters xi and yi.

25 Apr 2011

Added two new peakshapes: 6=Equal-width Gaussians, and 7=Equal-width Lorentzians.

04 May 2011

Expanded the description

25 Jul 2011

Several bug fixes. Reformatted peak table on graph.

10 Aug 2011

Version 1.8: Aug, 2011. Takes AUTOZERO setting as 9th input argument; improved data input argument flexibility.

24 Aug 2011

Bug fix

28 Sep 2011

Version 2.1: Sept, 2011. Accepts AUTOZERO 0 (none), 1 (linear), or 2 (quadratic).

20 Oct 2011

Version 2.2: October, 2011. Adds exponential pulse and sigmoid models

19 Jan 2012

Version 2.3: January, 2012. Bug fixes in background subtraction modes and
in handlng very small data sets.

19 Jan 2012

Version 2.3: January, 2012. Bug fixes in background subtraction modes and
in handlng very small data sets.

18 May 2012

Version 2.4: May, 2012, Exponential broadening uses normal rather than circular convolution.

18 May 2012

Version 2.4: May, 2012, Exponential broadening uses normal rather than circular convolution.

08 Jun 2012

Version 2.5: June, 2012, Allows zeros as placeholders for unspecified input arguments.

18 Jun 2012

Version 2.6: June, 2012. Added fixed-width Gaussian and Lorentzian peak shapes (shape numbers 11 and 12).

04 Sep 2012

Version 3: September, 2012, adds the ability to estimate the uncertainty of peak parameters using the bootstrap sampling method.

12 Sep 2012

Version 3.1: September, 2012. Unlimited Number of peaks. Bug fixes.

20 Sep 2012

Version 3.3: September, 2012. Added 11th input argument ('plots') to turn plotting OFF (plots=0) or on (plots=1); added Gaussian/Lorentzian blend, bifurcated (asymmetrical) Gaussian and Lorentzian.

15 Oct 2012

Version 3.4: October, 2012. Works in Matlab or Octave 3.6.1

02 Nov 2012

Version 3.4.2; Slight improvement in speed of exponentially-broadened shapes. Works in Matlab or Octave 3.6.1

15 Jan 2013

Version 3.51: January, 2013. Improved accuracy of linear autozero calculation. Improved calculation of default "start" guess when not specified in the input arguments.

21 Feb 2013

Version 3.6: February, 2013. Addition of fixed-position Gaussian shape (16) and fixed-position Lorentzian shape (17).