Code covered by the BSD License

# Interactive Peak Fitter (Version 9.7)

27 Mar 2009 (Updated 11 Nov 2013)

Keyboard operated peak fitting function for time-series signals.

DemoPeakfitBootstrap
```function DemoPeakfitBootstrap
% Demonstration script for peakfit.m, with built-in signal generator.
% Demonstrates bootstrap error estimation for Position, Height, and Width.
% You can change the character of the simulated signal in lines 10-16.
% (c) T. C. O'Haver (toh@umd.edu).  September 2012.
format short g
format compact
warning off all
% Generate simulated signal
n=2000;  % maximum x-value
increment=2; % Difference between adjacent x values
x=1:increment:n;
amp=[1 2 3 1 4 1 3 2 1 2 1];  % Amplitudes of the peaks (one entry for each peak)
pos=[100 130 200 400 600 630 800 850 900 1200 1700];   % Positions of the peaks (one entry for each peak)
wid=[30 40 50 30 20 20 30 40 50 100 200];   % Widths of the peaks (one entry for each peak)
noise=.05 ;
% A = matrix containing one of the unit-amplidude functions in each of its rows
clear A
A=zeros(length(pos),length(x));
for k=1:length(pos),
A(k,:)=exp(-((x-pos(k))./(0.6005612.*wid(k))) .^2);  % Gaussian function
area(k)=amp(k).*trapz(x,A(k,:)); % Computes peak areas via trapaziod method
end
TrueY=amp*A;  % Multiplies each row by the corresponding amplitude and adds them up
% TrueY=TrueY+lorentzian(x,0,1000); % Adds sloping curved baseline
y = TrueY + noise .* randn(size(x));

% Arrange the peak parameters of the simulated signal into a matrix
ColumnLabels='           Peak#    Position       Height      Width       Area';
NumPeaks=length(pos);
for m=1:NumPeaks,
if m==1,
ActualParameters=[m pos(m) amp(m) wid(m) area(m) ];
else
ActualParameters=[ActualParameters; [m pos(m) amp(m) wid(m) area(m) ]];
end
end
disp('-----------------------------------------------------------------')
disp('Demonstration  of the peakfit function with error estimation')
disp('by the bootstrap method. (You can change the character of the')
disp('simulated signal in lines 10 - 16) ')

% Display the peak parameters of  the simulated signal
disp('True peak parameters of the simulated signal:')
ColumnLabels
ActualParameters
clf
plotfit(x,y);
fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
clf
disp(' ')
% Execute peakfit
disp('Dividing the signal into 5 groups of overlapping peaks and ')
disp('applying peakfit.m to each group: ')
disp(' ')
disp('GROUP 1 (3 peaks)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],158,314,3,1,0,1,[99 29 129 41 200 50],0,0); ')
disp('            Peak      Position    Height      Width       Area')
TrueParameters=ActualParameters(1:3,:)
FitResults=peakfit([x;y],158,314,3,1,1,1,[99 29 129 41 200 50],0,0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],158,314,3,1,1,1,[99 29 129 41 200 50],0,0);
disp(' ')
drawnow
fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
clf

disp('GROUP 2 (1 peak)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],403,144,1,1,0,1,0,0); ')
disp('            Peak      Position    Height      Width       Area')
TrueParameters=ActualParameters(4,:)
FitResults=peakfit([x;y],403,144,1,1,0,1,0,0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],403,144,1,1,0,1,0,0);
disp(' ')

fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
disp('GROUP 3 (2 peaks)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],612,110,2,1,0,1,0,0);')
disp('            Peak      Position    Height      Width       Area')
TrueParameters=ActualParameters(5:6,:)
FitResults=peakfit([x;y],612,110,2,1,0,1,0,0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],612,110,2,1,0,1,0,0);
disp(' ')

fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
clf
disp('GROUP 4 (3 peaks)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],847,224,3,1,0,1,[800 30 849 39 899 52],0);')
TrueParameters=ActualParameters(7:9,:)
disp('            Peak      Position    Height      Width       Area')
FitResults=peakfit([x;y],847,224,3,1,0,1,[800 30 849 39 899 52],0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],847,224,3,1,0,1,[800 30 849 39 899 52],0);
disp(' ')

fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
clf
disp('GROUP 5 (1 peak)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],1196,358,1,1,0,1,0,0);')
disp('            Peak      Position    Height      Width       Area')
TrueParameters=ActualParameters(10,:)
FitResults=peakfit([x;y],1196,358,1,1,0,1,0,0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],1196,358,1,1,0,1,0,0);
disp('>> ')

fprintf(2,'Press any key to continue.....or press Ctrl-C to abort.\n')
pause
clf
disp('GROUP 6 (1 peak)')
disp('>> [FitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],1702,594,1,1,0,1,0,0);')
disp('            Peak      Position    Height      Width       Area')
TrueParameters=ActualParameters(11,:)
FitResults=peakfit([x;y],1702,594,1,1,0,1,0,0)
drawnow
[bFitResults,MeanFitError,BestStart,xi,yi,BootResults]=peakfit([x;y],1702,594,1,1,0,1,0,0);

disp('  ')
disp('HOW TO USE THE BOOTSTRAP RESULTS ')
disp('When the 6th output argument is provided (BootResults), peakfit')
disp('computes the expected precision of the measured peak parameters')
disp('by the bootstrap method and places the standard deviations of')
disp('the peak position, height, width, and area in BootResults(5:8).')
disp('In most cases, the measured peak parameters should fall within')
disp('plus of minus 2 standard deviations (STD) of the correct result.')
disp('  ')
disp(['For example, in the last fit, the true value of the peak position is ' num2str(TrueParameters(2)) ',' ])
disp(['and the peakfit function measured the peak position as ' num2str(FitResults(2) ) ', with a  ' ] )
disp(['predicted precision of plus or minus ' num2str(2.*BootResults(5) ) ' (twice the standard deviation).' ])
disp('  ')
disp(['The true value of the peak height is ' num2str(TrueParameters(3)) ' and the peakfit function,' ])
disp(['measured the peak height as ' num2str(FitResults(3) ) ', with a predicted precision of ' ] )
disp(['plus or minus ' num2str(2.*BootResults(6) ) '.' ])
disp('  ')
disp(['The true value of the peak width is ' num2str(TrueParameters(4)) ' and the peakfit function,' ])
disp(['measured the peak width as ' num2str(FitResults(4) ) ', with a predicted precision of ' ] )
disp(['plus or minus ' num2str(2.*BootResults(7) ) '.' ])
disp('  ')
disp(['The true value of the peak area is ' num2str(TrueParameters(5)) ' and the peakfit function,' ])
disp(['measured the peak area as ' num2str(FitResults(5) ) ', with a predicted precision of ' ] )
disp(['plus or minus ' num2str(2.*BootResults(8) ) '.' ])
disp('  ')
disp('Note: If the reported value of the standard deviation (STD) is ')
disp('much greater than the IQR (InterQuartile Range), then a better')
disp('estimate of the standard deviation is IQR/1.34896. ')
disp('  ')
fprintf(2,'End of demo.\n')
disp('If you run this again, you''ll get slightly different results')
disp('because random noise is added to the signal each time. You may')
disp('change the nature of the simulated signal in lines 10 - 16.')

% Internal functions
% ----------------------------------------------------------------------
function [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedwidth,plots)
% Version 3.4: October, 2012. Works in Matlab or Octave 3.6.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

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;
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';
plots=1;

% 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;
case 2
NumPeaks=1;
peakshape=1;
extra=0;
NumTrials=1;
xx=signal;yy=center;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
case 3
NumPeaks=1;
peakshape=1;
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
case 4
peakshape=1;
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
case 5
extra=0;
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
case 6
NumTrials=1;
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
case 7
start=calcstart(xx,NumPeaks,xoffset);
AUTOZERO=1;
FIXEDWIDTH=0;
case 8
AUTOZERO=1;
FIXEDWIDTH=0;
case 10
FIXEDWIDTH=fixedwidth;
case 11
FIXEDWIDTH=fixedwidth;
plots=0;
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

% Remove linear baseline from data segment if AUTOZERO==1
X1=min(xx);
X2=max(xx);
bkgsize=round(length(xx)/10);
if bkgsize<2,bkgsize=2;end
if AUTOZERO==1, % linear autozero operation
Y1=mean(yy(1:bkgsize));
Y2=mean(yy((length(xx)-bkgsize):length(xx)));
yy=yy-((Y2-Y1)/(X2-X1)*(xx-X1)+Y1);
end % if
lxx=length(xx);
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';
otherwise
end % switch peakshape

% 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,
% 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);
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);
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);
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
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);
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
for marker=1:NumPeaks,
markx=BestStart((2*marker)-1);
subplot(2,1,1);plot([markx+xoffset markx+xoffset],[0 max(yy)],'m--')
end % for
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.3 Autozero OFF.')
case 1
title('Peakfit 3.3 Linear autozero.')
case 2
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
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 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
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

% 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,fixedwidth)
% Based on peakfit Version 3: June, 2012.
global PEAKHEIGHTS FIXEDWIDTH
format short g
format compact
warning off all
FIXEDWIDTH=fixedwidth;
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);
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);
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/5];
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 = FitFWLorentzian(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) = 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);
% ----------------------------------------------------------------------
% ExpBroaden(y,t) convolutes y by an exponential decay of time constant t
% by multiplying Fourier transforms and inverse transforming the result.
ly=length(y);
ey=[zeros(size(y));y;zeros(size(y))];
fy=fft(ey);
a=exp(-(1:length(fy))./t);
fa=fft(a);
fy1=fy.*fa';
ybz=real(ifft(fy1))./sum(a);
yb=ybz(ly+2:length(ybz)-ly+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))';
% ----------------------------------------------------------------------
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,:));```