Version 8.5 adds two new peak shapes and changes IQR (Interquatrile range) reported in the bootstrap method to
IQR/1.34896, which equals the relative standard deviation without outliers for a normal distribution.
Updated summary

Version 7.9: April, 2016, Added shapes 38=variable time constant ExpLorentzian; 40=sine wave; 41=rectangle;
42=flattened Gaussian; 43=Gompertz function (3 variable logistic); 44=1-exp(-k*x). "extra" defaults to 1 if not specified in input arguments.
Edited discription

August, 2015, Allows fixed width models (shapes 11 ands 12)
to have different fixed widths for each peak (specified by a vector);
also adds a fixed-height Gaussian model (shape 34) with heights specified
by a vector (input argument 10)

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.

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.

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

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.

Editor's Note: This file was selected as MATLAB Central Pick of the Week

peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,autozero,fixedparameters,plots,bipolar,minwidth,DELTA,clipheight)
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 44 different shapes, including models with multiple 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.

Cite As

Tom O'Haver (2021). peakfit.m (https://www.mathworks.com/matlabcentral/fileexchange/23611-peakfit-m), MATLAB Central File Exchange.
Retrieved .

MIng Jiang, I think what you want is covered on https://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm#findpeaksb. There are several functions described there that do just that: find successive multiple peaks and fit each to a model. Generally, these work best when the basic shape of all the peaks is the same.

Brilliant code! I would like to ask if we can use this code to perform fittings for lots of spectra automatically by just roughly assigning the peak positions to one of my spectra, i.e. it automatically recognises the where my peaks are located in each of my spectra and execute the fitting by the selected function? maybe this has been included in the latest version but I didn't find it......

Johannes,
Honestly, it's hard to choose between those two approaches. A fixed-position Voigt will be more complicated to code because it has a variable shape parameter as well as specific peak center positions. I've never tried that particular combination.

Dear Tom,
thank you very much for this great programm. Exactly what I was looking for and simple to use.I used voigt shapes to fit my nmr data and it worked very good.
In addition to the models already available, I would like to fit one NMR peak which in theory consisty of 3 Voigt peaks with different positions. Therefore, I am interested in a Voigt model using fix positions or even better fix positions in given areas. Even if this lowers the fit accuracy, the model / interpretation would be more meaningful in my case. The fitting algorithm should work like peakshape 16 or 17, but using Voigt profiles instead of lorentzian or gaussian respectively. I also noticed the model "peakshape 34 fix-width Voigt", which could be quite similar to my suggested model. Do you think it is easier to start from model 34 than from 16 or 17 to realize a "fix-position" Voigt model?
Thank you in advance for your help.
Best,
Johannes

This function is really well designed and allows numerous fits, however is it possible to implement a sinc-Gaussian convoluted function to adjust such a theoretical function to diffraction signal (sinc = sin(x)/x)?

I have an issue with running the script and wondered if you could possibly give me any advice/if it was a compatibility issue with more recent versions of MTEX/Matlab. I get the error 'Index exceeds the number of array elements (2).' on line 1388 of peakfit, which states 'Heights=PEAKHEIGHTS(2:1+NumPeaks);' As PEAKHEIGHTS is defined earlier on as length(NumPeaks) I can see why this error is happening, I'm just not sure how to fix it...
Any advice would be greatly appreciated!

Generally, when that happens, it's telling you that a small peak between 5 and 6 is not a good match for your data.

You could modify the fitting algorithm, the same way that the peak width is constrained to a minimum, but of course that will lead to much poorer fits with the peak position often ending up at the extreme of the range that you set.

Dear Tom
Thank you very much for this great script!

I was wondering if it was possible to set stricter criteria for the peak positions? For example, I am using shape 5 and I would really like that it finds a peak between x=5 and x=6 (it is a very small one). If I just enter position1 to be 5.5, I typically end up with peak1 being somewhere completely else - sometimes even outside of the data range. Any suggestions?

It depends on what you mean by "reciprocal of a gaussian". If you mean literally 1/gaussian, then you could code this as y=1./(exp(-((x-xo)./(0.60056120439323.*w)) .^2)) where xo is the peak center position and w is the full width at half maximum.

Dear Tom,
thank you for your answer. I was wrong then. Could you please advice what is the best solution to fit the reciprocal of a gaussian? Then, I will try to add it to the peakfit.m.
Anyway, thank you very much for your help.

Dear Tom,
awesome script!!!
I have some problems to understand the fitting with shape #5 (exponentionally broadened Gaussian).
Please, see the scipte below. What I am trying to do is fit the reciprocal of a gaussian to this shape, but I am lost in which value should be the "extra". If 0, the script takes forever, around 100 seesm that work, but the fitting is not that good. Could you please help me to improve the fitting?
Thank you very much.

Dear Tom,
Does peakfit allow to use constraints in general, similar to input 13 (minwidth), e.g. for max witdh, and - what I am searching for - to limit the peak center position inside the data range?

Shape= Lorentzian % Fitting Error= 2.4399% R2= 0.98718
Peak# Position Height Width Area
1 47.191 1793.8 0.51283 1377.3

Alternatively, if you must fit the entire signal in one go, you'll have to assign an off-screen peak to fit the baseline, then provide approximate first guesses for all 4 peaks.

Hello,
I would so happy to be able to use this peakfit function. But for now, i have totally failed despite spending a significant time to understand the different possibilities. I have a defined number of peaks (of typical shape pseudo-voigt = mix of Lorentzian and gaussian) and i already know quite well their position and width. So everything should be easy BUT since my background has some complex shape, i am unable to know how to manage this. My original peak always get lost since the program use one or two of them to fit the background.

please help me to understand how i can manage to fit those peaks?

Tom,
Sorry for the lack of clarity. I'm talking about the variable type double: https://fr.mathworks.com/help/matlab/ref/double.html.
As for the dimensions of PEAKHEIGHT: in the fitlorentzian function, there is a line which does "If BaselineMode==3,A=[ones(size(y))' A];end". The next instruction then assigns "PEAKHEIGTS=abs(A\y');end", which means that in general PEAKHEIGHT is a vector with only one value, but when BaselineMode is set to 3, PEAKHEIGHT is a vector storing two values.
I hope this is clearer now,
Loïc

Loic Toraille, I'd like to help but I don't understand your questions. What is a "double variable"? Also, what modification has been added to the dimensions of PeakHeights in the fitlorentzian function? Sorry. I just don't understand. Tom O'Haver

