% Demonstration script (for peakfit.m, version 3.1). Generates an % overlapping peak signal, adds noise, fits it with peakfit.m , then % computes the accuracy and precision of peak paramater measurements. % You can change any of the initial values in lines 13-30. T. C. O'Haver % (toh@umd.edu) Ver. 6: April, 2011, for use with peakfit ver. 3.1 (2012) format long g format compact warning off all disp('Computing......... may take several seconds.') % >>> Initial settings of user-adjustable parameters <<< n=10; % maximum x-value increment=.1; % Difference between adjacent x values x=[1:increment:n]; % Creates x-vector amp=[1 .5]; % Amplitude s of the peaks (one entry for each peak) pos=[4 6]; % Positions of the peaks (in increasing order) wid=[1.666 1.666]; % Widths of the peaks (full width at half-maximum) noise=.05; % Amplitude of random noise added to signal NumRepeats=100; % Number of repeat samples with independent noise NumTrials=10; % Number of trail fits per sample PeakShape=1; % Gaussian=1, Lorentzian=2, Logistic=3, Pearson=4 % ("extra" determines the shape of the Pearson function: as extra % goes from 1 -> large, shape goes from Lorentzian -> Gaussian). % Exponentially-broadened Gaussian = 5; Equal-widths Gaussians = 6; % Equal-widths Lorentzians = 7. extra=2; center=5; % x-value in center of fitted range window=10; % range of x-values in fitted range startvector=[5 2 5 2]; % First guess vector [pos1 wid1 pos2 wid2 ...] Must be two values for each peak % Generate simulated signal peaks % 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), switch PeakShape case 1 A(k,:)=exp(-((x-pos(k))./(0.6006.*wid(k))) .^2); % Gaussian function case 2 A(k,:)=ones(size(x))./(1+((x-pos(k))./(0.5.*wid(k))).^2); % Lorentzian function case 3 L = exp(-((x-pos(k))/(.477.*wid(k))) .^2); A(k,:)= (2.*L)./(1+L); case 4 % pearson function (shape determined by value of extra in line 21 A(k,:)=ones(size(x))./(1+((x-pos(k))./((0.5.^(2/extra)).*4.62.*wid(k))).^2).^extra; case 6 A(k,:)=exp(-((x-pos(k))./(0.6006.*wid(k))) .^2); % Gaussian function case 7 A(k,:)=ones(size(x))./(1+((x-pos(k))./(0.5.*wid(k))).^2); % Laorentzian function end % switch area(k)=amp(k).*trapz(x,A(k,:)); % Computes peak areas via trapaziod method end % for TrueY=amp*A; % Multiplies each row by the corresponding amplitude and adds them up % To test the autozero baseline subtraction, un-comment next line. % TrueY=TrueY+ones(size(x))./(1+((x+100)./(300)).^2); % Adds sloping curved baseline % 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 % Generate and fit NumRepeats separate samples of the simulated signal with noise. AllFitResults=zeros(NumPeaks,5,NumRepeats); AllFittingErrors=zeros(1,NumRepeats); AllPercentParameterErrors=zeros(NumPeaks,5,NumRepeats); for repeat=1:NumRepeats, % Add normally-distributed random noise to simulated signal y = TrueY + noise .* randn(size(x)); % Execute the peakfit script [FitResults,LowestError]=peakfit([x' y'],center,window,NumPeaks,PeakShape,extra,NumTrials,startvector); FitResults=SORTROWS(FitResults,2); % Sort results by increasing peak positions (column 2) AllFitResults(:,:,repeat)=FitResults; AllFittingErrors(repeat)=LowestError(1); AllPercentParameterErrors(:,:,repeat)=100.*(FitResults-ActualParameters)./ActualParameters; end % Calculate the mean and relative standard deviaiton of peak parameters % measured by curve fitting. \ if NumRepeats>1, for peak=1:NumPeaks, MeanParameters(peak,:)=mean(squeeze(AllFitResults(peak,:,:))'); RSDparameters(peak,:)=100.*std(squeeze(AllFitResults(peak,:,:))')./mean(squeeze(AllFitResults(peak,:,:))'); end else MeanParameters=AllFitResults; end disp('-------------------------------------------------------') switch PeakShape case 1 ShapeString='1 Gaussian'; case 2 ShapeString='2 Lorentzian'; case 3 ShapeString='3 Logistic'; case 4 ShapeString='4 Pearson'; case 5 ShapeString='5 ExpGaussian'; case 6 ShapeString='6 Equal width Gaussians'; case 7 ShapeString='7 Equal width Lorentzians'; otherwise end % switch peakshape disp([ 'Peak shape: ' ShapeString ]) disp(['Number of peaks: ' num2str(NumPeaks)]) disp(['Noise amplitude: ' num2str(noise)]) disp(['Number of data points in signal: ' num2str(length(x))]) disp(['Number of repeat signals: ' num2str(NumRepeats)]) disp(['Number of trial fits per sample signal: ' num2str(NumTrials)]) disp(' ') disp('Actual parameters:'); disp('Position Height Width Area'); disp(num2str(ActualParameters(:,2:5))) disp(' ') disp(['Average percent fitting error: ' num2str(mean(AllFittingErrors))]) disp('Average measured parameters:'); disp('Position Height Width Area'); if NumRepeats>1, disp(num2str(MeanParameters(:,2:5))) else disp(num2str(AllFitResults(:,2:5))) end % NumRepeats>1, disp(' ') RSDparameters(:,1)=ActualParameters(:,1); % Replace first column with peak numbers if NumRepeats>1, disp('Percent RSDs of measured parameters:'); disp('Position Height Width Area'); disp(num2str(RSDparameters(:,2:5))) meanRSD=mean(RSDparameters); meanRSD=meanRSD(2:5); disp(' ') disp('Average %RSD of measured parameters for all peaks: '); disp(num2str(meanRSD)) end % if NumRepeats>1 if NumRepeats>1, PercentParameterErrors=100.*(MeanParameters-ActualParameters)./ActualParameters; else PercentParameterErrors=100.*(MeanParameters-ActualParameters)./ActualParameters; end % NumRepeats>1, PercentParameterErrors(:,1)=ActualParameters(:,1); % Replace first column with peak numbers disp(' ') disp('Percent errors of measured parameters:'); disp('Position Height Width Area'); disp(num2str(PercentParameterErrors(:,2:5))) AverageParameterErrors=mean(abs(PercentParameterErrors(:,2:5))); disp(' ') disp('Average Percent Parameter Error for all peaks: '); disp(num2str(AverageParameterErrors)) disp(' ') disp('If you have peakfit version 3.1 installed, you can compare ') disp('the standard deviations and percent errors of the measured ') disp('parameters (based on 100 repeat signals) to the bootstrap') disp('estimates derived from a single sample of the signal: ') disp('>> [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit([x'' y''],0,0,2);') [FitResults,LowestError,BestStart,xi,yi,BootResults]=peakfit([x' y'],0,0,2); disp(' ') FitResults disp('Almost always, the percent errors of the measured peak parameters ') disp('are within one standard deviation of the true values. ')