| [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedparameters,plots)
|
function [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedparameters,plots)
% Version 3.6: February, 2013. Addition of fixed-position
% Gaussian shape (16) and fixed-position Lorentzian shape (17).
% A command-line peak fitting program for time-series signals,
% written as a self-contained Matlab function in a single m-file.
% Uses an non-linear optimization algorithm to decompose a complex,
% overlapping-peak signal into its component parts. The objective
% is to determine wpeakhether your signal can be represented as the sum of
% fundamental underlying peaks shapes. Accepts signals of any length,
% including those with non-integer and non-uniform x-values. Fits any
% number of peaks of Gaussian, Lorentzian, equal-width Gaussian and
% Lorentzian, fixed-width Gaussian and Lorentzian, biburfated Gaussian
% and Lorentzian, exponentially-broadened Gaussian, Pearson, Logistic,
% exponential pulse, sigmoid, and Gaussian/Lorentzian blend. 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. (c) Tom O'Haver
%
% peakfit(signal);
% Performs an iterative least-squares fit of a single Gaussian
% peak to the data matrix "signal", which has x values
% in column 1 and Y values in column 2 (e.g. [x y])
%
% peakfit(signal,center,window);
% Fits a single Gaussian peak to a portion of the
% matrix "signal". The portion is centered on the
% x-value "center" and has width "window" (in x units).
%
% peakfit(signal,center,window,NumPeaks);
% "NumPeaks" = number of peaks in the model (default is 1 if not
% specified). No limit to maximum number of peaks in version 3.1
%
% peakfit(signal,center,window,NumPeaks,peakshape);
% Specifies the peak shape of the model: "peakshape" = 1-15. (1=Gaussian
% (default), 2=Lorentzian, 3=logistic, 4=Pearson, 5=exponentionally
% broadened Gaussian; 6=equal-width Gaussians; 7=Equal-width Lorentzians;
% 8=exponentionally broadened equal-width Gaussian, 9=exponential pulse,
% 10=sigmoid, 11=Fixed-width Gaussian, 12=Fixed-width Lorentzian;
% 13=Gaussian/Lorentzian blend; 14=BiGaussian, 15=BiLorentzian,
% 16=Fixed-position Gaussians; 17=Fixed-position Lorentzians;
%
% peakfit(signal,center,window,NumPeaks,peakshape,extra) Specifies the
% value of 'extra', used in the Pearson, exponentionally broadened
% Gaussian, Gaussina/:orentzian blend, and the bifurcated Gaussian and
% Lorentzian shapes to fine-tune the peak shape.
%
% peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials);
% Performs "NumTrials" trial fits and selects the best one (with lowest
% fitting error). NumTrials can be any positive integer (default is 1).
%
% peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start)
% Specifies the first guesses vector "firstguess" for the peak positions
% and widths, e.g. start=[position1 width1 position2 width2 ...]
%
% [FitResults,MeanFitError]=peakfit(signal,center,window...)
% Returns the FitResults vector in the order peak number, peak
% position, peak height, peak width, and peak area), and the MeanFitError
% (the percent RMS difference between the data and the model in the
% selected segment of that data) of the best fit.
%
% [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit(signal,...)
% Prints out parameter error estimates for each peak fit.
%
% Optional output parameters
% 1. FitResults: a table of model peak parameters, one row for each peak,
% listing Peak number, Peak position, Height, Width, and Peak area.
% 2. LowestError: The rms fitting error of the best trial fit.
% 3. BestStart: the starting guesses that gave the best fit.
% 4. xi: vector containing 100 interploated x-values for the model peaks.
% 5. yi: matrix containing the y values of each model peak at each xi.
% Type plot(xi,yi(1,:)) to plot peak 1 or plot(xi,yi) to plot all peaks
% 6. BootResults: a table of bootstrap precision results for a each peak
% and peak parameter.
%
% Example 1:
% >> x=[0:.1:10]';y=exp(-(x-5).^2);peakfit([x y])
% Fits exp(-x)^2 with a single Gaussian peak model.
%
% Peak number Peak position Height Width Peak area
% 1 5 1 1.665 1.7725
%
% >> y=[0 1 2 4 6 7 6 4 2 1 0 ];x=1:length(y);
% >> peakfit([x;y],length(y)/2,length(y),0,0,0,0,0,0)
% Fits small set of manually entered y data to a single Gaussian peak model.
%
% Example 2:
% x=[0:.1:10]';y=exp(-(x-5).^2)+.1*randn(size(x));peakfit([x y])
% Like Example 1, except that random noise is added to the y data.
% ans =
% 1 5.0279 0.9272 1.7948 1.7716
%
% Example 3:
% x=[0:.1:10];y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.1*randn(size(x));
% peakfit([x' y'],0,0,2)
% Fits a noisy two-peak signal with a double Gaussian model (NumPeaks=2).
% ans =
% 1 3.0001 0.49489 1.642 0.86504
% 2 4.9927 1.0016 1.6597 1.7696
%
% Example 4:
% >> x=1:100;y=ones(size(x))./(1+(x-50).^2);peakfit(y,0,0,1,2)
% Fit Lorentzian (peakshape=2) located at x=50, height=1, width=2.
% ans =
% 1 50 0.99974 1.9971 3.1079
%
% Example 5:
% >> x=[0:.005:1];y=humps(x);peakfit([x' y'],.3,.7,1,4,3);
% Fits a portion of the humps function, 0.7 units wide and centered on
% x=0.3, with a single (NumPeaks=1) Pearson function (peakshape=4) with
% extra=3 (controls shape of Pearson function).
%
% Example 6:
% >> x=[0:.005:1];y=(humps(x)+humps(x-.13)).^3;smatrix=[x' y'];
% >> [FitResults,MeanFitError]=peakfit(smatrix,.4,.7,2,1,0,10)
% Creates a data matrix 'smatrix', fits a portion to a two-peak Gaussian
% model, takes the best of 10 trials. Returns FitResults and MeanFitError.
% FitResults =
% 1 0.31056 2.0125e+006 0.11057 2.3689e+005
% 2 0.41529 2.2403e+006 0.12033 2.8696e+005
% MeanFitError =
% 1.1899
%
% Example 7:
% >> peakfit([x' y'],.4,.7,2,1,0,10,[.3 .1 .5 .1]);
% As above, but specifies the first-guess position and width of the two
% peaks, in the order [position1 width1 position2 width2]
%
% Example 8:
% >>x=[0:.1:10]';y=10-x+exp(-(x-5).^2);peakfit([x y],5,8,0,0,0,0,0,1)
% Fitting single Gaussian on a linear background, using linear autozero (1)
% AUROZERO=0 does not subtract baseline from data segment.
% AUROZERO=1 (default) subtracts linear baseline from data segment.
% AUROZERO=2, subtracts quadratic baseline from data segment.
%
% >> x=[0:.1:10]';y=1./(1+x.^2)+exp(-(x-5).^2);peakfit([x y],5,5.5,0,0,0,0,0,2)
% ans =
% 1 5.0078 0.97001 1.6079 1.6598
% Fitting single Gaussian on a curved background, using quadratic autozero
% (2) and specifying center and window, but using placeholders (zeros) to
% use default values for NumPeaks, peakshape, extra, NumTrials, and start.
%
% Example 9:
% >> x=[0:.1:10];y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.1*randn(size(x));
% [FitResults,MeanFitError]=peakfit([x' y'],0,0,2,11,0,0,0,0,1.666)
% Same as example 3, fit with fixed-width Gaussian (shape 11), width=1.666
%
% Example 10: (Version 3 or later; Prints out parameter error estimates)
% >> x=0:.05:9;y=exp(-(x-5).^2)+.5*exp(-(x-3).^2)+.01*randn(1,length(x));
% >> [FitResults,LowestError,BestStart,xi,yi,BootstrapErrors]=peakfit([x;y],0,0,2,6,0,1,0,0,0);
%
% Example 11: (Version 3.2 or later)
% >> x=[0:.005:1];y=humps(x);[FitResults,MeanFitError]=peakfit([x' y'],0.54,0.93,2,13,15,10,0,0,0)
% FitResults =
% 1 0.30078 190.41 0.19131 23.064
% 2 0.89788 39.552 0.33448 6.1999
% MeanFitError =
% 0.34502
% Fits both peaks of the Humps function with a Gaussian/Lorentzian blend
% (shape 13) that is 15% Gaussian (Extra=15).
%
% Example 12: (Version 3.2 or later)
% >> x=[0:.1:10];y=exp(-(x-4).^2)+.5*exp(-(x-5).^2)+.01*randn(size(x));
% >> [FitResults,MeanFitError]=peakfit([x' y'],0,0,1,14,45,10,0,0,0)
% FitResults =
% 1 4.2028 1.2315 4.077 2.6723
% MeanFitError =
% 0.84461
% Fit a slightly asymmetrical peak with a bifurcated Gaussian (shape 14)
%
% Example 13: (Version 3.3 or later)
% >> x=[0:.1:10]';y=exp(-(x-5).^2);peakfit([x y],0,0,1,1,0,0,0,0,0,0)
% Example 1 without plotting (11th input argument = 0, default is 1)
%
% For more details, see
% http://terpconnect.umd.edu/~toh/spectrum/CurveFittingC.html and
% http://terpconnect.umd.edu/~toh/spectrum/InteractivePeakFitter.htm
%
global AA xxx PEAKHEIGHTS FIXEDWIDTH FIXEDPOSITIONS
format short g
format compact
warning off all
NumArgOut=nargout;
datasize=size(signal);
if datasize(1)<datasize(2),signal=signal';end
datasize=size(signal);
if datasize(2)==1, % Must be isignal(Y-vector)
X=1:length(signal); % Create an independent variable vector
Y=signal;
else
% Must be isignal(DataMatrix)
X=signal(:,1); % Split matrix argument
Y=signal(:,2);
end
X=reshape(X,1,length(X)); % Adjust X and Y vector shape to 1 x n (rather than n x 1)
Y=reshape(Y,1,length(Y));
% If necessary, flip the data vectors so that X increases
if X(1)>X(length(X)),
disp('X-axis flipped.')
X=fliplr(X);
Y=fliplr(Y);
end
% Isolate desired segment from data set for curve fitting
if nargin==1 || nargin==2,center=(max(X)-min(X))/2;window=max(X)-min(X);end
xoffset=center-window/2;
xoffset=0;
n1=val2ind(X,center-window/2);
n2=val2ind(X,center+window/2);
if window==0,n1=1;n2=length(X);end
xx=X(n1:n2)-xoffset;
yy=Y(n1:n2);
ShapeString='Gaussian';
% Define values of any missing arguments
switch nargin
case 1
NumPeaks=1;
peakshape=1;
extra=0;
NumTrials=1;
xx=X;yy=Y;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
plots=1;
case 2
NumPeaks=1;
peakshape=1;
extra=0;
NumTrials=1;
xx=signal;yy=center;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
plots=1;
case 3
NumPeaks=1;
peakshape=1;
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 4
peakshape=1;
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 5
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 6
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 7
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 8
AUTOZERO=1;
FIXEDWIDTH=0;
plots=1;
case 9
FIXEDWIDTH=0;
plots=1;
case 10
FIXEDWIDTH=fixedparameters;
plots=1;
case 11
FIXEDWIDTH=fixedparameters;
otherwise
end % switch nargin
% Default values for placeholder zeros
if NumTrials==0;NumTrials=1;end
if peakshape==0;peakshape=1;end
if NumPeaks==0;NumPeaks=1;end
if start==0;start=calcstart(xx,NumPeaks,xoffset);end
if FIXEDWIDTH==0, FIXEDWIDTH=length(xx)/10;end
if peakshape==16;FIXEDPOSITIONS=fixedparameters;end
% Remove linear baseline from data segment if AUTOZERO==1
bkgsize=round(length(xx)/10);
if bkgsize<2,bkgsize=2;end
lxx=length(xx);
if AUTOZERO==1, % linear autozero operation
XX1=xx(1:round(lxx/bkgsize));
XX2=xx((lxx-round(lxx/bkgsize)):lxx);
Y1=yy(1:(round(length(xx)/bkgsize)));
Y2=yy((lxx-round(lxx/bkgsize)):lxx);
bkgcoef=polyfit([XX1,XX2],[Y1,Y2],1); % Fit straight line to sub-group of points
bkg=polyval(bkgcoef,xx);
yy=yy-bkg;
end % if
if AUTOZERO==2, % Quadratic autozero operation
XX1=xx(1:round(lxx/bkgsize));
XX2=xx((lxx-round(lxx/bkgsize)):lxx);
Y1=yy(1:round(length(xx)/bkgsize));
Y2=yy((lxx-round(lxx/bkgsize)):lxx);
bkgcoef=polyfit([XX1,XX2],[Y1,Y2],2); % Fit parabola to sub-group of points
bkg=polyval(bkgcoef,xx);
yy=yy-bkg;
end % if autozero
PEAKHEIGHTS=zeros(1,NumPeaks);
n=length(xx);
newstart=start;
% for peaks=1:NumPeaks,
% peakindex=2*peaks-1;
% newstart(peakindex)=start(peakindex)-xoffset;
% end
% Assign ShapStrings
switch peakshape
case 1
ShapeString='Gaussian';
case 2
ShapeString='Lorentzian';
case 3
ShapeString='Logistic';
case 4
ShapeString='Pearson';
case 5
ShapeString='ExpGaussian';
case 6
ShapeString='Equal width Gaussians';
case 7
ShapeString='Equal width Lorentzians';
case 8
ShapeString='Exp. equal width Gaussians';
case 9
ShapeString='Exponential Pulse';
case 10
ShapeString='Sigmoid';
case 11
ShapeString='Fixed-width Gaussian';
case 12
ShapeString='Fixed-width Lorentzian';
case 13
ShapeString='Gaussian/Lorentzian blend';
case 14
ShapeString='BiGaussian';
case 15
ShapeString='BiLorentzian';
case 16
ShapeString='Fixed-position Gaussians';
case 17
ShapeString='Fixed-position Lorentzians';
otherwise
end % switch peakshape
% Perform peak fitting for selected peak shape using fminsearch function
options = optimset('TolX',.001,'Display','off' );
LowestError=1000; % or any big number greater than largest error expected
FitParameters=zeros(1,NumPeaks.*2);
BestStart=zeros(1,NumPeaks.*2);
height=zeros(1,NumPeaks);
bestmodel=zeros(size(yy));
for k=1:NumTrials,
% disp(['Trial number ' num2str(k) ] ) % optionally prints the current trial number as progress indicator
switch peakshape
case 1
TrialParameters=fminsearch(@(lambda)(fitgaussian(lambda,xx,yy)),newstart,options);
case 2
TrialParameters=fminsearch(@(lambda)(fitlorentzian(lambda,xx,yy)),newstart,options);
case 3
TrialParameters=fminsearch(@(lambda)(fitlogistic(lambda,xx,yy)),newstart,options);
case 4
TrialParameters=fminsearch(@(lambda)(fitpearson(lambda,xx,yy,extra)),newstart,options);
case 5
zxx=[zeros(size(xx)) xx zeros(size(xx)) ];
zyy=[zeros(size(yy)) yy zeros(size(yy)) ];
TrialParameters=fminsearch(@(lambda)(fitexpgaussian(lambda,zxx,zyy,-extra)),newstart,options);
case 6
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitewgaussian(lambda,xx,yy)),cwnewstart,options);
case 7
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitlorentziancw(lambda,xx,yy)),cwnewstart,options);
case 8
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitexpewgaussian(lambda,xx,yy,-extra)),cwnewstart,options);
case 9
TrialParameters=fminsearch(@(lambda)(fitexppulse(lambda,xx,yy)),newstart,options);
case 10
TrialParameters=fminsearch(@(lambda)(fitsigmoid(lambda,xx,yy)),newstart,options);
case 11
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFWGaussian(lambda,xx,yy)),fixedstart,options);
case 12
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFWLorentzian(lambda,xx,yy)),fixedstart,options);
case 13
TrialParameters=fminsearch(@(lambda)(fitGL(lambda,xx,yy,extra)),newstart,options);
case 14
TrialParameters=fminsearch(@(lambda)(fitBiGaussian(lambda,xx,yy,extra)),newstart,options);
case 15
TrialParameters=fminsearch(@(lambda)(fitBiLorentzian(lambda,xx,yy,extra)),newstart,options);
case 16
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=(max(xx)-min(xx))./(NumPeaks+1);
fixedstart(pc)=fixedstart(pc)+.1*randn().*fixedstart(pc);
end
TrialParameters=fminsearch(@(lambda)(FitFPGaussian(lambda,xx,yy)),fixedstart,options);
case 17
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=(max(xx)-min(xx))./(NumPeaks+1);
fixedstart(pc)=fixedstart(pc)+.1*randn().*fixedstart(pc);
end
TrialParameters=fminsearch(@(lambda)(FitFPLorentzian(lambda,xx,yy)),fixedstart,options);
otherwise
end % switch peakshape
% Construct model from Trial parameters
A=zeros(NumPeaks,n);
for m=1:NumPeaks,
switch peakshape
case 1
A(m,:)=gaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 2
A(m,:)=lorentzian(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 3
A(m,:)=logistic(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 4
A(m,:)=pearson(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 5
A(m,:)=expgaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m),-extra)';
case 6
A(m,:)=gaussian(xx,TrialParameters(m),TrialParameters(NumPeaks+1));
case 7
A(m,:)=lorentzian(xx,TrialParameters(m),TrialParameters(NumPeaks+1));
case 8
A(m,:)=expgaussian(xx,TrialParameters(m),TrialParameters(NumPeaks+1),-extra)';
case 9
A(m,:)=exppulse(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 10
A(m,:)=sigmoid(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 11
A(m,:)=gaussian(xx,TrialParameters(m),FIXEDWIDTH);
case 12
A(m,:)=lorentzian(xx,TrialParameters(m),FIXEDWIDTH);
case 13
A(m,:)=GL(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 14
A(m,:)=BiGaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 15
A(m,:)=BiLorentzian(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 16
A(m,:)=gaussian(xx,FIXEDPOSITIONS(m),TrialParameters(m));
case 17
A(m,:)=lorentzian(xx,FIXEDPOSITIONS(m),TrialParameters(m));
otherwise
end % switch
% for parameter=1:2:2*NumPeaks,
% newstart(parameter)=newstart(parameter)*(1+randn/50);
% newstart(parameter+1)=newstart(parameter+1)*(1+randn/10);
% end
end % for
% Multiplies each row by the corresponding amplitude and adds them up
model=PEAKHEIGHTS'*A;
% Compare trial model to data segment and compute the fit error
MeanFitError=100*norm(yy-model)./(sqrt(n)*max(yy));
% Take only the single fit that has the lowest MeanFitError
if MeanFitError<LowestError,
if min(PEAKHEIGHTS)>0, % Consider only fits with positive peak heights
LowestError=MeanFitError; % Assign LowestError to the lowest MeanFitError
FitParameters=TrialParameters; % Assign FitParameters to the fit with the lowest MeanFitError
BestStart=newstart; % Assign BestStart to the start with the lowest MeanFitError
height=PEAKHEIGHTS; % Assign height to the PEAKHEIGHTS with the lowest MeanFitError
bestmodel=model; % Assign bestmodel to the model with the lowest MeanFitError
end % if min(PEAKHEIGHTS)>0
end % if MeanFitError<LowestError
end % for k (NumTrials)
%
% Construct model from best-fit parameters
AA=zeros(NumPeaks,200);
xxx=linspace(min(xx),max(xx),200);
for m=1:NumPeaks,
switch peakshape
case 1
AA(m,:)=gaussian(xxx,FitParameters(2*m-1),FitParameters(2*m));
case 2
AA(m,:)=lorentzian(xxx,FitParameters(2*m-1),FitParameters(2*m));
case 3
AA(m,:)=logistic(xxx,FitParameters(2*m-1),FitParameters(2*m));
case 4
AA(m,:)=pearson(xxx,FitParameters(2*m-1),FitParameters(2*m),extra);
case 5
AA(m,:)=expgaussian(xxx,FitParameters(2*m-1),FitParameters(2*m),-extra*length(xxx)./length(xx))';
case 6
AA(m,:)=gaussian(xxx,FitParameters(m),FitParameters(NumPeaks+1));
case 7
AA(m,:)=lorentzian(xxx,FitParameters(m),FitParameters(NumPeaks+1));
case 8
AA(m,:)=expgaussian(xxx,FitParameters(m),FitParameters(NumPeaks+1),-extra*length(xxx)./length(xx))';
case 9
AA(m,:)=exppulse(xxx,FitParameters(2*m-1),FitParameters(2*m));
case 10
AA(m,:)=sigmoid(xxx,FitParameters(2*m-1),FitParameters(2*m));
case 11
AA(m,:)=gaussian(xxx,FitParameters(m),FIXEDWIDTH);
case 12
AA(m,:)=lorentzian(xxx,FitParameters(m),FIXEDWIDTH);
case 13
AA(m,:)=GL(xxx,FitParameters(2*m-1),FitParameters(2*m),extra);
case 14
AA(m,:)=BiGaussian(xxx,FitParameters(2*m-1),FitParameters(2*m),extra);
case 15
AA(m,:)=BiLorentzian(xxx,FitParameters(2*m-1),FitParameters(2*m),extra);
case 16
AA(m,:)=gaussian(xxx,FIXEDPOSITIONS(m),FitParameters(m));
case 17
AA(m,:)=lorentzian(xxx,FIXEDPOSITIONS(m),FitParameters(m));
otherwise
end % switch
end % for
% Multiplies each row by the corresponding amplitude and adds them up
heightsize=size(height');
AAsize=size(AA);
if heightsize(2)==AAsize(1),
mmodel=height'*AA;
else
mmodel=height*AA;
end
% Top half of the figure shows original signal and the fitted model.
if plots,
subplot(2,1,1);plot(xx+xoffset,yy,'b.'); % Plot the original signal in blue dots
hold on
end
for m=1:NumPeaks,
if plots, plot(xxx+xoffset,height(m)*AA(m,:),'g'),end % Plot the individual component peaks in green lines
area(m)=trapz(xxx+xoffset,height(m)*AA(m,:)); % Compute the area of each component peak using trapezoidal method
yi(m,:)=height(m)*AA(m,:); % (NEW) Place y values of individual model peaks into matrix yi
end
xi=xxx+xoffset; % (NEW) Place the x-values of the individual model peaks into xi
if plots,
% Mark starting peak positions with vertical dashed lines
if peakshape==16||peakshape==17
else
for marker=1:NumPeaks,
markx=BestStart((2*marker)-1);
subplot(2,1,1);plot([markx+xoffset markx+xoffset],[0 max(yy)],'m--')
end % for
end % if peakshape
plot(xxx+xoffset,mmodel,'r'); % Plot the total model (sum of component peaks) in red lines
hold off;
axis([min(xx)+xoffset max(xx)+xoffset min(yy) max(yy)]);
switch AUTOZERO,
case 0
title('Peakfit 3.6 Autozero OFF.')
case 1
title('Peakfit 3.6 Linear autozero.')
case 2
title('Peakfit 3.6 Quadratic autozero.')
end
if peakshape==4||peakshape==5||peakshape==8||peakshape==13||peakshape==14||peakshape==15, % Shapes with Extra factor
xlabel(['Peaks = ' num2str(NumPeaks) ' Shape = ' ShapeString ' Error = ' num2str(round(100*LowestError)/100) '% Extra = ' num2str(extra) ] )
else
xlabel(['Peaks = ' num2str(NumPeaks) ' Shape = ' ShapeString ' Error = ' num2str(round(100*LowestError)/100) '%' ] )
end
% Bottom half of the figure shows the residuals and displays RMS error
% between original signal and model
residual=yy-bestmodel;
subplot(2,1,2);plot(xx+xoffset,residual,'b.')
axis([min(xx)+xoffset max(xx)+xoffset min(residual) max(residual)]);
xlabel('Residual Plot')
end % if plots
% Put results into a matrix, one row for each peak, showing peak index number,
% position, amplitude, and width.
for m=1:NumPeaks,
if m==1,
if peakshape==6||peakshape==7||peakshape==8, % equal-width peak models only
FitResults=[[round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]];
else
if peakshape==11||peakshape==12, % Fixed-width shapes only
FitResults=[[round(m) FitParameters(m)+xoffset height(m) FIXEDWIDTH area(m)]];
else
if peakshape==16||peakshape==17, % Fixed-position shapes only
FitResults=[round(m) FIXEDPOSITIONS(m) height(m) FitParameters(m) area(m)];
else
FitResults=[round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)];
end
end
end % if peakshape
else
if peakshape==6||peakshape==7||peakshape==8, % equal-width peak models only
FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]];
else
if peakshape==11||peakshape==12, % Fixed-width shapes only
FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) FIXEDWIDTH area(m)]];
else
if peakshape==16||peakshape==17, % Fixed-position shapes only
FitResults=[FitResults ; [round(m) FIXEDPOSITIONS(m) height(m) FitParameters(m) area(m)]];
else
FitResults=[FitResults ; [round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]];
end
end
end % if peakshape
end % m==1
end % for m=1:NumPeaks
% Display Fit Results on upper graph
if plots,
subplot(2,1,1);
startx=min(xx)+xoffset+(max(xx)-min(xx))./20;
dxx=(max(xx)-min(xx))./10;
dyy=(max(yy)-min(yy))./10;
starty=max(yy)-dyy;
FigureSize=get(gcf,'Position');
if peakshape==9||peakshape==10, % Pulse and sigmoid shapes only
text(startx,starty+dyy/2,['Peak # tau1 Height tau2 Area'] );
else
text(startx,starty+dyy/2,['Peak # Position Height Width Area'] );
end
% Display FitResults using sprintf
for peaknumber=1:NumPeaks,
for column=1:5,
itemstring=sprintf('%0.4g',FitResults(peaknumber,column));
xposition=startx+(1.7.*dxx.*(column-1).*(600./FigureSize(3)));
yposition=starty-peaknumber.*dyy.*(400./FigureSize(4));
text(xposition,yposition,itemstring);
end
end
end % if plots
if NumArgOut==6,
if plots,disp('Computing bootstrap sampling statistics.....'),end
BootstrapResultsMatrix=zeros(5,100,NumPeaks);
BootstrapErrorMatrix=zeros(1,100,NumPeaks);
clear bx by
tic;
for trial=1:100,
n=1;
bx=xx;
by=yy;
while n<length(xx)-1,
if rand>.5,
bx(n)=xx(n+1);
by(n)=yy(n+1);
end
n=n+1;
end
bx=bx+xoffset;
[FitResults,BootFitError]=fitpeaks(bx,by,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,FIXEDWIDTH);
for peak=1:NumPeaks,
BootstrapResultsMatrix(:,trial,peak)=FitResults(peak,:);
BootstrapErrorMatrix(:,trial,peak)=BootFitError;
end
end
if plots,toc;end
for peak=1:NumPeaks,
if plots,
disp(' ')
disp(['Peak #',num2str(peak) ' Position Height Width Area']);
end % if plots
BootstrapMean=mean(real(BootstrapResultsMatrix(:,:,peak)'));
BootstrapSTD=std(BootstrapResultsMatrix(:,:,peak)');
BootstrapIQR=iqr(BootstrapResultsMatrix(:,:,peak)');
PercentRSD=100.*BootstrapSTD./BootstrapMean;
PercentIQR=100.*BootstrapIQR./BootstrapMean;
BootstrapMean=BootstrapMean(2:5);
BootstrapSTD=BootstrapSTD(2:5);
BootstrapIQR=BootstrapIQR(2:5);
PercentRSD=PercentRSD(2:5);
PercentIQR=PercentIQR(2:5);
if plots,
disp(['Bootstrap Mean: ', num2str(BootstrapMean)])
disp(['Bootstrap STD: ', num2str(BootstrapSTD)])
disp(['Bootstrap IQR: ', num2str(BootstrapIQR)])
disp(['Percent RSD: ', num2str(PercentRSD)])
disp(['Percent IQR: ', num2str(PercentIQR)])
end % if plots
BootResults(peak,:)=[BootstrapMean BootstrapSTD PercentRSD BootstrapIQR PercentIQR];
end % peak=1:NumPeaks,
end % if NumArgOut==6,
% ----------------------------------------------------------------------
function [FitResults,LowestError]=fitpeaks(xx,yy,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedparameters)
% Based on peakfit Version 3: June, 2012.
global PEAKHEIGHTS FIXEDWIDTH FIXEDPOSITIONS
format short g
format compact
warning off all
FIXEDWIDTH=fixedparameters;
xoffset=0;
if start==0;start=calcstart(xx,NumPeaks,xoffset);end
PEAKHEIGHTS=zeros(1,NumPeaks);
n=length(xx);
newstart=start;
for peaks=1:NumPeaks,
peakindex=2*peaks-1;
newstart(peakindex)=start(peakindex)-xoffset;
end
% Perform peak fitting for selected peak shape using fminsearch function
options = optimset('TolX',.00001,'Display','off' );
LowestError=1000; % or any big number greater than largest error expected
FitParameters=zeros(1,NumPeaks.*2);
BestStart=zeros(1,NumPeaks.*2);
height=zeros(1,NumPeaks);
bestmodel=zeros(size(yy));
for k=1:NumTrials,
switch peakshape
case 1
TrialParameters=fminsearch(@(lambda)(fitgaussian(lambda,xx,yy)),newstart);
case 2
TrialParameters=fminsearch(@(lambda)(fitlorentzian(lambda,xx,yy)),newstart);
case 3
TrialParameters=fminsearch(@(lambda)(fitlogistic(lambda,xx,yy)),newstart);
case 4
TrialParameters=fminsearch(@(lambda)(fitpearson(lambda,xx,yy,extra)),newstart);
case 5
zxx=[zeros(size(xx)) xx zeros(size(xx)) ];
zyy=[zeros(size(yy)) yy zeros(size(yy)) ];
TrialParameters=fminsearch(@(lambda)(fitexpgaussian(lambda,zxx,zyy,-extra)),newstart);
case 6
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitewgaussian(lambda,xx,yy)),cwnewstart);
case 7
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitlorentziancw(lambda,xx,yy)),cwnewstart);
case 8
cwnewstart(1)=newstart(1);
for pc=2:NumPeaks,
cwnewstart(pc)=newstart(2.*pc-1);
end
cwnewstart(NumPeaks+1)=(max(xx)-min(xx))/5;
TrialParameters=fminsearch(@(lambda)(fitexpewgaussian(lambda,xx,yy,-extra)),cwnewstart);
case 9
TrialParameters=fminsearch(@(lambda)(fitexppulse(lambda,xx,yy)),newstart);
case 10
TrialParameters=fminsearch(@(lambda)(fitsigmoid(lambda,xx,yy)),newstart);
case 11
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFWGaussian(lambda,xx,yy)),fixedstart);
case 12
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFWLorentzian(lambda,xx,yy)),fixedstart);
case 13
TrialParameters=fminsearch(@(lambda)(fitGL(lambda,xx,yy,extra)),newstart);
case 14
TrialParameters=fminsearch(@(lambda)(fitBiGaussian(lambda,xx,yy,extra)),newstart);
case 15
TrialParameters=fminsearch(@(lambda)(fitBiLorentzian(lambda,xx,yy,extra)),newstart);
case 16
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFPGaussian(lambda,xx,yy)),fixedstart,options);
case 17
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@(lambda)(FitFPLorentzian(lambda,xx,yy)),fixedstart,options);
otherwise
end % switch peakshape
% Construct model from Trial parameters
A=zeros(NumPeaks,n);
for m=1:NumPeaks,
switch peakshape
case 1
A(m,:)=gaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 2
A(m,:)=lorentzian(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 3
A(m,:)=logistic(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 4
A(m,:)=pearson(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 5
A(m,:)=expgaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m),-extra)';
case 6
A(m,:)=gaussian(xx,TrialParameters(m),TrialParameters(NumPeaks+1));
case 7
A(m,:)=lorentzian(xx,TrialParameters(m),TrialParameters(NumPeaks+1));
case 8
A(m,:)=expgaussian(xx,TrialParameters(m),TrialParameters(NumPeaks+1),-extra)';
case 9
A(m,:)=exppulse(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 10
A(m,:)=sigmoid(xx,TrialParameters(2*m-1),TrialParameters(2*m));
case 11
A(m,:)=gaussian(xx,TrialParameters(m),FIXEDWIDTH);
case 12
A(m,:)=lorentzian(xx,TrialParameters(m),FIXEDWIDTH);
case 13
A(m,:)=GL(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 14
A(m,:)=BiGaussian(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 15
A(m,:)=BiLorentzian(xx,TrialParameters(2*m-1),TrialParameters(2*m),extra);
case 16
A(m,:)=gaussian(xx,FIXEDPOSITIONS(m),TrialParameters(m));
case 17
A(m,:)=lorentzian(xx,FIXEDPOSITIONS(m),TrialParameters(m));
end % switch
for parameter=1:2:2*NumPeaks,
newstart(parameter)=newstart(parameter)*(1+randn/50);
newstart(parameter+1)=newstart(parameter+1)*(1+randn/10);
end
end % for
% Multiplies each row by the corresponding amplitude and adds them up
model=PEAKHEIGHTS'*A;
% Compare trial model to data segment and compute the fit error
MeanFitError=100*norm(yy-model)./(sqrt(n)*max(yy));
% Take only the single fit that has the lowest MeanFitError
if MeanFitError<LowestError,
if min(PEAKHEIGHTS)>0, % Consider only fits with positive peak heights
LowestError=MeanFitError; % Assign LowestError to the lowest MeanFitError
FitParameters=TrialParameters; % Assign FitParameters to the fit with the lowest MeanFitError
% BestStart=newstart; % Assign BestStart to the start with the lowest MeanFitError
height=PEAKHEIGHTS; % Assign height to the PEAKHEIGHTS with the lowest MeanFitError
% bestmodel=model; % Assign bestmodel to the model with the lowest MeanFitError
end % if min(PEAKHEIGHTS)>0
end % if MeanFitError<LowestError
end % for k (NumTrials)
for m=1:NumPeaks,
area(m)=trapz(xx+xoffset,height(m)*A(m,:)); % Compute the area of each component peak using trapezoidal method
end
for m=1:NumPeaks,
if m==1,
if peakshape==6||peakshape==7||peakshape==8, % equal-width peak models
FitResults=[[round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]];
else
if peakshape==11||peakshape==12, % Fixed-width shapes only
FitResults=[[round(m) FitParameters(m)+xoffset height(m) FIXEDWIDTH area(m)]];
else
FitResults=[[round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]];
end
end % if peakshape
else
if peakshape==6||peakshape==7||peakshape==8, % equal-width peak models
FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]];
else
if peakshape==11||peakshape==12, % Fixed-width shapes only
FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) FIXEDWIDTH area(m)]];
else
FitResults=[FitResults ; [round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]];
end
end % if peakshape
end % m==1
end % for m=1:NumPeaks
% ----------------------------------------------------------------------
function start=calcstart(xx,NumPeaks,xoffset)
n=max(xx)-min(xx);
start=[];
startpos=[n/(NumPeaks+1):n/(NumPeaks+1):n-(n/(NumPeaks+1))]+min(xx);
for marker=1:NumPeaks,
markx=startpos(marker)+ xoffset;
start=[start markx n/ (3.*NumPeaks)];
end % for marker
% ----------------------------------------------------------------------
function [index,closestval]=val2ind(x,val)
% Returns the index and the value of the element of vector x that is closest to val
% If more than one element is equally close, returns vectors of indicies and values
% Tom O'Haver (toh@umd.edu) October 2006
% Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);
% ----------------------------------------------------------------------
function err = fitgaussian(lambda,t,y)
% Fitting function for a Gaussian band signal.
global PEAKHEIGHTS
numpeaks=round(length(lambda)/2);
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = gaussian(t,lambda(2*j-1),lambda(2*j))';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = fitewgaussian(lambda,t,y)
% Fitting function for a Gaussian band signal with equal peak widths.
global PEAKHEIGHTS
numpeaks=round(length(lambda)-1);
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = gaussian(t,lambda(j),lambda(numpeaks+1))';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = FitFWGaussian(lambda,t,y)
% Fitting function for a fixed width Gaussian
global PEAKHEIGHTS FIXEDWIDTH
numpeaks=round(length(lambda));
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = gaussian(t,lambda(j),FIXEDWIDTH)';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = FitFPGaussian(lambda,t,y)
% Fitting function for fixed-position Gaussians
global PEAKHEIGHTS FIXEDPOSITIONS
numpeaks=round(length(lambda));
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = gaussian(t,FIXEDPOSITIONS(j), lambda(j))';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = FitFPLorentzian(lambda,t,y)
% Fitting function for fixed-position Lorentzians
global PEAKHEIGHTS FIXEDPOSITIONS
numpeaks=round(length(lambda));
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = lorentzian(t,FIXEDPOSITIONS(j), lambda(j))';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = FitFWLorentzian(lambda,t,y)
% Fitting function for fixed width Gaussians
global PEAKHEIGHTS FIXEDWIDTH
numpeaks=round(length(lambda));
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = lorentzian(t,lambda(j),FIXEDWIDTH)';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = fitlorentziancw(lambda,t,y)
% Fitting function for a Lorentzian band signal with equal peak widths.
global PEAKHEIGHTS
numpeaks=round(length(lambda)-1);
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = lorentzian(t,lambda(j),lambda(numpeaks+1))';
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = gaussian(x,pos,wid)
% gaussian(X,pos,wid) = gaussian peak centered on pos, half-width=wid
% X may be scalar, vector, or matrix, pos and wid both scalar
% Examples: gaussian([0 1 2],1,2) gives result [0.5000 1.0000 0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)).^2);
% ----------------------------------------------------------------------
function err = fitlorentzian(lambda,t,y)
% Fitting function for single lorentzian, lambda(1)=position, lambda(2)=width
% Fitgauss assumes a lorentzian function
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = lorentzian(t,lambda(2*j-1),lambda(2*j))';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = lorentzian(x,position,width)
% lorentzian(x,position,width) Lorentzian function.
% where x may be scalar, vector, or matrix
% position and width scalar
% T. C. O'Haver, 1988
% Example: lorentzian([1 2 3],2,2) gives result [0.5 1 0.5]
g=ones(size(x))./(1+((x-position)./(0.5.*width)).^2);
% ----------------------------------------------------------------------
function err = fitlogistic(lambda,t,y)
% Fitting function for logistic, lambda(1)=position, lambda(2)=width
% between the data and the values computed by the current
% function of lambda. Fitlogistic assumes a logistic function
% T. C. O'Haver, May 2006
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = logistic(t,lambda(2*j-1),lambda(2*j))';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = logistic(x,pos,wid)
% logistic function. pos=position; wid=half-width (both scalar)
% logistic(x,pos,wid), where x may be scalar, vector, or matrix
% pos=position; wid=half-width (both scalar)
% T. C. O'Haver, 1991
n = exp(-((x-pos)/(.477.*wid)) .^2);
g = (2.*n)./(1+n);
% ----------------------------------------------------------------------
function err = fitlognormal(lambda,t,y)
% Fitting function for lognormal, lambda(1)=position, lambda(2)=width
% between the data and the values computed by the current
% function of lambda. Fitlognormal assumes a lognormal function
% T. C. O'Haver, May 2006
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = lognormal(t,lambda(2*j-1),lambda(2*j))';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = lognormal(x,pos,wid)
% lognormal function. pos=position; wid=half-width (both scalar)
% lognormal(x,pos,wid), where x may be scalar, vector, or matrix
% pos=position; wid=half-width (both scalar)
% T. C. O'Haver, 1991
g = exp(-(log(x/pos)/(0.01.*wid)) .^2);
% ----------------------------------------------------------------------
function err = fitpearson(lambda,t,y,shapeconstant)
% Fitting functions for a Pearson 7 band signal.
% T. C. O'Haver (toh@umd.edu), Version 1.3, October 23, 2006.
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = pearson(t,lambda(2*j-1),lambda(2*j),shapeconstant)';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = pearson(x,pos,wid,m)
% Pearson VII function.
% g = pearson7(x,pos,wid,m) where x may be scalar, vector, or matrix
% pos=position; wid=half-width (both scalar)
% m=some number
% T. C. O'Haver, 1990
g=ones(size(x))./(1+((x-pos)./((0.5.^(2/m)).*wid)).^2).^m;
% ----------------------------------------------------------------------
function err = fitexpgaussian(lambda,t,y,timeconstant)
% Fitting functions for a exponentially-broadened Gaussian band signal.
% T. C. O'Haver, October 23, 2006.
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = expgaussian(t,lambda(2*j-1),lambda(2*j),timeconstant);
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function err = fitexpewgaussian(lambda,t,y,timeconstant)
% Fitting function for exponentially-broadened Gaussian bands with equal peak widths.
global PEAKHEIGHTS
numpeaks=round(length(lambda)-1);
A = zeros(length(t),numpeaks);
for j = 1:numpeaks,
A(:,j) = expgaussian(t,lambda(j),lambda(numpeaks+1),timeconstant);
end
PEAKHEIGHTS = abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = expgaussian(x,pos,wid,timeconstant)
% Exponentially-broadened gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
% x may be scalar, vector, or matrix, pos and wid both scalar
% T. C. O'Haver, 2006
g = exp(-((x-pos)./(0.6005615.*wid)) .^2);
g = ExpBroaden(g',timeconstant);
% ----------------------------------------------------------------------
function yb = ExpBroaden(y,t)
% ExpBroaden(y,t) zero pads y and convolutes result by an exponential decay
% of time constant t by multiplying Fourier transforms and inverse
% transforming the result.
hly=round(length(y)./2);
ey=[zeros(1,hly)';y;zeros(1,hly)'];
fy=fft(ey);
a=exp(-(1:length(fy))./t);
fa=fft(a);
fy1=fy.*fa';
ybz=real(ifft(fy1))./sum(a);
yb=ybz(hly+2:length(ybz)-hly+1);
% ----------------------------------------------------------------------
function err = fitexppulse(tau,x,y)
% Iterative fit of the sum of exponental pulses
% of the form Height.*exp(-tau1.*x).*(1-exp(-tau2.*x)))
global PEAKHEIGHTS
A = zeros(length(x),round(length(tau)/2));
for j = 1:length(tau)/2,
A(:,j) = exppulse(x,tau(2*j-1),tau(2*j));
end
PEAKHEIGHTS =abs(A\y');
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = exppulse(x,t1,t2)
% Exponential pulse of the form
% Height.*exp(-tau1.*x).*(1-exp(-tau2.*x)))
e=(x-t1)./t2;
p = 4*exp(-e).*(1-exp(-e));
p=p .* [p>0];
g = p';
% ----------------------------------------------------------------------
function err = fitsigmoid(tau,x,y)
% Fitting function for iterative fit to the sum of
% sigmiods of the form Height./(1 + exp((t1 - t)/t2))
global PEAKHEIGHTS
A = zeros(length(x),round(length(tau)/2));
for j = 1:length(tau)/2,
A(:,j) = sigmoid(x,tau(2*j-1),tau(2*j));
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g=sigmoid(x,t1,t2)
% g=1./(1 + exp((t1 - x)./t2))'; % old version of sigmoid
g=1/2 + 1/2* erf(real((x-t1)/sqrt(2*t2))); % Modified sigmoid
% Bifurcated sigmoid
% lx=length(x);
% hx=val2ind(x,t1);
% g(1:hx)=1./(1 + exp((t1 - x(1:hx))./t2));
% g(hx+1:lx)=1./(1 + exp((t1 - x(hx+1:lx))./(1.3*t2)));
% ----------------------------------------------------------------------
function err = fitGL(lambda,t,y,shapeconstant)
% Fitting functions for Gaussian/Lorentzian blend.
% T. C. O'Haver (toh@umd.edu), 2012.
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = GL(t,lambda(2*j-1),lambda(2*j),shapeconstant)';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = GL(x,pos,wid,m)
% Gaussian/Lorentzian blend. m = percent Gaussian character
% pos=position; wid=half-width
% m = percent Gaussian character.
% T. C. O'Haver, 2012
g=((m/100)*gaussian(x,pos,wid)+(1-(m/100))*lorentzian(x,pos,wid))/2;
% ----------------------------------------------------------------------
function err = fitBiGaussian(lambda,t,y,shapeconstant)
% Fitting functions for BiGaussian.
% T. C. O'Haver (toh@umd.edu), 2012.
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = BiGaussian(t,lambda(2*j-1),lambda(2*j),shapeconstant)';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = BiGaussian(x,pos,wid,m)
% BiGaussian (different widths on leading edge and trailing edge).
% pos=position; wid=width
% m determines shape; symmetrical if m=50.
% T. C. O'Haver, 2012
lx=length(x);
hx=val2ind(x,pos);
g(1:hx)=gaussian(x(1:hx),pos,wid*(m/100));
g(hx+1:lx)=gaussian(x(hx+1:lx),pos,wid*(1-m/100));
% ----------------------------------------------------------------------
function err = fitBiLorentzian(lambda,t,y,shapeconstant)
% Fitting functions for BiGaussian.
% T. C. O'Haver (toh@umd.edu), 2012.
global PEAKHEIGHTS
A = zeros(length(t),round(length(lambda)/2));
for j = 1:length(lambda)/2,
A(:,j) = BiLorentzian(t,lambda(2*j-1),lambda(2*j),shapeconstant)';
end
PEAKHEIGHTS = A\y';
z = A*PEAKHEIGHTS;
err = norm(z-y');
% ----------------------------------------------------------------------
function g = BiLorentzian(x,pos,wid,m)
% BiLorentzian (different widths on leading edge and trailing edge).
% pos=position; wid=width
% m determines shape; symmetrical if m=50.
% T. C. O'Haver, 2012
lx=length(x);
hx=val2ind(x,pos);
g(1:hx)=lorentzian(x(1:hx),pos,wid*(m/100));
g(hx+1:lx)=lorentzian(x(hx+1:lx),pos,wid*(1-m/100));
% ----------------------------------------------------------------------
function b=iqr(a)
% b = IQR(a) returns the interquartile range of the values in a. For
% vector input, b is the difference between the 75th and 25th percentiles
% of a. For matrix input, b is a row vector containing the interquartile
% range of each column of a.
% T. C. O'Haver, 2012
mina=min(a);
sizea=size(a);
NumCols=sizea(2);
for n=1:NumCols,b(:,n)=a(:,n)-mina(n);end
Sorteda=sort(b);
lx=length(Sorteda);
SecondQuartile=round(lx/4);
FourthQuartile=3*round(lx/4);
b=abs(Sorteda(FourthQuartile,:)-Sorteda(SecondQuartile,:)); |
|