I'm slowly getting it to work, and I've encountered two problems so far:

- The text function which prints out the results seems to need double variables for the positions, whereas in my use of the code the input was single variables. To keep it simple, I added a double conversion just in the text function and it solves the problem, but I wonder why I'm the only one encountering it.

- I had some issues when alternating between different BaselineMode, (more specifically the any mode and the 3 mode I think), and after some digging I finally got to the core of the problem: when using BaselineMode 3, there is a modification added to the dimensions of PeakHeights in the fitlorentzian function. However, that modification was not being performed when I switched the BaselineMode back and forth, even though it worked fine using your ipf.m code. I finally understood that it was because the BaselineMode variable is called as a global variable in the fitlorentzian, but it cannot work since it is an input variable of the peakfit code. To solve this, I've added BaselineMode to the list of global variables of peakfit, changed the name of the input variable to InputBaselineMode, and declared at the beginning of peakfit that BaselineMode = InputBaselineMode.

Hope this helps, let me know if I did anything wrong which resulted in these problems appearing.

Hana, the problem might be caused by the global variables AA, xxx, PEAKHEIGHTS, FIXEDPARAMETERS, delta, BIPOLAR, and CLIPHEIGHT, Make sure they don't interfere with the variables in your GUI and rename them if necessary.

Hello, I'm running into an error. The first time I call peakfit in my GUI (making sure it has a figure of its own to work with) - it fails to run the bootstrapping, and it returns x and y variables that are of different length, which causes my program to crash.

If I run a separate script that does not use a GUI, it runs the bootstrapping fine, and when I run in the GUI, it will run for the rest of the day, with no problem, until I restart my Matlab. In the separate script, the code is exactly the same for the peakfit command... so I don't see what is being initialized there that is not working properly in my GUI on the first run. Thanks for any help. For now we will run a dummy program at the beginning of the day, but would be great if I could figure it out!

Adolfo, I'm sorry I missed this. We were traveling in Mongolia at the time and I was not checking in here. Better to contact me by email (toh@umd.edu). With this much offset in your baseline, it's best to subtract out most of it beforehand. Then:

Hello all and thankyou @Tom O'Haver for sharing this code. I have landed here looking for an automated way of fitting peaks in Laser-Induced Breakdown Spectroscopy (LIBS) spectra. I think that it should be straightforward, as the useful peaks I want to fit are well defined, with low noise, and isolated. The only problem is that some of them have a large baseline. But I have found that sometimes the fitting is very bad, and I don't what to try to improve. I know the expected central wavelength and width, and even with those starting values, the fit is not Ok.

If the only variable is the peak height, you don't need an iterative algorithm for that; you can just use the much simpler and faster classical least squares computation.

Hi Tom,
I am trying to figurate how can I modify the code to get a Voigt profile fit with known peak position, width and alpha and just have the height as the variable, similar to the recently added shape 50. Could you give me some guidelines? maybe be using shape 50 and modify the fitting function.
Thanks!

Alexandra, one way is to set bipolar=1 and use baseline mode 3.
For example, consider a single negative Lorentzian on a basline of 100:
x=1:100;
y=100-30*lorentzian(x,50,10);
datamatrix=[x;y];
[FitResults,FitError]=peakfit(datamatrix,0,0,1,2,1,1,0,3,0,1,1)

I'm trying to fit a Lorentzian to data that is offset in the positive-y direction, but which has a negative amplitude. Is there any way to specify negative amplitudes (I'm using peakshape = 2).

@Tom O'Haver, just to conclude on the harmonic oscillator topic, if you would like to understand what kind of data I was talking about, you can now check the paper I was working on when I first contacted you (link below). Figure 4 and Figure 7 are the sort of peaks I was trying to fit. Fitting that data to a Lorentzian profile was not physically correct, so I asked you about the idea of adding also the harmonic oscillator in your script. https://onlinelibrary.wiley.com/doi/full/10.1002/smll.201803774 Cheers,
A.

Teodoro,
1. You can cite my website at https://terpconnect.umd.edu/~toh/spectrum/ or my book "Pragmatic Introduction to Signal Processing 2019: Applications in scientific measurement", on Amazon at https://www.amazon.com/dp/1792867123 2. I don't see this behavior. What version of peakfit.m are you using? Can you send me a screen image or an example segment of data that demonstrates the problem? email: toh@umd.edu.

Dear Tom, This is a great contribution. Have two questions:
1. How to correctly cite this in a paper.
2. The version I have works well for fitting two Lorentzians to overlapped NMR signals; however, when I try to fit a third Lorentzian, it gives me wrong parameters for the third curve (on the Matlab workspace; surprisingly not on the Figure). Have you heared about this detail?
Finally, have a wonderful 2019!

@Tom O'Haver, following up on the harmonic oscillator, I don't know exactly what combination of parameters you're trying, but I think none of them can be <=0. That also includes all values of x, because you cannot have zero frequency or negative frequency. The condition that c is approximately equal to the FWHM only holds when c << b. That PDF form Harvard on the harmonic oscillator helps you understand that the equation is valid for a special case of damped/driven harmonic oscillator, and contains a few indications of the boundary conditions, such as the one I just mentioned. Meanwhile, I wrote my own fitting code for the harmonic oscillator, using the fitnlm function from Matlab. I could get the standard error from there as well, so my problem is solved. But I still think that this equation would be a very useful addition to your script. Cheers.

Attempting to fit the entire signal in one fitting operation fails, because the least-squares method tries to minimize the overall residual error for the whole signal. Since the left-hand two peaks are so much weaker than the main peak on the right, the algorithm ignores the two smaller peaks because their contribution to the total error is insignificant compared to the two larger peaks. In other words, it can get a lower total error by adjusting the fit to the larger peaks only.

