Code covered by the BSD License  

Highlights from
peak fitting program for time-series signals

4.6087
4.6 | 24 ratings Rate this file 199 Downloads (last 30 days) File Size: 108 KB File ID: #23611
image thumbnail

peak fitting program for time-series signals

by

Tom O'Haver (view profile)

 

09 Apr 2009 (Updated )

Command-line peak fitter for time-series signals. Version 7.1, April 2015

| 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 28 different 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.

Acknowledgements

Interactive Peak Fitter (Version 11) inspired this file.

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

MATLAB release MATLAB 7.8 (R2009a)
MATLAB Search Path
/
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (104)
24 Apr 2015 Tom O'Haver

Tom O'Haver (view profile)

Nikola, first use the Matlab File > Import Data... command to import your spreadsheet data into Matlab's variable space. Then, once you have the independent and dependent variables (for example x and y) as vector variables or both together as a 2xn matrix, you can try the syntax example in the Help file (type "help peakfit") or on http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

Comment only
24 Apr 2015 Tom O'Haver

Tom O'Haver (view profile)

Victoria,

Example 15 of the internal examples shows how to fit data (the humps function in this case) with two Voigt profiles, flat baseline mode:
[FitResults,FitError]=peakfit(humps(0:.01:2),71,140,2,20,1.7,1,[31 4.7 90 8.8],3)

Example 25 shows how to fit that same data with two variable-alpha Voigt functions, shape 30.

x=[0:.005:1.3];y=humps(x);[FitResults,FitError]=peakfit([x' y'], 0, 0, 2, 30, 15, 10)

You will need version 7 or later of peakfit.m

Comment only
16 Mar 2015 Victoria

This seems like a great program, but unfortunately I am getting an error. I can fit my data to a single Voigt peak successfully, but it crashes when I ask for two peaks. Is there a quick fix for this? (it looks like "alpha" has size [1 2] in this call to voigt)

Error using /
Matrix dimensions must agree.
Error in peakfit>voigt (line 2128)
y = abs(xx-pos)/gV;
Error in peakfit>fitvoigt (line 2107)
A(:,j) = voigt(t,lambda(2*j-1),lambda(2*j),shapeconstant)';
Error in peakfit>@(lambda)(fitvoigt(lambda,xx,yy,extra)) (line 794)
TrialParameters=fminsearch(@(lambda)(fitvoigt(lambda,xx,yy,extra)),newstart,options);
Error in fminsearch (line 191)
fv(:,1) = funfcn(x,varargin{:});
Error in peakfit (line 794)
TrialParameters=fminsearch(@(lambda)(fitvoigt(lambda,xx,yy,extra)),newstart,options);
Error in RunVoigt2 (line 11)
fitres=peakfit([log10(q) I.*q.^2],0,0,2,20);

Comment only
13 Feb 2015 Nikola

Nikola (view profile)

Can any1 explain what i should put as input paramater? i have excel file of lifetime measruments of photoluminescence life time. and thats it...

Comment only
04 Feb 2015 Christian Stone

What lines must I alter? Sorry I am pretty bad

Comment only
02 Feb 2015 CKanellas

Is it possible for somebody to modify the code in order to fit Bézier curves? Or Better, is it in you plans?

Comment only
30 Jan 2015 CKanellas

wow!!

24 Jan 2015 Tom O'Haver

Tom O'Haver (view profile)

Any internal variable that want to "get out" can simply replace another output arguments in the first line of the function (FitResults, GOF, baseline, etc.) For example the residual data is the variable 'residual' and the fitted curve is 'bestmodel'. Just replace one of the existing output arguments that you don't need with the name of a variable that you do need.

Comment only
23 Jan 2015 Christian Stone

I am pretty new at this, how do I get the Residual data out to analyse? And the fitted curve data would also be useful
Thanks

Comment only
14 Dec 2014 xubiker  
28 Nov 2014 Tom O'Haver

Tom O'Haver (view profile)

Nazmul, the first step is to import your data into the form of Matlab variables, then download the peakfit.m function and place it in your Matlab path. The documentation for peakfit.m is on http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm, but you might want to start with ipf.m, which allows you to change parameters just by pressing single keystrokes, without typing a bunch of input arguments. See http://terpconnect.umd.edu/~toh/spectrum/ifpinstructions.html for animated instructions.

