% Script to test segmented peak finder findpeaksSG. Requires findpeaksSG.m % and gaussian.m (or other peak shape) function. Creates 100 noisy signals % with 3 peaks with different widths, then measures the peak positions % heights and widths, and computes the percent relative standard deviation % of 100 measurements. The segmented peak finder findpeaksSG, with % 3-segment peak detection parameters, reliably detects amd measures all % three peaks. In contrast, findpeaksG tuned to the middle peak (using line % 26 instead of line 25) measures the first and last peaks poorly. You can % change any of the variables in lines 10-18. clear clf p=[200 400 800]; % Peak positions of the 3 peaks. h=[1 1 1]; % Peak heights of the 3 peaks. w=[10 50 200]; % Peak widths of the 3 peaks. Noise=.04; % Standard deviation of the random noise % Define the peak detection parameters for each segment SlopeThresholds=[.001 .00001 .000001]; AmpThresholds=[.1 .1 .1]; smoothwidths=[10 50 80]; peakgroups=[8 50 250]; NumTrials=100; % Number of independent signals measured x=1:1000; for trial=1:NumTrials, % Create the noisy signal y=h(1).*gaussian(x,p(1),w(1))+h(2).*gaussian(x,p(2),w(2))+h(3).*gaussian(x,p(3),w(3))+Noise.*randn(size(x)); % Detect and measure peaks in the signal P=findpeaksSG(x,y,SlopeThresholds,AmpThresholds,smoothwidths,peakgroups,3); % P=findpeaksG(x,y,SlopeThresholds(2),AmpThresholds(2),smoothwidths(2),peakgroups(2),3); sizeP=size(P); NumPeaks(trial)=sizeP(1); if trial==1; % Plot the first signal plot(x,y,'c'), end hold on for peak=1:NumPeaks(trial), text(P(peak,2),P(peak,3),num2str(peak)) % Number the peaks % Record the measured peak parameters of each peak in each trial PeakPositions(trial,peak)=P(peak,2); PeakHeights(trial,peak)=P(peak,3); PeakWidths(trial,peak)=P(peak,4); PeakAreas(trial,peak)=P(peak,5); end end hold off xlabel('Each detected peak position and height is marked with peak number 1, 2 or 3') ylabel('Y') title('TestPrecisionFindpeaskSG.m: 100 repeat measurements of noisy signal with 3 peaks') disp(' ') disp(['Average number of peaks: ' num2str(mean(NumPeaks)) ]) disp('Percent Relative Standard Deviations (RSD):') disp(' Peak 1 Peak 2 Peak 3 Peak 4') RSDPeakPositions=100.*std(PeakPositions)./mean(PeakPositions) RSDPeakHeights=100.*std(PeakHeights)./mean(PeakHeights) RSDPeakWidths=100.*std(PeakWidths)./mean(PeakWidths) RSDPeakAreas=100.*std(PeakAreas)./mean(PeakAreas)