The solution to this is to fit the signal in two separate sections. I used the interactive peak fitter ipf.m to zoom in to different sections and try various fits before settling on these two (using the W key to print out the peakfit syntax. Note that in the first fit, a third model peak is added just to account for the background cased bt the off-screen peak to its right.

Dear Tom... I´m uploading both figure and the data. i would like to "impose" 4 peaks (two at the left hand side, and two inside the main peak). I´ve tried many parameters without sucess. If you could help me, would be great.

It does look like your data should be well fit by a 3-Gaussian model using peakfit.m. I can not demonstrate it because my function fits numerical data, not graphical data.

I´m working on a NMR experiment and the answer can be understood as a Gaussian mixture. So, I´m trying to use your code to decompose the signal. The answer I´m getting can be seem here: https://www.dropbox.com/s/2qsemu0rdzcesu0/Capture.PNG?dl=0

I´m using a "trying and error" approach to fit the components and I would like to know if there is a way to fine tunning the peaks and interval. Note some peak shifts (mostly in the comp.2) and the tail on the right hand side.

Andre Sartori, for me, that expression, y=a/(((b*2*pi)^2-(x*2*pi)^2)^2+(c*x*4*pi*pi)^2)^0.5, or in Matlab syntax, y=a./(((b.*2.*pi).^2-(x.*2.*pi).^2).^2+(c.*x*4.*pi.*pi).^2).^0.5;, results in a variety of shapes depending on the values of b, and c, including double peak shapes that are symmetrical on x=0. Moreover, c is the FWHM only if C=1. Is that what I should expect? I'm just not sure I have the expression right.

This is the equation that I find useful, i.e. with the frequency in Hz instead of radians, because typically the measurements are done that way, and also presented that way.

Then the output becomes directly comparable to the output of the Gaussian or Lorentzian profiles, with c being the FWHM (in Hz) in case the oscillator is under-damped (which is up to the user to know), b being the peak center (in Hz) and a equal to F0/m which is not quite a height, but which can be treated as such.

Then the calculation of the important Q-factor is simply b/c. If I get the std of these parameters, with simple error propagation I have the uncertainty of Q, which is how I currently do with the Lorentzian.

@Tom O'Haver, you can check the equation for the amplitude, A(omega) here: http://ipl.physics.harvard.edu/wp-uploads/2013/03/15c_s07_4.pdf I guess you can probably treat F0/m as constant a, omega_zero as constant b, and gamma as constant c.
But I don't think there will be a problem with infinite solutions due to 2*pi radians, because there's no sin or cos function in the equation. It's not periodic. The shape of the equation is similar to a Lorentzian profile, centered at omega_zero.

Right now I'm working on some resonance measurements of MEMS devices, and I've been applying the Lorentzian profile to fit the peaks and obtain the so-called Q-factor, which is the peak center divided by the FWHM, but this is not the correct curve from a physics point-of-view, because the devices resonate according to the damped/driven harmonic oscillator. The shapes don't quite match and the Q-factors are not quite right. So I need to re-fit everything to the harmonic oscillator equation. But my problem is that I can't easily fit the equation and at the same time get the standard deviation, which is a very handy feature of your peakfit function. I want to have a number for the Q-factor and a plus/minus something.

Good idea, Andre Sartori. I'll give it a try. Presumably there would be three nonlinear variables, the frequency and phase of the oscillation and the exponential decay constant, right? I anticipate that there might be a problem with the phase variable because there should be an infinite number of solutions with phases spaced 2*pi radians (360 degrees) apart.

Calculating the precision of the peak parameters:
[FitResults, GOF, baseline, coeff, residuals, xi, yi, BootstrapErrors] = peakfit([x;y],0,0,2,6,0,1,0,0,0);
Displays parameter error estimates by the bootstrap method; See DemoPeakfitBootstrap for a self-contained demo of this function.

I have a problem when I try to fit a fixed position gaussian peak and 2 lorentzian peak (variable) together.
My call is this: [FitResults,FitError]=peakfit(M1, ... % datamatrix
521, ... % center
100, ... % window
4, ... % Number of peaks
[16 2 1 ], ... % peakshape
[1 0 0], ... % extra
[10 1 1], ... % numtrials
[521 3.0 505 3.0 517 3.0 ], ... % first guess
0, ... % autozero
[521] ... % fixedparameters
)

The problem is, the 521 isn't fixed at all. May you have a hint for me.

Also whe I have the Peakshape in another order like [1 16 2] it will not fit at all. It just give me this back:
Index exceeds matrix dimensions.
Error in peakfit>fitmultiple (line 3539)
A(j,:)=gaussian(xx,FIXEDPARAMETERS(j),lambda(j));
Error in peakfit>@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)) (line 985)
TrialParameters=fminsearch(@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)),newstart,options);
Error in fminsearch (line 189)
fv(:,1) = funfcn(x,varargin{:});
Error in peakfit (line 985)
TrialParameters=fminsearch(@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)),newstart,options);
Error in excel1 (line 53)
[FitResults,FitError]=peakfit(M1, ... % datamatrix
>>
May you can help me, cause I need it to fit a huge set of data.

desert fish, the 4th column of the FitResults matrix returned by the peakfit function is the full width at half maximum. The left and right bounds depend on how you define them: e.g. at 1% of peak height, 0.1% or peak height, etc. You would have to calculate that yourself depending on your choice. See https://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm#PeakStartAndEnd

Hi Tom, thanks for sharing this great tool. Is peakfit a defined function? I got the error: Undefined function or variable 'peakfit'. Could you please tell me how to define it？ Thanks.

I finally manage.
In order to retrieve data (array) values from the peakfit plot, follow these steps:

1) Run the script,
2) Type the following commands:
h = findobj(gca,'Type','line'); x=get(h,'Xdata');
y=get(h,'Ydata');
3) Since we deal with several curves, to plot a unique curve you have to use indexing. For instance for the first green curve:
plot(x{1},y{1})

Dear Tom,
First, thanks for this amazing tool.
I'm not smart enough to get rid of a (simple?) problem.
Reading the comments, I can plot the green curves. Right.
However, I have no clue about getting the (xi) & (yi) values (same dimensions) to continue the game ! Indeed, the matlab plot function allows plotting a curve from a scalar (xi) and a matrix (yi). Is there any way to get the related xi & yi points ?!
This might be really useful for coming users...
Thanks for your consideration.

Tom O'Haver
That so very good day.
Let me turn to you asking for your kind help in order to get the code in matlab for the fourier transform of continuous-time (FFT) directly and inversely, since my knowledge in matlab are not very good.
I appreciate who I can help and collaborate.
Excuse my bad English and my native language is Spanish

Yes. Just call peakfit with the output arguments [FitResults,GOF,baseline,coeff,residual,xi,yi],
then plot(xi,yi) will plot the model components separately (the "green" line) and plot(x,y-residual) will plot the combined model (the "red" line). You can make then whatever color and style you want.