Comment only
28 Nov 2014 Nazmul Arefin

Tom, I am trying to fit multiple Gaussian peaks to my photoluminescence spectra. These are x-y data sets with luminescence in y-axis with wavelength in x-axis. Honestly, I don't have any strong background in Matlab scripting, I have been trying a lot to figure out how to apply your peakfit.m with my data. Is it possible you can suggest me where to start from? Many thanks for such a nice submission.

03 Oct 2014 Tom O'Haver

Tom O'Haver (view profile)

Thanks, Ismael. Both the command-line peakfit.m function and its interactive cousin ipf.m are described on http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm. There is no dead-trees publication. Thanks again.

Comment only
03 Oct 2014 Ismael

Ismael (view profile)

Hi Tom, First of all many thanks for this peak fitter, I have been using it quite a lot and it worked perfectly. Indeed I'm finishing one article where I have used it to fit some Raman spectra and I'm wondering how I can cite it, any clues?

01 Sep 2014 Tom O'Haver

Tom O'Haver (view profile)

Great! Thanks for testing it. I'll upload that version to the File Exchange.

Comment only
01 Sep 2014 Hao Wang

Tom, it works in my case. Thanks for your timely response and great work !

Comment only
28 Aug 2014 Hao Wang

Thanks Tom. I will experiment on it and give you some feedback asap.

Comment only
27 Aug 2014 Tom O'Haver

Tom O'Haver (view profile)

On the suggestion by Hao Wang, I've created a new version of peakfit.m (version 5.7) that includes a minimum width constraint as the 13th input argument. Download from http://terpconnect.umd.edu/~toh/spectrum/beta/peakfit.m
See example 19 in this file for an example of this feature. The default, when not specified, is equal to the independent variable (x) increment. Feedback please.

Comment only
24 Aug 2014 Tom O'Haver

Tom O'Haver (view profile)

Excellent idea, Hao Wang; I'll work on that.

Comment only
24 Aug 2014 Hao Wang

If the gaussian fit can have a lower bound that can help eliminate those very narrow fitted peaks. Is there a way to set such kind of a constraint (a lower or upper bound of the peak width) for any of the peak fit types? Other than that, works like a magic!

21 Aug 2014 Tom O'Haver

Tom O'Haver (view profile)

Ali, I can work on that. What's the specifications of your peaks: shape, height ratio, x-positions, and widths? How much noise? Are both peaks important in your work, or only the smaller one?

Comment only
21 Aug 2014 Ali

Ali (view profile)

Not smart enough to recognize real peaks. If you have very small peak before a very large peak, this function will fits the first one, not the second. Thanks for your hard work and sharing.

28 Jul 2014 Jade He

Thank you for your reply. I really want to thank you, but there is not that much I can do. I think I could share with you what I do with your code as a way to thank you.

I am working with MRI images. For my project, I need to calculate the lean muscle and fat within ROI. One method stated in "Distribution and Orientation of Bone in the Human Lumbar Vertebral Centrum" by T.S. Keller. This method requires Gaussian Fit in order to calculate the optimal value for thresholding.

Comment only
24 Jul 2014 Tom O'Haver

Tom O'Haver (view profile)

Yes, it's the full width at half maximum (at least for the Gaussian and Lorentzian shapes).

Comment only
24 Jul 2014 Jade He

Hello Tom,

I have a question about the output "width" of your code. Is this "width" the full width at half maximum?

Jade

Comment only
22 Jul 2014 Jade He

Thank you

22 Jun 2014 Ben

Ben (view profile)

Thanks Tom!

Comment only
21 Jun 2014 Tom O'Haver

Tom O'Haver (view profile)

Here's what I have in my library:

Fitting Experimental Data, Chris Brown, Randal Nelson

Non-Linear Regression, Samuel L.Baker

Probabilistic peak detection for first-order chromatographic data,
by Lopatka, Vivó-Truyols, and Sjerps

Modeling and Simulation: a Comprehensive and Integrative View, Tuncer I. Ören

Curve Fitting, the Reliability of
Inductive Inference, and the Error - Statistical Approach, Aris Spanos

Comparison of public peak detection algorithms for MALDI mass spectrometry data analysis, Chao Yang, Zengyou He and Weichuan Yu

Analysis of First-Derivative Based QRS Detection Algorithms, Natalia M. Arzeno, Zhi-De Deng, and Chi-Sang Poon,

