% Demonstration script for peakfit.m. Generates an overlapping peak signal, % adds noise, fits it with peakfit.m (in line 73), then computes the % accuracy and precision of peak paramater measurements. You can change % any of the initial values in lines 13-29. % T. C. O'Haver (toh@umd.edu). Version 5: March, 2011. clear close format long g format compact warning off all % Initial settings of users-adjustable parameters n=300; % maximum x-value increment=1; % Difference between adjacent x values x=[1:increment:n]; % Creates x-vector amp=[1 2 3]; % Amplitudes of the peaks (one entry for each peak) pos=[110 140 170]; % Positions of the peaks (in increasing order) wid=[30 30 30]; % Widths of the peaks (full width at half-maximum) noise=.05; % Amplitude of random noise added to signal NumRepeats=10; % Number of repeat samples with independent noise NumTrials=5; % Number of trail fits per sample PeakShape=1; % Gaussian=1, Lorentzian=2, Logistic=3, Pearson=4 % "extra" determines the shape of the Pearson function. When extra=1 % the shape is Lorentzian; when extra is large (>20), the shape % approaches Gaussian. At small values (<1), the shape is a pointy "cusp". extra=2; center=150; % x-value in center of fitted range window=200; % range of x-values in fitted range startvector=[100 40 150 40 200 40]; % 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; 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; AllPercentParameterErrors(:,:,repeat)=100.*(FitResults-ActualParameters)./ActualParameters; end % Calculate the mean and relative standard deviaiton of peak parameters % measured by curve fitting. for peak=1:NumPeaks, MeanParameters(peak,:)=mean(squeeze(AllFitResults(peak,:,:))'); RSDparameters(peak,:)=100.*std(squeeze(AllFitResults(peak,:,:))')./mean(squeeze(AllFitResults(peak,:,:))'); end disp('-------------------------------------------------------') 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('Actual parameters:'); disp('Position Height Width Area'); disp(num2str(ActualParameters(:,2:5))) disp(['Average percent fitting error: ' num2str(mean(AllFittingErrors))]) disp('Average measured parameters:'); disp('Position Height Width Area'); disp(num2str(MeanParameters(:,2:5))) RSDparameters(:,1)=ActualParameters(:,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('Average %RSD of measured parameters: '); disp(num2str(meanRSD)) PercentParameterErrors=100.*(MeanParameters-ActualParameters)./ActualParameters; PercentParameterErrors(:,1)=ActualParameters(:,1); disp('Percent errors of measured parameters:'); disp('Position Height Width Area'); disp(num2str(PercentParameterErrors(:,2:5))) AverageParameterErrors=mean(abs(PercentParameterErrors(:,2:5))); disp('Average Percent Parameter Error: '); disp(num2str(AverageParameterErrors))