Dear Tom,
is there way to simply extract the curves of the models?
So in principle, how can I plot the green and red line in a separate graph?
Thanks a lot! Great tool!

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

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.

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

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.

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.

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?

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.

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!

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?

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.

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.

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

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:

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.

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)

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

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?

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

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.

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.

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

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.

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.

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.

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.

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!

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?

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.

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

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

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

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?

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?

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.

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.

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.

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.

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.

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.

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.

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

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

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.

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.

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.

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!

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?

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.

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

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.

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

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?

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' );

You can also select a web site from the following list:

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Alex SchetininIt's just a wonderful, great piece of work. It's a pleasure to use and explore it. Thanks a lot, Tom!

Tom O'HaverMIng Jiang, I think what you want is covered on https://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm#findpeaksb. There are several functions described there that do just that: find successive multiple peaks and fit each to a model. Generally, these work best when the basic shape of all the peaks is the same.

Ming JiangBrilliant code! I would like to ask if we can use this code to perform fittings for lots of spectra automatically by just roughly assigning the peak positions to one of my spectra, i.e. it automatically recognises the where my peaks are located in each of my spectra and execute the fitting by the selected function? maybe this has been included in the latest version but I didn't find it......

Tom O'HaverJohannes,

Honestly, it's hard to choose between those two approaches. A fixed-position Voigt will be more complicated to code because it has a variable shape parameter as well as specific peak center positions. I've never tried that particular combination.

Johannes ThienenkampDear Tom,

thank you very much for this great programm. Exactly what I was looking for and simple to use.I used voigt shapes to fit my nmr data and it worked very good.

In addition to the models already available, I would like to fit one NMR peak which in theory consisty of 3 Voigt peaks with different positions. Therefore, I am interested in a Voigt model using fix positions or even better fix positions in given areas. Even if this lowers the fit accuracy, the model / interpretation would be more meaningful in my case. The fitting algorithm should work like peakshape 16 or 17, but using Voigt profiles instead of lorentzian or gaussian respectively. I also noticed the model "peakshape 34 fix-width Voigt", which could be quite similar to my suggested model. Do you think it is easier to start from model 34 than from 16 or 17 to realize a "fix-position" Voigt model?

Thank you in advance for your help.

Best,

Johannes

arnoldThanks Tom, that was easier than I thought. This function is so powerful and simple to use... kudos once more!

Tom O'HaverJust adjust "center" and/or "window" so that the truncated regions fall outside the window.

arnoldGreat function. While I've used this for years, I have on question: Is there a way to fit truncated peaks that are cut off on either side?

명욱 박Tom O'Haverstéphane béchu, sorry for the delay - I can perhaps in this case blame the Christmas holidays.

But if you look at https://terpconnect.umd.edu/~toh/spectrum/peakfitsinc.m, you will find a modified version of peakfit with the sinc funcitons replacing the triangular funciton (shape 21) :

example:

Create sinc funciton data vector

>> x=0:.001:20; w=1; p=10; x=0:.1:20;

>> y=sincfunction(x,p,w));

Fit sinc funciton to x,y

>> peakfitsinc([x;y],0,0,1,21)

Results:

1 10 1 1 0.97979

Instructions are on page 410 of mybook.

Tom O'HaverStéphane, that's a good idea; I'll try to implement that when I have the time.

Tom

stéphane béchuThis function is really well designed and allows numerous fits, however is it possible to implement a sinc-Gaussian convoluted function to adjust such a theoretical function to diffraction signal (sinc = sin(x)/x)?

Tobias Winckler-CarlsenCatherine GoddardHi again!

I've actually managed to get the code working now so please ignore my last message!

Best,

Rellie

Catherine GoddardDear Tom,

I have an issue with running the script and wondered if you could possibly give me any advice/if it was a compatibility issue with more recent versions of MTEX/Matlab. I get the error 'Index exceeds the number of array elements (2).' on line 1388 of peakfit, which states 'Heights=PEAKHEIGHTS(2:1+NumPeaks);' As PEAKHEIGHTS is defined earlier on as length(NumPeaks) I can see why this error is happening, I'm just not sure how to fix it...

Any advice would be greatly appreciated!

Best,

Rellie Goddard

Tom O'HaverGenerally, when that happens, it's telling you that a small peak between 5 and 6 is not a good match for your data.

You could modify the fitting algorithm, the same way that the peak width is constrained to a minimum, but of course that will lead to much poorer fits with the peak position often ending up at the extreme of the range that you set.

Dragana RistanovicDear Tom

Thank you very much for this great script!

I was wondering if it was possible to set stricter criteria for the peak positions? For example, I am using shape 5 and I would really like that it finds a peak between x=5 and x=6 (it is a very small one). If I just enter position1 to be 5.5, I typically end up with peak1 being somewhere completely else - sometimes even outside of the data range. Any suggestions?

Tom O'HaverIt depends on what you mean by "reciprocal of a gaussian". If you mean literally 1/gaussian, then you could code this as y=1./(exp(-((x-xo)./(0.60056120439323.*w)) .^2)) where xo is the peak center position and w is the full width at half maximum.

Esteban Pedrueza VillalmanzoDear Tom,

thank you for your answer. I was wrong then. Could you please advice what is the best solution to fit the reciprocal of a gaussian? Then, I will try to add it to the peakfit.m.

Anyway, thank you very much for your help.

Tom O'HaverEsteban Pedrueza Villalmanzo,

The exponentially broadened Gaussian will not fit the reciprocal of a Gaussian. That is not one of the available shapes, but you can add that as a new shape yourself, For instructions, see https://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm#NewShape

Esteban Pedrueza VillalmanzoDear Tom,

awesome script!!!

I have some problems to understand the fitting with shape #5 (exponentionally broadened Gaussian).

Please, see the scipte below. What I am trying to do is fit the reciprocal of a gaussian to this shape, but I am lost in which value should be the "extra". If 0, the script takes forever, around 100 seesm that work, but the fitting is not that good. Could you please help me to improve the fitting?

Thank you very much.

clear all

clc

x1 = 200:0.5:1000;

y1 = (exp(-((x1-600).^2)/(2*(100)^2)))';

x2 = (1240./(200:0.5:1000))';

center=2;

Wind=6;

NrPeaks=1;

Peak_shape=5;

extra=100; %for the exp-broad Gaussian

Numtrial=1;

