% 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 11-24. T. C. O'Haver % (toh@umd.edu) Ver. 7: November, 2017. format short 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.5 6.5]; % Positions of the peaks (in increasing order) wid=[1.666 1.666]; % Widths of the peaks (full width at half-maximum) noise=.02; % Standard deviation of random noise added to signal NumRepeats=100; % Number of repeat samples with independent noise NumTrials=3; % Number of trail fits per sample PeakShape=1; % Gaussian=1, Lorentzian=2, etc. extra=2; center=5; % x-value in center of fitted range window=10; % range of x-values in fitted range startvector=[4 1 7 1]; % First guess vector [pos1 wid1 pos2 wid2 ...] Must be two values for each peak baselinemode=0; % 0=no baseline correction; 1=linear baseline correction % 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 % Generate and fit NumRepeats separate samples of the simulated signal with noise. NumPeaks=length(amp); AllFitResults=zeros(NumPeaks,5,NumRepeats); AllFittingErrors=zeros(1,NumRepeats); AllPercentParameterErrors=zeros(NumPeaks,5,NumRepeats); for repeat=1:NumRepeats % Add different white noise sample to each repeat of simulated signal y = TrueY + noise .* randn(size(x)); % Execute the peakfit function [FitResults,LowestError]=peakfit([x' y'],center,window,NumPeaks,PeakShape,extra,NumTrials,startvector,baselinemode); FitResults=sortrows(FitResults,2); % Sort results by increasing peak positions (column 2) Heights1(repeat)=FitResults(1,3); Positions1(repeat)=FitResults(1,2); Widths1(repeat)=FitResults(1,4); Areas1(repeat)=FitResults(1,5); Heights2(repeat)=FitResults(2,3); Positions2(repeat)=FitResults(2,2); Widths2(repeat)=FitResults(2,4); Areas2(repeat)=FitResults(2,5); end % Calculate the mean and relative standard deviaiton of peak parameters % measured by curve fitting. disp('Mean Fit Results') disp(' Position Height Width Area') disp(['Peak 1 ' num2str([mean(Positions1) mean(Heights1) mean(Widths1) mean(Areas1)]) ] ) disp(['Peak 2 ' num2str([mean(Positions2) mean(Heights2) mean(Widths2) mean(Areas2)]) ] ) disp(' ') disp('Percent relative standard deviations') disp(' Position Height Width Area') disp(['Peak 1 ' num2str([rsd(Positions1).*100 rsd(Heights1).*100 rsd(Widths1).*100 rsd(Areas1).*100]) ] ) disp(['Peak 2 ' num2str([rsd(Positions2).*100 rsd(Heights2).*100 rsd(Widths2).*100 rsd(Areas2).*100]) ] ) % Required functions: peakfit.m, rsd.m (shown below) % function relstddev=rsd(x) % % Relative standard deviation of vector x % relstddev=std(x)./mean(x);