Code covered by the BSD License  

Highlights from
Peak Fitter

4.5

4.5 | 14 ratings Rate this file 165 Downloads (last 30 days) File Size: 14.1 KB File ID: #23611
image thumbnail

Peak Fitter

by

 

09 Apr 2009 (Updated )

Command-line peak fitter for time-series signals. Version 5, February 2014

| Watch this File

File Information
Description

A command-line peak fitting program for time-series signals, written as a self-contained Matlab function in a single m-file. Uses a 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 any number of peaks of Gaussian, Lorentzian, equal-width Gaussian and Lorentzian, fixed-width Gaussian and Lorentzian, bifurcated Gaussian and Lorentzian, exponentially-broadened Gaussian, Pearson, Logistic, exponential pulse, sigmoid, Gaussian/Lorentzian blend, and Voigt 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, and it can estimate the standard deviation of peak parameters from a single signal using the bootstrap method.
Version 5: Feb, 2014. Adds multiple-shape models, designated by using a vector as the 5th input argument. See examples 17 and 18.

Acknowledgements

Interactive Peak Fitter (Version 10) inspired this file.

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

MATLAB release MATLAB 7.8 (R2009a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (65)
14 Feb 2014 Tom O'Haver

As suggested by Alex, I've uploaded version 5 of peakfit.m that includes the option to specify a different peak-shape for each peak, designated by using vectors rather than scalars as the 5th and 6th input arguments (peakshape and the 'extra' parameter for variable-shape types).

08 Feb 2014 Tom O'Haver

Prompted by Alex's suggestion, I am working on version 5 of peakfit.m that includes the option to specify a different line-shape for each peak. I'll upload it here when it's 'completely' debugged, but it mostly works now and you can download it from http://terpconnect.umd.edu/~toh/spectrum/peakfit.m
Bug reports welcome.

01 Feb 2014 Tom O'Haver

I have tried and failed so far to find a clean way to specify the peak shape of each peak in a group of multiple peaks. I'll look up "super-Lorentzian" and see what I can make of it. Thanks.

01 Feb 2014 Alex

@ Tom, thank you very much for sharing this very nice work! A few suggestions though. It would be great to have an option to specify line-shape for each peak to be fitted when fitting multiple peaks. Also, can you add Super-Lorentzian function?

11 Dec 2013 Tom O'Haver

Matthias,
Probably what's happening is that you are inadvertently triggering the bootstrap error estimation routine, which happens when you have 6 output arguments (line 787). To avoid this, modify the peakfit definition to read:
[FitResults,LowestError,bkgcoef]=peakfit(....
then call the modified peakfit function with just these 3 output arguments.

10 Dec 2013 Matthias

Hi Tom!

Thank you so much for this wonderful work, your function helped me a lot!

I'm using it to fit Raman spectra and it works great.
I have just one question:

Why is the speed of the function decreasing when I add the variable 'bkgcoeff' to the output values? I solved the problem by adding a 'save bkgcoef' to peak fit and loading it afterwards, but thats not a nice solution.

Best Regards
Matthias

08 Dec 2013 Fan D.Chen

This program is perfect, I have to say, thank you

22 Jun 2013 Spencer

A couple of errors in the post below:

sqrt(2 log(2)) = 1.177

sqrt(log(2)) = 0.832

22 Jun 2013 Spencer

Hi Tom-

This is a great tool.

I have a question about the constant '0.600516' that appears in the function 'gaussian'. In the comment for that function, you write "half-width=wid". If you define your gaussian as:

exp(-x^2/(2*sigma)^2)

then the half width = sqrt(2log(2))*sigma = 0.776*sigma. If you define your gaussian as:

exp(-x^2/sigma^2)

then the half width = sqrt(log(2))*sigma = 0.549*sigma.

It appears that you chose the latter definition but with "wid=full width at half max". That way wid = 2*sqrt(log(2))*sigma or alternatively sigma = 0.600516*wid.

In this case, the comment should read "FWHM=wid".

Thanks
Spencer

25 May 2013 Tom O'Haver

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 Avigdor

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 Ben

Hi Tom,

Thanks for your reply!

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.

Thanks in advance!
Ben

26 Apr 2013 Tom O'Haver

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 Ben

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

21 Feb 2013 Tom O'Haver

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 Tom O'Haver

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 Rems

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 Tom O'Haver

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 Rems

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 Tom O'Haver

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 Tom O'Haver

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

07 Sep 2012 Himmst

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 NZN NZN

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 Tom O'Haver

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

12 Jun 2012 Hussein

Thank you for the valuable information!

12 Jun 2012 Tom O'Haver

Hussein, please continue this discussion via email. I have sent a message to your email address as stated in your profile.

12 Jun 2012 Hussein

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 Tom O'Haver

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 Hussein

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 Tom O'Haver

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 Hussein

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 Tom O'Haver

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 Bogdan Dzyubak

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 Tom O'Haver

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 Hussein

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 Tom O'Haver

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 Tom O'Haver

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

23 Mar 2012 Kalamaya

@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:

function yb = ExpBroaden(y,t)
% 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 Tom O'Haver

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 Bogdan Dzyubak

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.
Looking forward to further updates,
Thanks for the work

20 Feb 2012 Tom O'Haver

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 Kresten

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

07 Jan 2012 Tom O'Haver

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 Bogdan Dzyubak

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 AM

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 Tom O'Haver

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 Tom O'Haver

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 AM

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 Sahrul Hidayat

nice program. thanks for sharing

12 Sep 2011 Sadik

Good job!

31 Aug 2011 Tom O'Haver

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 mohammad

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

10 Aug 2011 Tom O'Haver

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 Tom O'Haver

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 Tom O'Haver

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 John Bowen

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 Tom O'Haver

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 Tom O'Haver

No, no toolboxes required.

19 Apr 2011 Victor Zhou

Is it require any toolbox to run the program?

16 Apr 2011 Tom O'Haver

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

14 Apr 2011 Tom O'Haver

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 MB46

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 Joe K

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 Tom O'Haver

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

21 Mar 2010 Jake

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!

Updates
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).

13 Sep 2013

Version 3.9 adds exponentially broadened Lorentzian(peak shape 18); and alpha function (peak shape 19).

11 Nov 2013

Version 4.2 ,corrects some bugs, adds an additional autozero mode that subtracts a flat baseline without requiring that the signal return to the baseline at both ends of the signal segment, and adds a Voigt profile peak shape.

03 Feb 2014

Version 4.31: Jan, 2014. Adds 12th input argument, for + or +/- peak mode.

11 Feb 2014

Version 5: Feb, 2014. Adds multiple-shape models, designated by using a vector as the 5th input argument. See examples 17 and 18.

Contact us