start=0;

%start=[2 1];

Baseline=0;

[FitResults,GOF,baseline,coeff,residual,xi,yi,BootResults]=peakfit([x2 y1],center,Wind,NrPeaks,Peak_shape,extra,Numtrial,start,Baseline);

Andre ZeugDear Tom,

Does peakfit allow to use constraints in general, similar to input 13 (minwidth), e.g. for max witdh, and - what I am searching for - to limit the peak center position inside the data range?

Tom O'HaverThe key factors in your case are (1) divide up the problem into three groups, and (2) set baseline mode 1 (tilted linear)

#1

>> [FitResults1,FitError1]=peakfit([x y],23.204,2.7847,1,2,0,0,0,1)

Shape= Lorentzian % Fitting Error= 7.6062% R2= 0.84139

Peak# Position Height Width Area

1 23.144 650.94 0.37617 372.88

#2

>> [FitResults2,FitError2]=peakfit([x y],34.0907,7.8497,2,2,0,0,0,1)

Shape= Lorentzian % Fitting Error= 5.3057% R2= 0.96875

Peak# Position Height Width Area

1 32.711 926.45 0.55781 776.19

2 34.894 875.59 1.093 1373.4

#3

>> [FitResults3,FitError3]=peakfit([x y],47.0961,7.0224,1,2,0,0,0,1)

Shape= Lorentzian % Fitting Error= 2.4399% R2= 0.98718

Peak# Position Height Width Area

1 47.191 1793.8 0.51283 1377.3

Alternatively, if you must fit the entire signal in one go, you'll have to assign an off-screen peak to fit the baseline, then provide approximate first guesses for all 4 peaks.

The best way to explore all these options is to use the interactive peak fitter ipf.m. See https://terpconnect.umd.edu/~toh/spectrum/ifpinstructions.html

Tom O'Haver

guilhem dezanneauHello,

I would so happy to be able to use this peakfit function. But for now, i have totally failed despite spending a significant time to understand the different possibilities. I have a defined number of peaks (of typical shape pseudo-voigt = mix of Lorentzian and gaussian) and i already know quite well their position and width. So everything should be easy BUT since my background has some complex shape, i am unable to know how to manage this. My original peak always get lost since the program use one or two of them to fit the background.

please help me to understand how i can manage to fit those peaks?

Loïc TorailleTom,

Sorry for the lack of clarity. I'm talking about the variable type double: https://fr.mathworks.com/help/matlab/ref/double.html.

As for the dimensions of PEAKHEIGHT: in the fitlorentzian function, there is a line which does "If BaselineMode==3,A=[ones(size(y))' A];end". The next instruction then assigns "PEAKHEIGTS=abs(A\y');end", which means that in general PEAKHEIGHT is a vector with only one value, but when BaselineMode is set to 3, PEAKHEIGHT is a vector storing two values.

I hope this is clearer now,

Loïc

Tom O'HaverLoic Toraille, I'd like to help but I don't understand your questions. What is a "double variable"? Also, what modification has been added to the dimensions of PeakHeights in the fitlorentzian function? Sorry. I just don't understand. Tom O'Haver

Loïc TorailleHello Tom, and thank you for sharing your code!

I'm slowly getting it to work, and I've encountered two problems so far:

- The text function which prints out the results seems to need double variables for the positions, whereas in my use of the code the input was single variables. To keep it simple, I added a double conversion just in the text function and it solves the problem, but I wonder why I'm the only one encountering it.

- I had some issues when alternating between different BaselineMode, (more specifically the any mode and the 3 mode I think), and after some digging I finally got to the core of the problem: when using BaselineMode 3, there is a modification added to the dimensions of PeakHeights in the fitlorentzian function. However, that modification was not being performed when I switched the BaselineMode back and forth, even though it worked fine using your ipf.m code. I finally understood that it was because the BaselineMode variable is called as a global variable in the fitlorentzian, but it cannot work since it is an input variable of the peakfit code. To solve this, I've added BaselineMode to the list of global variables of peakfit, changed the name of the input variable to InputBaselineMode, and declared at the beginning of peakfit that BaselineMode = InputBaselineMode.

Hope this helps, let me know if I did anything wrong which resulted in these problems appearing.

Loïc

Tom O'HaverHana, the problem might be caused by the global variables AA, xxx, PEAKHEIGHTS, FIXEDPARAMETERS, delta, BIPOLAR, and CLIPHEIGHT, Make sure they don't interfere with the variables in your GUI and rename them if necessary.

Hana KalinovaThanks. Please let me know if you need any more details. It could be something I am doing wrong, but it is reproducible and I'm stuck.

Tom O'HaverHana, I've not seen that one before. I'll have a look and see if I get any ideas. Tom.

Hana KalinovaHello, I'm running into an error. The first time I call peakfit in my GUI (making sure it has a figure of its own to work with) - it fails to run the bootstrapping, and it returns x and y variables that are of different length, which causes my program to crash.

If I run a separate script that does not use a GUI, it runs the bootstrapping fine, and when I run in the GUI, it will run for the rest of the day, with no problem, until I restart my Matlab. In the separate script, the code is exactly the same for the peakfit command... so I don't see what is being initialized there that is not working properly in my GUI on the first run. Thanks for any help. For now we will run a dummy program at the beginning of the day, but would be great if I could figure it out!

Tom O'HaverAdolfo, I'm sorry I missed this. We were traveling in Mongolia at the time and I was not checking in here. Better to contact me by email (toh@umd.edu). With this much offset in your baseline, it's best to subtract out most of it beforehand. Then:

[FitResults,FitError]=peakfit([xdata ydata-4000],0,0,2,1,0,10)

FitResults =

1 650.52 1981.6 0.82016 981.27

2 650.73 5264.7 0.099534 557.8

FitError =

2.7354 0.98746

Adolfo CoboHello all and thankyou @Tom O'Haver for sharing this code. I have landed here looking for an automated way of fitting peaks in Laser-Induced Breakdown Spectroscopy (LIBS) spectra. I think that it should be straightforward, as the useful peaks I want to fit are well defined, with low noise, and isolated. The only problem is that some of them have a large baseline. But I have found that sometimes the fitting is very bad, and I don't what to try to improve. I know the expected central wavelength and width, and even with those starting values, the fit is not Ok.

One example, this is the data of the peak: https://www.dropbox.com/s/ttgx9ngdl0qqn4x/peakfitvars.mat?dl=0