You can find these and others online using Google search or Google Scholar.

Comment only
21 Jun 2014 Ben

Ben (view profile)

@Tom Very good function!

I was wondering whether you could recommend some review or classical papers for peak fitting and finding?

Thank you very much!

12 Jun 2014 Tobias

Tobias (view profile)

Hi Tom,

yes that works... Thanks a lot for your help!

Comment only
09 Jun 2014 Tom O'Haver

Tom O'Haver (view profile)

Try using a negative value for the Fano factor (the input argument called 'extra')

Comment only
09 Jun 2014 Tobias

Tobias (view profile)

Hi Tom,

you are perfectly right… Sorry for this, I simply never upgraded to version 5 before….

I tested the peak fitting (Raman). One thing I noticed is that the fitting curve falls steeper off towards the origin. This makes perfectly sense when the energy increases from the origin on.

In my case however, the energy becomes smaller and smaller (due to Raman scattering).

This means the BWF curve does not fit very well. In fact one can see that the steep falling slope should be on the other side.
(Unfortunately, I cannot add a screenshot here.)

Comment only
05 Jun 2014 Tom O'Haver

Tom O'Haver (view profile)

Ever since version 5, peakfit.m has been capable of fitting multiple peak types. Just enter the peak shape as a vector of peak shape numbers, and also use a vector for the "extra" variable. For example, to fit with a mode consisting of a Lorentzian plus a BWF peak with a Fano (asymmetry) factor of 5:

[FitResults,MeanFitError,Baseline]=peakfit([x;y],0,0,2,[2 15],[0 5],10,0,3)

Comment only
05 Jun 2014 Tobias

Tobias (view profile)

Hi Tom,

thanks a lot that was very fast!!

I will check for some data set to use data that allows to fit single BWF peaks.

At the moment however, I fit amorphous carbon which consists of two overlaying peaks. It can be fitted by a combination of two overlaying Gaussians and use the peak area (which works very well with your program) or, depending on the a-C nature it can be fitted using Lorentzian + BWF peak and use the peak height for further analysis.

Since it is a convolution of two differently shaped peaks I cannot compare both versions side by side for now.

I don't know how much work it would be for you to write the program such that I could use both peak types (Lorentzian+BWF).

If you have any time and inspiration to do that, I would be more than happy to compare both fitting options.

Again, many thanks for your help!
Toby

Comment only
05 Jun 2014 Tom O'Haver

Tom O'Haver (view profile)

OK, Toby, I have a beta version of peakfit.m and its interactive cousin ipf.m, both with the Breit-Wigner-Fano resonance peak (Shift-B key; Shape=15) replacing the bifurcated Lorentzian. They are located at http://terpconnect.umd.edu/~toh/spectrum/beta/. Let me know if this suits your needs and I'll upload it to the File Exchange.
Example:
>> x=1:250;
>> y=BWF(x,75,20,5)+.1.*BWF(x,175,20,5)+.00001*whitenoise(x);
>> [FitResults,MeanFitError,Baseline]=peakfit([x;y],0,0,2,15,5,10,0,3)

Comment only
04 Jun 2014 Tobias

Tobias (view profile)

Hi Tom,

great that you plan to implement the BWF peak!

Regarding the peak shape I found different approaches, one using the peak intensity itself

http://dx.doi.org/10.1103/PhysRevB.61.14095

and another one is provided by Origin, which uses a base line y0 in addition, but looks more similar to what you are suggesting

http://www.originlab.com/doc/Origin-Help/BWF-PAFunc

Probably, it depends on what is easier to implement.

Many thanks again,
Toby

Comment only
03 Jun 2014 Tom O'Haver

Tom O'Haver (view profile)

Excellent suggestion, Tobias. I'll add that peak shape to my program as soon as I get the chance. Let's make sure I have the right expression: y=((q*w/2+x-p).^2)./(((w/2).^2)+(x-p).^2) where q is the Fano parameter, w is the resonance line width, and p is the resonance peak energy, and x is energy (the independent variable).

Comment only
03 Jun 2014 Tobias

Tobias (view profile)

Hi Tom,

thanks a lot for the file! I use it on a regular basis to fit Raman peaks. Recently, I had to fit Raman peaks of amorphous carbon, which is sometimes done using a BWF (Breit-Wigner-Fano) peak. This peak shape however, is not included in your wonderful program.

