function P=findpeaksSb(x,y,SlopeThreshold,AmpThreshold,smoothwidth,peakgroup,smoothtype,window,PeakShape,extra,NumTrials,AUTOZERO,fixedparameters) % function P=findpeaksb(x,y,SlopeThreshold,AmpThreshold,smoothwidth, % peakgroup,smoothtype,window,PeakShape,extra,NumTrials,AUTOZERO) % Version 1, November 2016. Function to locate the positive peaks % in a noisy x-y time series data set and fit each peak to a selected peak % model with variable window width, peak shape, and optional background % subtraction. SlopeThreshold, AmpThreshold, smoothwidth, peakgroup, window % width, PeakShape, extra, NumTrials, AUTOZERO, and fixedparameters can all % optionally be scalars or vectors with one entry for each segment. % Returns a matrix P listing the peak number, position, height, width, area, % percent fitting error and "R2" of each detected peak. % See http://terpconnect.umd.edu/~toh/spectrum/Smoothing.html and % http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm % (c) T.C. O'Haver, 2016. Version 1.1, November, 2016 % % Example 1: Narrow Gaussian followed by much broader Gaussian, both with % peak height=1, zero baseline. % x=1:.2:100; % y=.1+gaussian(x,20,1.5)+gaussian(x,80,30)+.02.*randn(size(x)); % plot(x,y,'b.'); % disp(' Peak # Position Height Width Area') % findpeaksSb(x,y,[0.001 .00005],[.6 .6],[5 70],[8 150],3,200,1,0,1,3) % % Example 2: Narrow Lorentzian peak on sloping baseline, followed by much % broader Gaussian peak with no baseline, both height=1. % x=1:.2:100; % y=.2.*gaussian(x,0,30)+lorentzian(x,20,1.5)+gaussian(x,80,30)+.02.*randn(size(x)); % plot(x,y,'b.'); % disp(' Peak # Position Height Width Area') % findpeaksSb(x,y,[0.001 .00005],[.6 .6],[5 120],[8 150],3,[30 200],[2 1],0,1,[3 0]) % % Related functions: % findpeaksG.m, findpeaksSG.m, findvalleys.m, findpeaksL.m, findpeaksb.m, % findpeaksb3.m, findpeaksplot.m, peakstats.m, findpeakSNR.m, % findpeaksGSS.m, findpeaksLSS.m, findpeaksfit.m, findsteps.m, % findsquarepulse.m, idpeaks.m % Copyright (c) 2016 Thomas C. O'Haver % Permission is hereby granted, free of charge, to any person obtaining a copy % of this software and associated documentation files (the "Software"), to deal % in the Software without restriction, including without limitation the rights % to use, copy, modify, merge, publish, distribute, sublicense, and/or sell % copies of the Software, and to permit persons to whom the Software is % furnished to do so, subject to the following conditions: % % The above copyright notice and this permission notice shall be included in % all copies or substantial portions of the Software. % % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR % IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, % FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE % AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER % LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, % OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN % THE SOFTWARE. % sizex=size(x);if sizex(1)>sizex(2),x=x';end sizey=size(y);if sizey(1)>sizey(2),y=y';end NumPeaks=1; % Number of models shapes to you to fit each peak (usually 1) start=0; fixedparameters=0; % Used for fixed parameters shapes only (11,16,12,and 17) plots=0; % Supress plotting if nargin~=7;smoothtype=1;end % smoothtype=1 if not specified in argument if smoothtype>3;smoothtype=3;end if smoothtype<1;smoothtype=1;end if smoothwidth<1;smoothwidth=1;end if isscalar(AmpThreshold),AmpThreshold=AmpThreshold.*ones(size(SlopeThreshold));end if isscalar(smoothwidth),smoothwidth=smoothwidth.*ones(size(SlopeThreshold));end if isscalar(peakgroup),peakgroup=peakgroup.*ones(size(SlopeThreshold));end if isscalar(window),window=window.*ones(size(SlopeThreshold));end if isscalar(PeakShape),PeakShape=PeakShape.*ones(size(SlopeThreshold));end if isscalar(extra),extra=extra.*ones(size(SlopeThreshold));end if isscalar(NumTrials),NumTrials=NumTrials.*ones(size(SlopeThreshold));end if isscalar(AUTOZERO),AUTOZERO=AUTOZERO.*ones(size(SlopeThreshold));end if isscalar(fixedparameters),fixedparameters=fixedparameters.*ones(size(SlopeThreshold));end smoothwidth=round(smoothwidth); peakgroup=round(peakgroup); if smoothwidth>1, d=SegmentedSmooth(deriv(y),smoothwidth,smoothtype); else d=deriv(y); end P=[0 0 0 0 0 0 0]; vectorlength=length(y); NumSegs=length(SlopeThreshold); peak=1; for j=2*round(smoothwidth/2)-1:length(y)-smoothwidth-1, Seg=fix(1+NumSegs./(vectorlength./j)); n=round(peakgroup(Seg)/2+1); % [j Seg] if sign(d(j)) > sign (d(j+1)), % Detects zero-crossing if d(j)-d(j+1) > SlopeThreshold(Seg), % if slope of derivative is larger than SlopeThreshold if y(j) > AmpThreshold(Seg), % if height of peak is larger than AmpThreshold xx=zeros(size(peakgroup(Seg)));yy=zeros(size(peakgroup(Seg))); for k=1:peakgroup(Seg), % Create sub-group of points near peak groupindex=j+k-n+2; if groupindex<1, groupindex=1;end if groupindex>vectorlength, groupindex=vectorlength;end % groupindex=groupindex xx(k)=x(groupindex); yy(k)=y(groupindex); end if peakgroup(Seg)>2, startpoint=round(j-window(Seg)/2); if startpoint<1;startpoint=1;end endpoint=round(j+window(Seg)/2); if endpoint>length(x);endpoint=length(x); end XXX=x(startpoint:endpoint); YYY=y(startpoint:endpoint); % [startpoint endpoint] signal=[XXX',YYY']; % plots=1; % <<<<< Remove after testing % figure(1)% <<<<< Remove after testing % clf % plot(x(startpoint:endpoint),y(startpoint:endpoint)) % <<<<< Remove after testing % figure(2)% <<<<< Remove after testing [FitResults,Error]=peakfit(signal,x(j),window(Seg),NumPeaks,PeakShape(Seg),extra(Seg),NumTrials(Seg),start,AUTOZERO(Seg),fixedparameters(Seg),plots); title(['Peak number ' num2str(peak)]); % pause % <<<<< Remove after testing PeakX=FitResults(2); PeakY=FitResults(3); MeasuredWidth=FitResults(4); else PeakY=max(yy); pindex=val2ind(yy,PeakY); PeakX=xx(pindex(1)); MeasuredWidth=0; end % Construct matrix P. One row for each peak detected, % containing the peak number, peak position (x-value) and % peak height (y-value). If peak measurement fails and % results in NaN, or if the measured peak height is less % than AmpThreshold, skip this peak if isnan(PeakX) || isnan(PeakY) || PeakYly,endindix=ly;end indexrange=startindex:endindix; SmoothedSignal(indexrange)=SmoothSegment(indexrange,Segment); end function SmoothY=fastsmooth(Y,w,type,ends) % fastsmooth(Y,w,type,ends) smooths vector Y with smooth % of width w. Version 3.0, October 2016. % The argument "type" determines the smooth type: % If type=1, rectangular (sliding-average or boxcar) % If type=2, triangular (2 passes of sliding-average) % If type=3, pseudo-Gaussian (3 passes of sliding-average) % If type=4, pseudo-Gaussian (4 passes of same sliding-average) % If type=5, multiple-width (4 passes of different sliding-average) % The argument "ends" controls how the "ends" of the signal % (the first w/2 points and the last w/2 points) are handled. % If ends=0, the ends are zero. (In this mode the elapsed % time is independent of the smooth width). The fastest. % If ends=1, the ends are smoothed with progressively % smaller smooths the closer to the end. (In this mode the % elapsed time increases with increasing smooth widths). % fastsmooth(Y,w,type) smooths with ends=0. % fastsmooth(Y,w) smooths with type=1 and ends=0. % Examples: % fastsmooth([1 1 1 10 10 10 1 1 1 1],3)= [0 1 4 7 10 7 4 1 1 0] % % fastsmooth([1 1 1 10 10 10 1 1 1 1],3,1,1)= [1 1 4 7 10 7 4 1 1 1] % % x=1:100; % y=randn(size(x)); % plot(x,y,x,fastsmooth(y,5,3,1),'r') % xlabel('Blue: white noise. Red: smoothed white noise.') % % Copyright (c) 2012, Thomas C. O'Haver % switch type case 1 SmoothY=sa(Y,w,ends); case 2 SmoothY=sa(sa(Y,w,ends),w,ends); case 3 SmoothY=sa(sa(sa(Y,w,ends),w,ends),w,ends); case 4 SmoothY=sa(sa(sa(sa(Y,w,ends),w,ends),w,ends),w,ends); case 5 SmoothY=sa(sa(sa(sa(Y,round(1.6*w),ends),round(1.4*w),ends),round(1.2*w),ends),w,ends); end function SmoothY=sa(Y,smoothwidth,ends) w=round(smoothwidth); SumPoints=sum(Y(1:w)); s=zeros(size(Y)); halfw=round(w/2); L=length(Y); for k=1:L-w, s(k+halfw-1)=SumPoints; SumPoints=SumPoints-Y(k); SumPoints=SumPoints+Y(k+w); end s(k+halfw)=sum(Y(L-w+1:L)); SmoothY=s./w; % Taper the ends of the signal if ends=1. if ends==1, startpoint=(smoothwidth + 1)/2; SmoothY(1)=(Y(1)+Y(2))./2; for k=2:startpoint, SmoothY(k)=mean(Y(1:(2*k-1))); SmoothY(L-k+1)=mean(Y(L-2*k+2:L)); end SmoothY(L)=(Y(L)+Y(L-1))./2; end % ---------------------------------------------------------------------- function [FitResults,GOF,baseline,coeff,BestStart,xi,yi,BootResults]=peakfit(signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,autozero,fixedparameters,plots,bipolar,minwidth,DELTA,clipheight) % A command-line peak fitting program for time-series signals. % Version 7: March, 2015, adds peak shapes with three unconstrained % iterated variables: 30=voigt (variable alpha), 31=ExpGaussian (variable % time constant), 32=Pearson (variable shape factor), 34=Gaussian/ % Lorentzian blend (variable percent). % global AA xxx PEAKHEIGHTS FIXEDPARAMETERS AUTOZERO delta BIPOLAR CLIPHEIGHT format short g format compact warning off all NumArgOut=nargout; datasize=size(signal); if datasize(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 % Y=Y-min(Y); 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'; coeff=0; CLIPHEIGHT=max(Y); LOGPLOT=0; % Define values of any missing arguments % (signal,center,window,NumPeaks,peakshape,extra,NumTrials,start,autozero,fixedparameters,plots,bipolar,minwidth,DELTA) switch nargin case 1 NumPeaks=1; peakshape=1; extra=0; NumTrials=1; xx=X;yy=Y; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; plots=1; BIPOLAR=0; MINWIDTH=xx(2)-xx(1); delta=1; CLIPHEIGHT=max(Y); case 2 NumPeaks=1; peakshape=1; extra=0; NumTrials=1; xx=signal;yy=center; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; plots=1; BIPOLAR=0; MINWIDTH=xx(2)-xx(1); delta=1; CLIPHEIGHT=max(Y); case 3 NumPeaks=1; peakshape=1; extra=0; NumTrials=1; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=xx(2)-xx(1); delta=1; CLIPHEIGHT=max(Y); case 4 % Numpeaks specified in arguments peakshape=1; extra=0; NumTrials=1; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=xx(2)-xx(1); delta=1; CLIPHEIGHT=max(Y); case 5 % Numpeaks, peakshape specified in arguments extra=zeros(1,NumPeaks); NumTrials=1; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; CLIPHEIGHT=max(Y); case 6 NumTrials=1; start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; case 7 start=calcstart(xx,NumPeaks,xoffset); AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; CLIPHEIGHT=max(Y); case 8 AUTOZERO=0; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; CLIPHEIGHT=max(Y); case 9 AUTOZERO=autozero; FIXEDPARAMETERS=0; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; case 10 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; plots=1; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; case 11 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; BIPOLAR=0; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; CLIPHEIGHT=max(Y); case 12 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; BIPOLAR=bipolar; MINWIDTH=zeros(size(peakshape))+(xx(2)-xx(1)); delta=1; CLIPHEIGHT=max(Y); case 13 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; BIPOLAR=bipolar; MINWIDTH=minwidth; delta=1; case 14 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; BIPOLAR=bipolar; MINWIDTH=minwidth; delta=DELTA; CLIPHEIGHT=max(Y); case 15 AUTOZERO=autozero; FIXEDPARAMETERS=fixedparameters; BIPOLAR=bipolar; MINWIDTH=minwidth; delta=DELTA; CLIPHEIGHT=clipheight; otherwise end % switch nargin % Saturation Code, skips points greater than set maximum if CLIPHEIGHT3,AUTOZERO=3,end if AUTOZERO<0,AUTOZERO=0,end Heights=zeros(1,NumPeaks); FitResults=zeros(NumPeaks,6); % % Remove linear baseline from data segment if AUTOZERO==1 baseline=0; bkgcoef=0; 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; % Assign ShapStrings switch peakshape(1) 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='Up Sigmoid (logistic function)'; case 23 ShapeString='Down Sigmoid (logistic function)'; 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='Breit-Wigner-Fano'; case 16 ShapeString='Fixed-position Gaussians'; case 17 ShapeString='Fixed-position Lorentzians'; case 18 ShapeString='Exp. Lorentzian'; case 19 ShapeString='Alpha function'; case 20 ShapeString='Voigt (equal alphas)'; case 21 ShapeString='triangular'; case 22 ShapeString=num2str(shapesvector); case 24 ShapeString='Negative Binomial Distribution'; case 25 ShapeString='Lognormal Distribution'; case 26 ShapeString='slope'; case 27 ShapeString='First derivative'; case 28 ShapeString='Polynomial'; case 29 ShapeString='Segmented linear'; case 30 ShapeString='Voigt (variable alphas)'; case 31 ShapeString='ExpGaussian (var. time constant)'; case 32 ShapeString='Pearson (var. shape constant)'; case 33 ShapeString='Variable Gaussian/Lorentzian'; otherwise end % switch peakshape % Perform peak fitting for selected peak shape using fminsearch function options = optimset('TolX',.001,'Display','off','MaxFunEvals',1000 ); 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, % StartMatrix(k,:)=newstart; % disp(['Trial number ' num2str(k) ] ) % optionally prints the current trial number as progress indicator switch peakshape(1) case 1 TrialParameters=fminsearch(@(lambda)(fitgaussian(lambda,xx,yy)),newstart,options); for Peak=1:NumPeaks; if TrialParameters(2*Peak)=-BIPOLAR*10^100, % 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=Heights; % 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 MeanFitError1, title(['Best of ' num2str(NumTrials) ' fits']) else title(['Single fit']) end end % if plots % Put results into a matrix FitResults, one row for each peak, showing peak index number, % position, amplitude, and width. FitResults=zeros(NumPeaks,6); switch peakshape(1), case {6,7,8}, % equal-width peak models only for m=1:NumPeaks, if m==1, FitResults=[[round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]]; else FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) abs(FitParameters(NumPeaks+1)) area(m)]]; end end case {11,12}, % Fixed-width shapes only for m=1:NumPeaks, if m==1, FitResults=[[round(m) FitParameters(m)+xoffset height(m) FIXEDPARAMETERS area(m)]]; else FitResults=[FitResults ; [round(m) FitParameters(m)+xoffset height(m) FIXEDPARAMETERS area(m)]]; end end case {16,17}, % Fixed-position shapes only for m=1:NumPeaks, if m==1, FitResults=[round(m) FIXEDPOSITIONS(m) height(m) FitParameters(m) area(m)]; else FitResults=[FitResults ; [round(m) FIXEDPOSITIONS(m) height(m) FitParameters(m) area(m)]]; end end case 28, % Simple polynomial fit FitResults=PEAKHEIGHTS; case 29, % Segmented linear fit FitResults=PEAKHEIGHTS; case {30,31,32,33} % Special case of shapes with 3 iterated variables for m=1:NumPeaks, if m==1, FitResults=[round(m) FitParameters(3*m-2) height(m) abs(FitParameters(3*m-1)) area(m) FitParameters(3*m)]; else FitResults=[FitResults ; [round(m) FitParameters(3*m-2) height(m) abs(FitParameters(3*m-1)) area(m)] FitParameters(3*m)]; end end otherwise % Normal shapes with 2 iterated variables for m=1:NumPeaks, if m==1, FitResults=[round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]; else FitResults=[FitResults ; [round(m) FitParameters(2*m-1)+xoffset height(m) abs(FitParameters(2*m)) area(m)]]; end % if m==1 end % for m=1:NumPeaks, end % switch peakshape(1) % Display Fit Results on lower graph if plots, % Display Fit Results on lower graph subplot(2,1,2); startx=min(xx)+(max(xx)-min(xx))./20; dxx=(max(xx)-min(xx))./10; dyy=((max(residual)-min(residual))./10); starty=max(residual)-dyy; FigureSize=get(gcf,'Position'); switch peakshape(1) case {9,19,11} % Pulse and sigmoid shapes only text(startx,starty+dyy/2,['Peak # tau1 Height tau2 Area'] ); case 28, % Polynomial text(startx,starty+dyy/2,['Polynomial coefficients'] ); case 29 % Segmented linear text(startx,starty+dyy/2,['x-axis breakpoints'] ); case {30,31,32,33} % Special case of shapes with 3 iterated variables text(startx,starty+dyy/2,['Peak # Position Height Width Area Shape factor'] ); otherwise text(startx,starty+dyy/2,['Peak # Position Height Width Area '] ); end % Display FitResults using sprintf if peakshape(1)==28||peakshape(1)==29, % Polynomial or segmented linear for number=1:length(FitResults), column=1; itemstring=sprintf('%0.4g',FitResults(number)); xposition=startx+(1.7.*dxx.*(column-1).*(600./FigureSize(3))); yposition=starty-number.*dyy.*(400./FigureSize(4)); text(xposition,yposition,[' ' itemstring]); end else 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 xposition=startx; yposition=starty-(peaknumber+1).*dyy.*(400./FigureSize(4)); if AUTOZERO==3, text(xposition,yposition,[ 'Baseline= ' num2str(baseline) ]); end % if AUTOZERO end % if peakshape(1) if peakshape(1)==30 || peakshape(1)==31 || peakshape(1)==32 || peakshape(1)==33, for peaknumber=1:NumPeaks, column=6; itemstring=sprintf('%0.4g',FitParameters(3*peaknumber)); 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==8, if plots,disp('Computing bootstrap sampling statistics.....'),end BootstrapResultsMatrix=zeros(6,100,NumPeaks); BootstrapErrorMatrix=zeros(1,100,NumPeaks); clear bx by tic; for trial=1:100, n=1; bx=xx; by=yy; while n.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,FIXEDPARAMETERS); for peak=1:NumPeaks, switch peakshape(1) case {30,31,32,33} BootstrapResultsMatrix(1:6,trial,peak)=FitResults(peak,1:6); otherwise BootstrapResultsMatrix(1:5,trial,peak)=FitResults(peak,1:5); end 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 Shape Factor']); 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:6); BootstrapSTD=BootstrapSTD(2:6); BootstrapIQR=BootstrapIQR(2:6); PercentRSD=PercentRSD(2:6); PercentIQR=PercentIQR(2:6); 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==8, if AUTOZERO==3; else baseline=bkgcoef; end % ---------------------------------------------------------------------- function [FitResults,LowestError]=fitpeaks(xx,yy,NumPeaks,peakshape,extra,NumTrials,start,AUTOZERO,fixedparameters) % Based on peakfit Version 3: June, 2012. global PEAKHEIGHTS FIXEDPARAMETERS AUTOZERO BIPOLAR MINWIDTH coeff format short g format compact warning off all FIXEDPARAMETERS=fixedparameters; xoffset=0; if start==0;start=calcstart(xx,NumPeaks,xoffset);end PEAKHEIGHTS=zeros(1,NumPeaks); n=length(xx); newstart=start; coeff=0; LOGPLOT=0; % Perform peak fitting for selected peak shape using fminsearch function options = optimset('TolX',.001,'Display','off','MaxFunEvals',1000 ); 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, % StartVector=newstart switch peakshape(1) case 1 TrialParameters=fminsearch(@(lambda)(fitgaussian(lambda,xx,yy)),newstart,options); for Peak=1:NumPeaks; if TrialParameters(2*Peak)=-BIPOLAR*10^100, % 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=Heights; % Assign height to the PEAKHEIGHTS with the lowest MeanFitError end % if min(PEAKHEIGHTS)>0 end % if MeanFitError0); g = p'; % ---------------------------------------------------------------------- function err = fitalphafunction(tau,x,y) % Iterative fit of the sum of alpha funciton % of the form Height.*exp(-tau1.*x).*(1-exp(-tau2.*x))) global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = alphafunction(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function g = alphafunction(x,pos,spoint) % alpha function. pos=position; wid=half-width (both scalar) % alphafunction(x,pos,wid), where x may be scalar, vector, or matrix % pos=position; wid=half-width (both scalar) % Taekyung Kwon, July 2013 g = (x-spoint)./pos.*exp(1-(x-spoint)./pos); for m=1:length(x);if g(m)<0;g(m)=0;end;end % ---------------------------------------------------------------------- function err = fitdownsigmoid(tau,x,y) % Fitting function for iterative fit to the sum of % downward moving sigmiods global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = downsigmoid(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function err = fitupsigmoid(tau,x,y) % Fitting function for iterative fit to the sum of % upwards moving sigmiods global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = upsigmoid(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function g=downsigmoid(x,t1,t2) % down step sigmoid g=.5-.5*erf(real((x-t1)/sqrt(2*t2))); % ---------------------------------------------------------------------- function g=upsigmoid(x,t1,t2) % up step sigmoid g=1/2 + 1/2* erf(real((x-t1)/sqrt(2*t2))); % ---------------------------------------------------------------------- function err = fitGL(lambda,t,y,shapeconstant) % Fitting functions for Gaussian/Lorentzian blend. % T. C. O'Haver (toh@umd.edu), 2012. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT 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 if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function err = fitGLv(lambda,t,y) % Fitting functions for Gaussian/Lorentzian blend function with % independently variable percent Gaussian % T. C. O'Haver (toh@umd.edu), 2015. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(t),round(length(lambda)/3)); for j = 1:length(lambda)/3, A(:,j) = GL(t,lambda(3*j-2),lambda(3*j-1),lambda(3*j))'; end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- 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 % sizex=size(x) % sizepos=size(pos) % sizewid=size(wid) % sizem=size(m) g=2.*((m/100).*gaussian(x,pos,wid)+(1-(m(1)/100)).*lorentzian(x,pos,wid))/2; % ---------------------------------------------------------------------- function err = fitvoigt(lambda,t,y,shapeconstant) % Fitting functions for Voigt profile function % T. C. O'Haver (toh@umd.edu), 2013. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(t),round(length(lambda)/2)); for j = 1:length(lambda)/2, A(:,j) = voigt(t,lambda(2*j-1),lambda(2*j),shapeconstant)'; end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function err = fitvoigtv(lambda,t,y) % Fitting functions for Voigt profile function with independently variable % alphas % T. C. O'Haver (toh@umd.edu), 2015. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(t),round(length(lambda)/3)); for j = 1:length(lambda)/3, A(:,j) = voigt(t,lambda(3*j-2),lambda(3*j-1),lambda(3*j))'; end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function g=voigt(xx,pos,gD,alpha) % Voigt profile function. xx is the independent variable (energy, % wavelength, etc), gD is the Doppler (Gaussian) width, and alpha is the % shape constant (ratio of the Lorentzian width gL to the Doppler width gD. % Based on Chong Tao's "Voigt lineshape spectrum simulation", % File ID: #26707 % alpha=alpha gL=alpha.*gD; gV = 0.5346*gL + sqrt(0.2166*gL.^2 + gD.^2); x = gL/gV; y = abs(xx-pos)/gV; g = 1/(2*gV*(1.065 + 0.447*x + 0.058*x^2))*((1-x)*exp(-0.693.*y.^2) + (x./(1+y.^2)) + 0.016*(1-x)*x*(exp(-0.0841.*y.^2.25)-1./(1 + 0.021.*y.^2.25))); g=g./max(g); % ---------------------------------------------------------------------- function err = fitBiGaussian(lambda,t,y,shapeconstant) % Fitting functions for BiGaussian. % T. C. O'Haver (toh@umd.edu), 2012. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT 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 if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- 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 = fitBWF(lambda,t,y,shapeconstant) % Fitting function for Breit-Wigner-Fano. % T. C. O'Haver (toh@umd.edu), 2014. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(t),round(length(lambda)/2)); for j = 1:length(lambda)/2, A(:,j) = BWF(t,lambda(2*j-1),lambda(2*j),shapeconstant)'; end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function g = BWF(x,pos,wid,m) % BWF (Breit-Wigner-Fano) http://en.wikipedia.org/wiki/Fano_resonance % pos=position; wid=width; m=Fano factor % T. C. O'Haver, 2014 y=((m*wid/2+x-pos).^2)./(((wid/2).^2)+(x-pos).^2); % y=((1+(x-pos./(m.*wid))).^2)./(1+((x-pos)./wid).^2); g=y./max(y); % ---------------------------------------------------------------------- function err = fitnbinpdf(tau,x,y) % Fitting function for iterative fit to the sum of % Negative Binomial Distributions % (http://www.mathworks.com/help/stats/negative-binomial-distribution.html) global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = nbinpdf(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function err = fitlognpdf(tau,x,y) % Fitting function for iterative fit to the sum of % Lognormal Distributions % (http://www.mathworks.com/help/stats/lognormal-distribution.html) global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = lognormal(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- 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 = fitsine(tau,x,y) % Fitting function for iterative fit to the sum of % sine waves (alpha test, NRFPT) global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, A(:,j) = sine(x,tau(2*j-1),tau(2*j)); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function g=sine(x,f,phase) % Sine wave (alpha test) g=sin(2*pi*f*(x+phase)); % ---------------------------------------------------------------------- function err = fitd1gauss(lambda,t,y) % Fitting functions for the first derivative of a Gaussian % T. C. O'Haver, 2014 global PEAKHEIGHTS AUTOZERO BIPOLAR A = zeros(length(t),round(length(lambda)/2)); for j = 1:length(lambda)/2, A(:,j) = d1gauss(t,lambda(2*j-1),lambda(2*j))'; end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; err = norm(z-y'); % ---------------------------------------------------------------------- function y=d1gauss(x,p,w) % First derivative of Gaussian (alpha test) y=-(5.54518.*(x-p).*exp(-(2.77259.*(p-x).^2)./w^2))./w.^2; y=y./max(y); % ---------------------------------------------------------------------- function coeff = fitpolynomial(t,y,order) coeff=polyfit(t,y,order); % order=order % coeff=coeff % ---------------------------------------------------------------------- function y=polynomial(t,coeff) y=polyval(coeff,t); % ---------------------------------------------------------------------- function err = fitsegmented(lambda,t,y) % Fitting functions for articulated segmented linear % T. C. O'Haver, 2014 global LOGPLOT breakpoints=[t(1) lambda max(t)]; z = segmented(t,y,breakpoints); % lengthz=length(z); if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y); end % ---------------------------------------------------------------------- function yi=segmented(x,y,segs) global PEAKHEIGHTS clear yy for n=1:length(segs) yind=val2ind(x,segs(n)); yy(n)=y(yind(1)); end yi=INTERP1(segs,yy,x); PEAKHEIGHTS=segs; % ---------------------------------------------------------------------- function err = fitlinslope(tau,x,y) % Fitting function for iterative fit to linear function global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT A = zeros(length(x),round(length(tau)/2)); for j = 1:length(tau)/2, z = (x.*tau(2*j-1)+tau(2*j))'; A(:,j) = z./max(z); end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function y=linslope(x,slope,intercept) y=x.*slope+intercept; % y=y./max(y); % ---------------------------------------------------------------------- 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,:)); % ---------------------------------------------------------------------- function err = fitmultiple(lambda,t,y,shapesvector,m) % Fitting function for a multiple-shape band signal. % The sequence of peak shapes are defined by the vector "shape". % The vector "m" determines the shape of variable-shape peaks. global PEAKHEIGHTS AUTOZERO BIPOLAR LOGPLOT coeff numpeaks=round(length(lambda)/2); A = zeros(length(t),numpeaks); for j = 1:numpeaks, if shapesvector(j)==28, coeff=polyfit(t,y,m(j)); A(:,j) = polyval(coeff,t); else A(:,j) = peakfunction(shapesvector(j),t,lambda(2*j-1),lambda(2*j),m(j))'; end end if AUTOZERO==3,A=[ones(size(y))' A];end if BIPOLAR,PEAKHEIGHTS=A\y';else PEAKHEIGHTS=abs(A\y');end z = A*PEAKHEIGHTS; if LOGPLOT, err = norm(log10(z)-log10(y)'); else err = norm(z-y'); end % ---------------------------------------------------------------------- function p=peakfunction(shape,x,pos,wid,m,coeff) % function that generates any of 20 peak types specified by number. 'shape' % specifies the shape type of each peak in the signal: "peakshape" = 1-20. % 1=Gaussian 2=Lorentzian, 3=logistic, 4=Pearson, 5=exponentionally % broadened Gaussian; 9=exponential pulse, 10=up sigmoid, % 13=Gaussian/Lorentzian blend; 14=BiGaussian, 15=Breit-Wigner-Fano (BWF) , % 18=exponentionally broadened Lorentzian; 19=alpha function; 20=Voigt % profile; 21=triangular; 23=down sigmoid; 25=lognormal. "m" is required % for variable-shape peaks only. switch shape, case 1 p=gaussian(x,pos,wid); case 2 p=lorentzian(x,pos,wid); case 3 p=logistic(x,pos,wid); case 4 p=pearson(x,pos,wid,m); case 5 p=expgaussian(x,pos,wid,m); case 6 p=gaussian(x,pos,wid); case 7 p=lorentzian(x,pos,wid); case 8 p=expgaussian(x,pos,wid,m)'; case 9 p=exppulse(x,pos,wid); case 10 p=upsigmoid(x,pos,wid); case 11 p=gaussian(x,pos,wid); case 12 p=lorentzian(x,pos,wid); case 13 p=GL(x,pos,wid,m); case 14 p=BiGaussian(x,pos,wid,m); case 15 p=BWF(x,pos,wid,m); case 16 p=gaussian(x,pos,wid); case 17 p=lorentzian(x,pos,wid); case 18 p=explorentzian(x,pos,wid,m)'; case 19 p=alphafunction(x,pos,wid); case 20 p=voigt(x,pos,wid,m); case 21 p=triangular(x,pos,wid); case 23 p=downsigmoid(x,pos,wid); case 25 p=lognormal(x,pos,wid); case 26 p=linslope(x,pos,wid); case 27 p=d1gauss(x,pos,wid); case 28 p=polynomial(x,coeff); otherwise end % switch