And a plot: https://www.dropbox.com/s/1x3hjjnhk2uc9al/peakplot.PNG?dl=0

The fitting with these parameters (lorentzian fit):

peakfit([xdata ydata],0,0,1,2,0,99,[650.73 0.06],0,0,1)

is very poor: https://www.dropbox.com/s/990cu14o6gaezbh/fit1.PNG?dl=0

With baseline correction=1 (linear) there is no fit at all: https://www.dropbox.com/s/s52dtu8sz78ntax/fit_baseline1.PNG?dl=0

With baseline correction=2 (quadratic) there is an acceptable fitting, but the estimated height is about 3000 units, while the real one is closer to 6000.

https://www.dropbox.com/s/omd47bpo0smm0cw/fit_baseline2.PNG?dl=0

I would appreciate advice regarding which parameters to tweak in this case. Thanks in advance!

Tom O'HaverIf the only variable is the peak height, you don't need an iterative algorithm for that; you can just use the much simpler and faster classical least squares computation.

Felipe Ademir aleman hernandezHi Tom,

I am trying to figurate how can I modify the code to get a Voigt profile fit with known peak position, width and alpha and just have the height as the variable, similar to the recently added shape 50. Could you give me some guidelines? maybe be using shape 50 and modify the fitting function.

Thanks!

Tom O'HaverAlexandra, one way is to set bipolar=1 and use baseline mode 3.

For example, consider a single negative Lorentzian on a basline of 100:

x=1:100;

y=100-30*lorentzian(x,50,10);

datamatrix=[x;y];

[FitResults,FitError]=peakfit(datamatrix,0,0,1,2,1,1,0,3,0,1,1)

Alexandra HanselmanI'm trying to fit a Lorentzian to data that is offset in the positive-y direction, but which has a negative amplitude. Is there any way to specify negative amplitudes (I'm using peakshape = 2).

Tom O'HaverThanks, Andre, that looks very interesting. I will read it and see if I can figure out how to add that to my script.

Andre Sartori@Tom O'Haver, just to conclude on the harmonic oscillator topic, if you would like to understand what kind of data I was talking about, you can now check the paper I was working on when I first contacted you (link below). Figure 4 and Figure 7 are the sort of peaks I was trying to fit. Fitting that data to a Lorentzian profile was not physically correct, so I asked you about the idea of adding also the harmonic oscillator in your script. https://onlinelibrary.wiley.com/doi/full/10.1002/smll.201803774

Cheers,

A.

James DrewittTom O'HaverTeodoro,

1. You can cite my website at https://terpconnect.umd.edu/~toh/spectrum/ or my book "Pragmatic Introduction to Signal Processing 2019: Applications in scientific measurement", on Amazon at https://www.amazon.com/dp/1792867123

2. I don't see this behavior. What version of peakfit.m are you using? Can you send me a screen image or an example segment of data that demonstrates the problem? email: toh@umd.edu.

Teodoro KaufmanDear Tom, This is a great contribution. Have two questions:

1. How to correctly cite this in a paper.

2. The version I have works well for fitting two Lorentzians to overlapped NMR signals; however, when I try to fit a third Lorentzian, it gives me wrong parameters for the third curve (on the Matlab workspace; surprisingly not on the Figure). Have you heared about this detail?

Finally, have a wonderful 2019!

Tom O'HaverYes, I just posted the new version on that page.

Joao FigueiraThank you Tom. I realize I should probably have commented on the ipf page.

Calvin O'HaverJoao,

That's a good idea. I'll add that to a new version and upload it.

Tom O'Haver

Joao FigueiraI'm using the interactive function. Is there a more straightforward way to save the model data table, that to print to the terminal (keyboard d)?

Andre Sartori@Tom O'Haver, following up on the harmonic oscillator, I don't know exactly what combination of parameters you're trying, but I think none of them can be <=0. That also includes all values of x, because you cannot have zero frequency or negative frequency. The condition that c is approximately equal to the FWHM only holds when c << b. That PDF form Harvard on the harmonic oscillator helps you understand that the equation is valid for a special case of damped/driven harmonic oscillator, and contains a few indications of the boundary conditions, such as the one I just mentioned. Meanwhile, I wrote my own fitting code for the harmonic oscillator, using the fitnlm function from Matlab. I could get the standard error from there as well, so my problem is solved. But I still think that this equation would be a very useful addition to your script. Cheers.

Tom O'HaverOops, there was a typo in that first section fit. It should read:

>> [FitResults,FitError]=peakfit([x10;y],-0.25,3.5,3,1,0,10,[-0.1 0.5 1 0.5 2 1])

Tom O'HaverAttempting to fit the entire signal in one fitting operation fails, because the least-squares method tries to minimize the overall residual error for the whole signal. Since the left-hand two peaks are so much weaker than the main peak on the right, the algorithm ignores the two smaller peaks because their contribution to the total error is insignificant compared to the two larger peaks. In other words, it can get a lower total error by adjusting the fit to the larger peaks only.

