Peak Fitter

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

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
  Gaussian, equal-width Gaussians, exponentially-broadened Gaussian,
  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 ([e-mail]). 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)

[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

[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?

[Ben]

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

[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.

[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.

[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.

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

[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

[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

[Tom O'Haver]

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

[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?

[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

[Tom O'Haver]

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

[Hussein]

Thank you for the valuable information!

[Tom O'Haver]

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

[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?

[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.

[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.

[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.

[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.

[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.

[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.

[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

[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.

[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.

[Tom O'Haver]

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

[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?

[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.

[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

[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.

[Kresten]

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

[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.

[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.

[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!

[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?

[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.

[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

[Sahrul Hidayat]

nice program. thanks for sharing

[Sadik]

Good job!

[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

[mohammad]

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

[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.

[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

[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.

[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.

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

[Tom O'Haver]

No, no toolboxes required.

[Victor Zhou]

Is it require any toolbox to run the program?

[Tom O'Haver]

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

[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

[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?

[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 ?

[Tom O'Haver]

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

[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!