Code covered by the BSD License

Highlights from Demo functions for peak detection and fitting.

• DemoPeakfitBootstrap Demonstration script for peakfit.m, with built-in signal generator.
• PeakFitDemo11Lor Demonstration script for ipf.m, with built-in signal generator.
• idemo Self-contained demonstration function for comparing the ipeak and
• idemo1 Self-contained demonstration function for comparing the iPeak and
• idemo2 Self-contained demonstration function for comparing the iPeak and.
• View all files

Demo functions for peak detection and fitting.

Tom O'Haver (view profile)

04 Sep 2012 (Updated )

A collection of self-contained demonstration functions for iPeak.m and peakfit.m.

DemoPeakfitBootstrap
function DemoPeakfitBootstrap
% Demonstration script for peakfit.m, with built-in signal generator.
% Demonstrates bootstrap error estimation.
% 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.6006.*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 deviation (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 Noise level in line 16 if you wish.')

% ----------------------------------------------------------------------

function [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedwidth)
% Version 3: August, 2012. Added parameter error estimation by bootstrap.
% 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';

% 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;
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;
switch NumPeaks
case 1
newstart(1)=start(1)-xoffset;
case 2
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
case 3
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
case 4
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
case 5
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
newstart(9)=start(9)-xoffset;
case 6
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
newstart(9)=start(9)-xoffset;
newstart(11)=start(11)-xoffset;
otherwise
end % switch NumPeaks

% 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';
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(@fitgaussian,newstart,options,xx,yy);
case 2
TrialParameters=fminsearch(@fitlorentzian,newstart,options,xx,yy);
case 3
TrialParameters=fminsearch(@fitlogistic,newstart,options,xx,yy);
case 4
TrialParameters=fminsearch(@fitpearson,newstart,options,xx,yy,extra);
case 5
zxx=[zeros(size(xx)) xx zeros(size(xx)) ];
zyy=[zeros(size(yy)) yy zeros(size(yy)) ];
TrialParameters=fminsearch(@fitexpgaussian,newstart,options,zxx,zyy,-extra);n';
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(@fitewgaussian,cwnewstart,options,xx,yy);
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(@fitlorentziancw,cwnewstart,options,xx,yy);
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(@fitexpewgaussian,cwnewstart,options,xx,yy,-extra);
case 9
TrialParameters=fminsearch(@fitexppulse,newstart,options,xx,yy);
case 10
TrialParameters=fminsearch(@fitsigmoid,newstart,options,xx,yy);
case 11
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@FitFWGaussian,fixedstart,options,xx,yy);
case 12
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@FitFWLorentzian,fixedstart,options,xx,yy);
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);
otherwise
end % switch
switch NumPeaks % adds random variation to non-linear parameters
case 1
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10)];
case 2
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10)];
case 3
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)];
case 4
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)  newstart(7)*(1+randn/50) newstart(8)*(1+randn/10)];
case 5
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)  newstart(7)*(1+randn/50) newstart(8)*(1+randn/10)  newstart(9)*(1+randn/50) newstart(10)*(1+randn/10)];
otherwise
end % switch NumPeaks
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,100);
xxx=linspace(min(xx),max(xx));
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);
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.
subplot(2,1,1);plot(xx+xoffset,yy,'b.'); % Plot the original signal in blue dots
hold on
for m=1:NumPeaks,
plot(xxx+xoffset,height(m)*AA(m,:),'g')  % 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

% 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. Autozero OFF.')
case 1
title('Peakfit 3. Linear autozero.')
case 2
end
if peakshape==4||peakshape==5||peakshape==8, % 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')

% 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
FitResults=[[round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]];
else
if peakshape==11||peakshape==12,
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,
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
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,
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
if NumArgOut==6,
disp('Computing bootstrap sampling statistics.....')
BootstrapResultsMatrix=zeros(5,100,NumPeaks);
BootstrapErrorMatrix=zeros(1,100,NumPeaks);
clear bx by
cutoff=0.5;
% offset=xoffset
% sizexx=size(xx)
% sizeyy=size(yy)
tic;
for trial=1:100,
n=1;
bx=xx;
by=yy;
while n<length(xx)-1,
if rand>cutoff,
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
toc;
for peak=1:NumPeaks,
disp(' ')
disp(['Peak #',num2str(peak) '         Position    Height       Width       Area']);
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);
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)])
BootResults(peak,:)=[BootstrapMean BootstrapSTD PercentRSD BootstrapIQR PercentIQR];
end
end
% ----------------------------------------------------------------------
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;
switch NumPeaks
case 1
newstart(1)=start(1)-xoffset;
case 2
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
case 3
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
case 4
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
case 5
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
newstart(9)=start(9)-xoffset;
case 6
newstart(1)=start(1)-xoffset;
newstart(3)=start(3)-xoffset;
newstart(5)=start(5)-xoffset;
newstart(7)=start(7)-xoffset;
newstart(9)=start(9)-xoffset;
newstart(11)=start(11)-xoffset;
otherwise
end % switch NumPeaks

% 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(@fitgaussian,newstart,options,xx,yy);
case 2
TrialParameters=fminsearch(@fitlorentzian,newstart,options,xx,yy);
case 3
TrialParameters=fminsearch(@fitlogistic,newstart,options,xx,yy);
case 4
TrialParameters=fminsearch(@fitpearson,newstart,options,xx,yy,extra);
case 5
zxx=[zeros(size(xx)) xx zeros(size(xx)) ];
zyy=[zeros(size(yy)) yy zeros(size(yy)) ];
TrialParameters=fminsearch(@fitexpgaussian,newstart,options,zxx,zyy,-extra);n';
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(@fitewgaussian,cwnewstart,options,xx,yy);
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(@fitlorentziancw,cwnewstart,options,xx,yy);
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(@fitexpewgaussian,cwnewstart,options,xx,yy,-extra);
case 9
TrialParameters=fminsearch(@fitexppulse,newstart,options,xx,yy);
case 10
TrialParameters=fminsearch(@fitsigmoid,newstart,options,xx,yy);
case 11
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@FitFWGaussian,fixedstart,options,xx,yy);
case 12
fixedstart=[];
for pc=1:NumPeaks,
fixedstart(pc)=min(xx)+pc.*(max(xx)-min(xx))./(NumPeaks+1);
end
TrialParameters=fminsearch(@FitFWLorentzian,fixedstart,options,xx,yy);
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);
otherwise
end % switch
switch NumPeaks % adds random variation to non-linear parameters
case 1
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10)];
case 2
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10)];
case 3
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)];
case 4
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)  newstart(7)*(1+randn/50) newstart(8)*(1+randn/10)];
case 5
newstart=[newstart(1)*(1+randn/50) newstart(2)*(1+randn/10) newstart(3)*(1+randn/50) newstart(4)*(1+randn/10) newstart(5)*(1+randn/50) newstart(6)*(1+randn/10)  newstart(7)*(1+randn/50) newstart(8)*(1+randn/10)  newstart(9)*(1+randn/50) newstart(10)*(1+randn/10)];
otherwise
end % switch NumPeaks
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,
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,
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 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.
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,:));