The solution to this is to fit the signal in two separate sections. I used the interactive peak fitter ipf.m to zoom in to different sections and try various fits before settling on these two (using the W key to print out the peakfit syntax. Note that in the first fit, a third model peak is added just to account for the background cased bt the off-screen peak to its right.

FIRST SECTION

>> [FitResults,FitError]=peakfit92([x10;y],-0.25,3.5,3,1,0,10,[-0.1 0.5 1 0.5 2 1])

FitResults =

Peak# Position Height Width Area

1 -0.098159 148.9 0.48739 77.252

2 0.99263 144.2 0.49226 75.028

3 1.8245 647.92 0.90581 127.27

FitError =

0.98664 0.99824

-----------------------------------------------------------------------------------------------

SECOND SECTION

>> [FitResults,FitError]=peakfit([x10;y],2.4,2.1,2,1)

FitResults =

Peak# Position Height Width Area

1 2.5644 3619.8 1.1759 4310.1

2 2.9925 2366.3 0.77348 1775.2

FitError =

1.3845 0.99814

Bruno HonórioDear Tom... I´m uploading both figure and the data. i would like to "impose" 4 peaks (two at the left hand side, and two inside the main peak). I´ve tried many parameters without sucess. If you could help me, would be great.

Data:

https://www.dropbox.com/s/zqkd1kofv3ckpb5/2share.mat?dl=0

Reference figure: https://www.dropbox.com/s/e3qtup1blt0xog6/2share.png?dl=0

Failed attempt: https://www.dropbox.com/s/2khpy1tdl1qh4pw/2share_failed.png?dl=0

Thanks again..

Best wishes

Tom O'HaverBruno Honório,

It does look like your data should be well fit by a 3-Gaussian model using peakfit.m. I can not demonstrate it because my function fits numerical data, not graphical data.

Bruno HonórioDear Tom,

Thanks for sharing the codes.

I´m working on a NMR experiment and the answer can be understood as a Gaussian mixture. So, I´m trying to use your code to decompose the signal. The answer I´m getting can be seem here: https://www.dropbox.com/s/2qsemu0rdzcesu0/Capture.PNG?dl=0

I´m using a "trying and error" approach to fit the components and I would like to know if there is a way to fine tunning the peaks and interval. Note some peak shifts (mostly in the comp.2) and the tail on the right hand side.

Thanks in advance

Bruno

Tom O'HaverAndre Sartori, for me, that expression, y=a/(((b*2*pi)^2-(x*2*pi)^2)^2+(c*x*4*pi*pi)^2)^0.5, or in Matlab syntax, y=a./(((b.*2.*pi).^2-(x.*2.*pi).^2).^2+(c.*x*4.*pi.*pi).^2).^0.5;, results in a variety of shapes depending on the values of b, and c, including double peak shapes that are symmetrical on x=0. Moreover, c is the FWHM only if C=1. Is that what I should expect? I'm just not sure I have the expression right.

ThorBarraAndre SartoriAndre Sartori@Tom O'Haver as an example, I fitted the amplitude vs. frequency of a resonant damped/driven harmonic oscillator manually:

https://www.dropbox.com/s/r9hhp7dmzf73i0u/harmonic_oscillator.png

This is the equation that I find useful, i.e. with the frequency in Hz instead of radians, because typically the measurements are done that way, and also presented that way.

y=a/(((b*2*pi)^2-(x*2*pi)^2)^2+(c*x*4*pi*pi)^2)^0.5

Then the output becomes directly comparable to the output of the Gaussian or Lorentzian profiles, with c being the FWHM (in Hz) in case the oscillator is under-damped (which is up to the user to know), b being the peak center (in Hz) and a equal to F0/m which is not quite a height, but which can be treated as such.

Then the calculation of the important Q-factor is simply b/c. If I get the std of these parameters, with simple error propagation I have the uncertainty of Q, which is how I currently do with the Lorentzian.

Cheers,

Andre Sartori@Tom O'Haver, you can check the equation for the amplitude, A(omega) here: http://ipl.physics.harvard.edu/wp-uploads/2013/03/15c_s07_4.pdf

I guess you can probably treat F0/m as constant a, omega_zero as constant b, and gamma as constant c.

But I don't think there will be a problem with infinite solutions due to 2*pi radians, because there's no sin or cos function in the equation. It's not periodic. The shape of the equation is similar to a Lorentzian profile, centered at omega_zero.

Right now I'm working on some resonance measurements of MEMS devices, and I've been applying the Lorentzian profile to fit the peaks and obtain the so-called Q-factor, which is the peak center divided by the FWHM, but this is not the correct curve from a physics point-of-view, because the devices resonate according to the damped/driven harmonic oscillator. The shapes don't quite match and the Q-factors are not quite right. So I need to re-fit everything to the harmonic oscillator equation. But my problem is that I can't easily fit the equation and at the same time get the standard deviation, which is a very handy feature of your peakfit function. I want to have a number for the Q-factor and a plus/minus something.

Cheers.

Tom O'HaverGood idea, Andre Sartori. I'll give it a try. Presumably there would be three nonlinear variables, the frequency and phase of the oscillation and the exponential decay constant, right? I anticipate that there might be a problem with the phase variable because there should be an infinite number of solutions with phases spaced 2*pi radians (360 degrees) apart.

Andre SartoriDear Tom, thanks for the script. Very handy. Any chance of adding also a damped/driven harmonic oscillator function in the futuer? Cheers.

Tom O'HaverAnwesha Palit, yes:

Calculating the precision of the peak parameters:

[FitResults, GOF, baseline, coeff, residuals, xi, yi, BootstrapErrors] = peakfit([x;y],0,0,2,6,0,1,0,0,0);

Displays parameter error estimates by the bootstrap method; See DemoPeakfitBootstrap for a self-contained demo of this function.

Also see Example 15 on https://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm

Anwesha PalitCan peakfit.m provide me with an error analysis for every parameter it measures ?

Tim SieberMay you can add lorentian and gaussian models where I can set borders for the position and width?

Tim SieberHello Tom,

I have a problem when I try to fit a fixed position gaussian peak and 2 lorentzian peak (variable) together.

My call is this: [FitResults,FitError]=peakfit(M1, ... % datamatrix

521, ... % center

100, ... % window

4, ... % Number of peaks

[16 2 1 ], ... % peakshape

[1 0 0], ... % extra

[10 1 1], ... % numtrials

[521 3.0 505 3.0 517 3.0 ], ... % first guess

0, ... % autozero

[521] ... % fixedparameters

)

The problem is, the 521 isn't fixed at all. May you have a hint for me.

Also whe I have the Peakshape in another order like [1 16 2] it will not fit at all. It just give me this back:

Index exceeds matrix dimensions.

Error in peakfit>fitmultiple (line 3539)

A(j,:)=gaussian(xx,FIXEDPARAMETERS(j),lambda(j));

Error in peakfit>@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)) (line 985)

TrialParameters=fminsearch(@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)),newstart,options);

Error in fminsearch (line 189)

fv(:,1) = funfcn(x,varargin{:});

Error in peakfit (line 985)

TrialParameters=fminsearch(@(lambda)(fitmultiple(lambda,xx,yy,NumPeaks,shapesvector,extra)),newstart,options);

Error in excel1 (line 53)

