% Demonstration of findpeaks2d.m, which shows that this function can locate % peaks resulting in 'shoulders' that do not produce a distinct maximum in % the original signal. Detects peaks in the negative of the smoothed second % derivative of the signal (shown as the dotted line in the figure). % Requires gaussian.m, findpeaksG.m, findpeaksG2d.m, fastsmooth.m, and % peakfit.m in the path. Download from http://tinyurl.com/cey8rwh format compact figure(1); x=1:.1:50; noise=.1; % << change: Try larger or smaller amounts of noise AmpThreshold=0.4; SlopeThreshold=1e-7; smoothwidth=25; y1=gaussian(x,20,10); % Peak 1: Position=20 Height=1 Width=10 y2=.5.*gaussian(x,30,10); % Peak 2: Position=30 Height=0.5 Width=10 n=noise.*randn(size(x)); % random white noise y=y1+y2+n; % Simulated signal contains two peaks plus noise clf; sy=fastsmooth(y,smoothwidth,3); % Pre-smoothing before calling findpeaks ssd=-2000.*fastsmooth(deriv2(sy),smoothwidth,3); % -Smoothed second derivative plot(x,y1,x,y2,x,y,x,ssd,'-.'); disp(' ') noise=noise disp(' Peak Position Height Width Area') PfindpeaksG=findpeaksG(x,sy,SlopeThreshold,AmpThreshold,smoothwidth,25,3) PfindpeaksG2d=findpeaksG2d(x,sy,SlopeThreshold,AmpThreshold,smoothwidth,25,3) xlabel('Blue= Peak 1 Green= Peak 2 Red= Peak 1 + Peak 2 Dotted= -smoothed second derivative'); title('Comparison of findpeaksG.m and findpeaksG2d.m for locating a weak side shoulder') disp('findpeaksG detects only one peak, but findpeaksG2d detects both. ') disp(' ') disp('If the findpeaksG2d results are not accurate enough, you can use its ') disp('results as the "start" value for peakfit.m, which takes longer but gives') disp('more accurate results, especially for width and area:') % Comparison to peakfit sizeP=size(PfindpeaksG2d); NumPeaks=sizeP(1); start=[PfindpeaksG2d(1,2) PfindpeaksG2d(1,4) PfindpeaksG2d(2,2) PfindpeaksG2d(2,4)]; figure(2) disp('PeakfitResults=peakfit([x;y],0,0,NumPeaks,0,10,start)') PeakfitResults=peakfit([x;y],0,0,NumPeaks,0,10,start)