Before I find time, understand all the connections and writing it myself, I was wondering if you have any plans to implement this peak?

Many thanks and best regards,
Toby

07 May 2014 Feng

Feng (view profile)

Hi Tom,

Thanks a lot for sharing this great work!

14 Feb 2014 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
08 Feb 2014 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
01 Feb 2014 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
01 Feb 2014 Alex

Alex (view profile)

@ 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

Tom O'Haver (view profile)

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.

Comment only
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

Comment only
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

Tom O'Haver (view profile)

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.

Comment only
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!

Comment only
29 Apr 2013 Ben

Ben (view profile)

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

Comment only
26 Apr 2013 Tom O'Haver

Tom O'Haver (view profile)

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?

Comment only
26 Apr 2013 Ben

Ben (view profile)

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

21 Feb 2013 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
18 Feb 2013 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
18 Feb 2013 Rems

Rems (view profile)

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.

Comment only
18 Feb 2013 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
17 Feb 2013 Rems

Rems (view profile)

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

Tom O'Haver (view profile)

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

Comment only
19 Sep 2012 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
07 Sep 2012 Himmst

Himmst (view profile)

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

Tom O'Haver (view profile)

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

Comment only
12 Jun 2012 Hussein

Thank you for the valuable information!

Comment only
12 Jun 2012 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
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?

Comment only
12 Jun 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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.

Comment only
12 Jun 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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.

Comment only
08 Jun 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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.

Comment only
07 Jun 2012 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
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.

Comment only
17 May 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
31 Mar 2012 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
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?

Comment only
07 Mar 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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

Comment only
20 Feb 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
20 Feb 2012 Kresten

Kresten (view profile)

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

07 Jan 2012 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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

AM (view profile)

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!

Comment only
04 Nov 2011 Tom O'Haver

Tom O'Haver (view profile)

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?

Comment only
04 Nov 2011 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
03 Nov 2011 AM

AM (view profile)

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

Comment only
01 Nov 2011 Sahrul Hidayat

nice program. thanks for sharing

12 Sep 2011 Sadik

Sadik (view profile)

Good job!

31 Aug 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
31 Aug 2011 mohammad

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

Comment only
10 Aug 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
25 Jul 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
24 Jul 2011 Tom O'Haver

Tom O'Haver (view profile)

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.

Comment only
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.

Comment only
26 Apr 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
19 Apr 2011 Tom O'Haver

Tom O'Haver (view profile)

No, no toolboxes required.

Comment only
19 Apr 2011 Victor Zhou

Is it require any toolbox to run the program?

Comment only
16 Apr 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
14 Apr 2011 Tom O'Haver

Tom O'Haver (view profile)

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

Comment only
19 Nov 2010 MB46

MB46 (view profile)

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

Joe K (view profile)

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

Tom O'Haver (view profile)

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

Comment only
21 Mar 2010 Jake

Jake (view profile)

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!

Comment only
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.

20 Jun 2014

Version 5.4: June, 2014. Replaces bifurcated Lorentzian with Breit-Wigner-Fano resonance peak (Shape=15).

29 Jul 2014

Version 5.5: July, 2014. Adds shape 21 (triangular) and 25 (Lognormal distribution). Moves peak table to lower panel and plots residuals as red dots.

02 Sep 2014

Version 5.7: August, 2014. Adds minimum width constraint as 13th input argument (See example 19); Can be a vector for multiple peak shapes. The default if not specified is the independent variable (x) interval.

24 Jan 2015

Adds an additional optional input argument, new peak shapes, and new output arguments.

24 Jan 2015

Adds an additional optional input argument, new peak shapes, and new output arguments.

24 Apr 2015

Version 7.1: April, 2015, adds peak shapes with three unconstrained iterated variables: 30=Voigt (variable alpha), 31=ExpGaussian (variable time constant), 32=Pearson (variable shape factor), 34=Gaussian/Lorentzian blend (variable percent).

24 Apr 2015

Version 7.1: March, 2015, adds peak shapes with three iterated variables: 30=voigt (variable alpha), 31=ExpGaussian (variable time constant), 32=Pearson (variable shape factor), 34=Gaussian/Lorentzian blend (variable percent). See Examples 25-28.

Contact us