[FitResults,FitError]=peakfit(M1, ... % datamatrix

>>

May you can help me, cause I need it to fit a huge set of data.

Tom O'Haveryanhua yao, thanks for your rating! Best of luck in your work. Tom

yanhua yaoTom O'Haverdesert fish, the 4th column of the FitResults matrix returned by the peakfit function is the full width at half maximum. The left and right bounds depend on how you define them: e.g. at 1% of peak height, 0.1% or peak height, etc. You would have to calculate that yourself depending on your choice. See https://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm#PeakStartAndEnd

desert fishhow can i find the peaks' left bound and right bound?Also the halfheight width?

thank you.

Tian TianHi Tom, thanks for sharing this great tool. Is peakfit a defined function? I got the error: Undefined function or variable 'peakfit'. Could you please tell me how to define it？ Thanks.

Tom O'HaverSebastien, just include all the output arguments up to and includinf xi and yi. Here's an example:

x=7:.05:13;

y=x.^3/500+exp(-(x-9).^2)+exp(-(x-11).^2)+.0.*randn(size(x));

[FitResults,GOF,baseline,coeff,residual,xi,yi]=peakfit([x;y],0,0,3,[1 1 46],[1 1 1],10);

plot(xi,yi)

Tom O'HaverIf anyone else has a question, you will get faster response from me if you sent it to my email toh@umd.edu.

saiedeh kkplease send me this file per email , I could not see it

SebI finally manage.

In order to retrieve data (array) values from the peakfit plot, follow these steps:

1) Run the script,

2) Type the following commands:

h = findobj(gca,'Type','line'); x=get(h,'Xdata');

y=get(h,'Ydata');

3) Since we deal with several curves, to plot a unique curve you have to use indexing. For instance for the first green curve:

plot(x{1},y{1})

Hope this helps !

SebDear Tom,

First, thanks for this amazing tool.

I'm not smart enough to get rid of a (simple?) problem.

Reading the comments, I can plot the green curves. Right.

However, I have no clue about getting the (xi) & (yi) values (same dimensions) to continue the game ! Indeed, the matlab plot function allows plotting a curve from a scalar (xi) and a matrix (yi). Is there any way to get the related xi & yi points ?!

This might be really useful for coming users...

Thanks for your consideration.

SebSourinTom O'HaverAmir,

I don't know about parallel processing. Does this require special hardware?

Tom O'HaverJohan Manuel, do you mean computing the Fouruer transform NOT using Matlab's built-in FFT function? How about this? http://www.mathworks.com/matlabcentral/fileexchange/2271-numerical-methods-for-physics/content/edition1/matlab4/sft.m

amir afrashtehpourby the way have you thought about using parallel computing in anyways to speed up the code?

thanks

Death SaurerTom O'Haver

That so very good day.

Let me turn to you asking for your kind help in order to get the code in matlab for the fourier transform of continuous-time (FFT) directly and inversely, since my knowledge in matlab are not very good.

I appreciate who I can help and collaborate.

Excuse my bad English and my native language is Spanish

LongTom O'HaverKnut,

Yes. Just call peakfit with the output arguments [FitResults,GOF,baseline,coeff,residual,xi,yi],

then plot(xi,yi) will plot the model components separately (the "green" line) and plot(x,y-residual) will plot the combined model (the "red" line). You can make then whatever color and style you want.

Knut TutgutDear Tom,

is there way to simply extract the curves of the models?

So in principle, how can I plot the green and red line in a separate graph?

Thanks a lot! Great tool!

Tom O'HaverSanjay, the new version 7.4 includes fixed-height Gaussians (Shape 34). See Example 30 in the help file.

Sanjay SekaranHi Tom,

Would it be possible to specify fixed height gaussians during the fitting? Thanks.

Tom O'HaverNikola, 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

Tom O'HaverVictoria,

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

VictoriaThis 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);

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

Christian StoneWhat lines must I alter? Sorry I am pretty bad

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

CKanellaswow!!

Tom O'HaverAny 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.

Christian StoneI 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

xubikerTom O'HaverNazmul, 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.

Nazmul ArefinTom, 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.

Tom O'HaverThanks, 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.

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

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

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

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

Tom O'HaverOn 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.

Tom O'HaverExcellent idea, Hao Wang; I'll work on that.

Hao WangIf 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!

Tom O'HaverAli, 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?

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

Jade HeThank 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.

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

Jade HeHello Tom,

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

Jade

Jade HeThank you

BenThanks Tom!

Tom O'HaverHere'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.

Ben@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!

TobiasHi Tom,

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

Tom O'HaverTry using a negative value for the Fano factor (the input argument called 'extra')

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

Tom O'HaverEver 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)

TobiasHi 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

Tom O'HaverOK, 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)

TobiasHi 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

Tom O'HaverExcellent 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).

TobiasHi 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

FengHi Tom,

Thanks a lot for sharing this great work!

Tom O'HaverAs 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).

Tom O'HaverPrompted 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.

Tom O'HaverI 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.

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?

Tom O'HaverMatthias,

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.

MatthiasHi 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

Fan D.ChenThis program is perfect, I have to say, thank you

SpencerA couple of errors in the post below:

sqrt(2 log(2)) = 1.177

sqrt(log(2)) = 0.832

SpencerHi 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

Tom O'HaverThanks, 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.

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

BenHi 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'HaverBen, 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?

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

Tom O'HaverVersion 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'HaverSo 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.

RemsThanks 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'HaverThanks, 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).

RemsHi,

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'HaverVersion 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'HaverThat's a good suggestion, Himmst. I'll add an graphing ON/OFF switch to the next version. Thanks.

HimmstFirst 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 NZNThese are very good script. Thank you very much. The scripts are very useful to help to solve my problem, and my headache

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

HusseinThank you for the valuable information!

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

HusseinI 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'HaverIt 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.

Husseinbut 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'HaverNegative 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.

HusseinThank 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'HaverVersion 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 DzyubakUsing 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'HaverHussein, for that purpose you can use the much simpler technique of multilinear regression. For example, see http://terpconnect.umd.edu/~toh/spectrum/CurveFittingB.html

Husseinvery 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'HaverVersion 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'HaverYes, 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'HaverBogdan, 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 DzyubakHi, 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'HaverThat 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.

KrestenGreat program!

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

Tom O'HaverThanks 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 DzyubakDoes 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.

AMI 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'HaverAM, 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'HaverAM, 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.

AMWith 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 Hidayatnice program. thanks for sharing

SadikGood job!

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

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

thanks

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

Tom O'HaverThe 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'HaverSorry 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 B.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'HaverAs 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'HaverNo, no toolboxes required.

Victor ZhouIs it require any toolbox to run the program?

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

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

MB46This 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 KWonderful 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'HaverTry turning your 'signal" matrix into its transpose, signal'

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