Code covered by the BSD License  

Highlights from
1D Non-derivative Peak Finder

image thumbnail

1D Non-derivative Peak Finder

by

 

21 Feb 2011 (Updated )

Up-sample and filter noisy data to find peaks without using derivatives.

varargout =peakfind(x_data,y_data,varargin)
function varargout =peakfind(x_data,y_data,varargin)
%
%PEAKFIND general 1D peak finding algorithm
%Tristan Ursell, 2013.
%
%peakfind(x_data,y_data)
%peakfind(x_data,y_data,upsam)
%peakfind(x_data,y_data,upsam,gsize,gstd)
%peakfind(x_data,y_data,upsam,htcut,'cuttype')
%peakfind(x_data,y_data,upsam,gsize,gstd,htcut,'cuttype')
%
%[xpeaks]=peakfind()
%[xout,yout,peakspos]=peakfind()
%
%This function finds peaks without taking first or second derivatives,
%rather it uses local slope features in a given data set.  The function has
%four basic modes.
%
%   Mode 1: peakfind(x_data,y_data) simply finds all peaks in the data
%   given by 'xdata' and 'ydata'.
%
%   Mode 2: peakfind(x_data,y_data,upsam) finds peaks after up-sampling
%   the data by the integer factor 'upsam' -- this allows for higher
%   resolution peak finding.  The interpolation uses a cubic spline that
%   does not introduce fictitious peaks.
%
%   Mode 3: peakfind(x_data,y_data,upsam,gsize,gstd) up-samples and then
%   convolves the data with a Gaussian point spread vector of length
%   gsize (>=3) and standard deviation gstd (>0).
%
%   Mode 4: peakfind(x_data,y_data,upsam,htcut,'cuttype') up-samples the
%   data (upsam=1 analyzes the data unmodified).  The string 'cuttype' can
%   either be 'abs' (absolute) or 'rel' (relative), which specifies a peak 
%   height cutoff which is either:
%
%       'abs' - finds peaks that lie an absolute amount 'htcut' above the
%       lowest value in the data set.
%
%            for (htcut > 0) peaks are found if
%
%            peakheights > min(yout) + htcut
%       
%       'rel' - finds peaks that are an amount 'htcut' above the lowest
%       value in the data set, relative to the full change in y-input
%       values.
%
%           for (0 < htcut < 1) peaks are found if
%
%           (peakheights-min(yout))/(max(yout)-min(yout)) > htcut
%
%Upsampling and convolution allows one to find significant peaks in noisy
%data with sub-pixel resolution.  The algorithm also finds peaks in data 
%where the peak is surrounded by zero first derivatives, i.e. the peak is 
%actually a large plateau.  
%
%The function outputs the x-position of the peaks in 'xpeaks' or the
%processed input data in 'xout' and 'yout' with 'peakspos' as the indices
%of the peaks, i.e. xpeaks = xout(peakspos).  
%
%If you want the algorithm to find the position of minima, simply input 
%'-y_data'.  Peaks within half the convolution box size of the boundary 
%will be ignored (to avoid this, pad the data before processing).
%
% Example 1:
%
% x_data = -50:50;
% y_data =(sin(x_data)+0.000001)./(x_data+0.000001)+1+0.025*(2*rand(1,length(x_data))-1);
%
% [xout,yout,peakspos]=peakfind(x_data,y_data);
%
% plot(x_data,y_data,'r','linewidth',2)
% hold on
% plot(xout,yout,'b','linewidth',2)
% plot(xout(peakspos),yout(peakspos),'g.','Markersize',30)
% xlabel('x')
% ylabel('y')
% title(['Found ' num2str(length(peakspos)) ' peaks.'])
% box on
%
% Example 2:
%
% x_data = -50:50;
% y_data =(sin(x_data)+0.000001)./(x_data+0.000001)+1+0.025*(2*rand(1,length(x_data))-1);
%
% [xout,yout,peakspos]=peakfind(x_data,y_data,4,6,2,0.2,'rel');
%
% plot(x_data,y_data,'r','linewidth',2)
% hold on
% plot(xout,yout,'b','linewidth',2)
% plot(xout(peakspos),yout(peakspos),'g.','Markersize',30)
% xlabel('x')
% ylabel('y')
% title(['Found ' num2str(length(peakspos)) ' peaks.'])
% box on


%check if input data are same length
if length(x_data)~=length(y_data)
    error('Input data vectors must be the same length.')
end

%check to make sure data is all real
if or(sum(isreal(x_data)==0)>0,sum(isreal(y_data)==0)>0)
    error('All input data must be real valued.')
end

%check if input data is flat (no peaks)
if sum(diff(y_data)==0)==(length(y_data)-1)
    %disp('Warning: Input data is completely flat.')
    if nargout==1
        varargout(1)={[]};
    elseif nargout==3
        varargout(3)={[]};
        varargout(2)={y_data};
        varargout(1)={x_data};
    else
        error('Incorrect number of output arguments.')
    end   
    return
end

%check if input data is monotonic (no peaks)
if or(sum(diff(y_data)<=0)==(length(y_data)-1),sum(diff(y_data)>=0)==(length(y_data)-1))
    %disp('Warning: Input data is monotonic.')
    if nargout==1
        varargout(1)={[]};
    elseif nargout==3
        varargout(3)={[]};
        varargout(2)={y_data};
        varargout(1)={x_data};
    else
        error('Incorrect number of output arguments.')
    end   
    return
end

%parse varargin input
nargs=size(varargin,2);
if nargs==1
    upsam=varargin{1};
elseif nargs==3
    upsam=varargin{1};
    if ischar(varargin{3})
        htcut=varargin{2};
        cuttype=varargin{3};
    else
        gsize=varargin{2};
        gstd=varargin{3};
    end
elseif nargs==5
    upsam=varargin{1};
    gsize=varargin{2};
    gstd=varargin{3};
    htcut=varargin{4};
    cuttype=varargin{5};
elseif nargs~=0
    error('Incorrect number of input arguments.  See "help peakfind".')
end

%check value of upsam
if nargs>0
    if upsam<1
        error('Upsampling factor must be an integer greater than or equal to 1.')
    else
        if and(nargs>2,ischar(varargin{3}))
            upsam=round(upsam);
        end
    end
end

%upsample data
if nargs>0
    xup_data=linspace(x_data(1),x_data(end),upsam*length(x_data));
    yup_data=interp1(x_data,y_data,xup_data,'cubic');
else
    xup_data=x_data;
    yup_data=y_data;
end

if nargs>2
    %check string entries
    if or(and(nargs==3,ischar(varargin{3})),nargs==5)
        if ~or(strcmp('rel',cuttype),strcmp('abs',cuttype))
            error('Invalid string entry in function call.')
        end
    end
    
    %check value of gsize and gstd
    if or(nargs==5,and(nargs==3,~ischar(varargin{3})))
        if gsize<3
            error('Gaussian window must have at least 3 points.')
        end
        if gstd<=0
            error('Gaussian STD must be greater than zero.')
        end
    end
    
    if ischar(varargin{3})
        ycv_data=yup_data;
        xcv_data=xup_data;
        stpt=1;
    else
        %define convolution filter
        filt1 = exp(-linspace(-gsize/2,gsize/2,gsize).^2/(2*gstd^2));
        conv1 = filt1/sum(filt1);
        %convolve peak data with small STD gaussian
        ycv_data_unnorm=conv(yup_data,conv1,'same');
        ycv_data_ones=conv(ones(size(yup_data)),conv1,'same');
        ycv_data=ycv_data_unnorm./ycv_data_ones;
        xcv_data=xup_data;
        %xcv_data=xup_data+1/2*(xup_data(2)-xup_data(1));
        stpt=ceil(gsize/2);
        
        ycv_data(1:stpt)=yup_data(1:stpt);
        ycv_data(end-stpt+1:end)=yup_data(end-stpt+1:end);
    end
else
    ycv_data=yup_data;
    xcv_data=xup_data;
    stpt=1;
end

%number of upsampled data points
tipn=length(ycv_data);

%construct peak finding vectors
lpts=ycv_data(stpt:end-stpt-2);
mpts=ycv_data(stpt+1:end-stpt-1);
rpts=ycv_data(stpt+2:end-stpt);
%lpts=ycv_data(stpt:end-2);
%mpts=ycv_data(stpt+1:end-1);
%rpts=ycv_data(stpt+2:end);

%find peaks
Rpeaks=zeros(1,tipn);
Lpeaks=zeros(1,tipn);
minpeaks=zeros(1,tipn);
fpeaks=zeros(1,tipn);
Rpeaks(stpt+1:end-stpt-1)=((mpts-lpts)>=0).*((rpts-mpts)<0); %triggers just right of peak
Lpeaks(stpt+1:end-stpt-1)=((mpts-lpts)>0).*((rpts-mpts)<=0); %triggers just left of peak
minpeaks(stpt+1:end-stpt-1)=((mpts-lpts)<=0).*((rpts-mpts)>0);
%Rpeaks(stpt+1:end-stpt)=((mpts-lpts)>=0).*((rpts-mpts)<0); %triggers just right of peak
%Lpeaks(stpt+1:end-stpt)=((mpts-lpts)>0).*((rpts-mpts)<=0); %triggers just left of peak
%minpeaks(stpt+1:end-stpt)=((mpts-lpts)<=0).*((rpts-mpts)>0);

%use minimums and multiple maximums to find best guess of peak
minlist=find(minpeaks==1);

%triggers when there is a flat landscape with a single trough
if and(sum(Rpeaks)==0,sum(Lpeaks)==0)
    disp('Warning: found only a trough.')
    npeaks=0;
    peakspos=zeros(1,0);
    xcv_data=x_data;
    ycv_data=y_data;
    return
end

%find and average doubly-occupied peaks
for j=1:length(minlist)+1
    %find bounds between which to look for two pseudo-peaks
    if j==1
        lbound=1;
        if isempty(minlist)
            rbound=tipn;
        else
            rbound=minlist(j);
        end
    elseif j==length(minlist)+1
        lbound=minlist(j-1);
        rbound=length(ycv_data);
    else
        lbound=minlist(j-1);
        rbound=minlist(j);
    end
    
    %look for peaks
    peakmask=zeros(1,tipn);
    peakmask(lbound:rbound)=1;
    
    Lpeak=find(Lpeaks.*peakmask==1,1);
    Rpeak=find(Rpeaks.*peakmask==1,1);
    fpeaks(round((Lpeak+Rpeak)/2))=1;
end

%do a final check against the peaks in the original up-sampled data
peakspos=find(fpeaks==1);
npeaks=sum(fpeaks);

%remove peaks if they are not high enough
if and(npeaks>1,exist('cuttype','var'))
    ymin=min(ycv_data);
    ymax=max(ycv_data);
    
    if strcmp('rel',cuttype)
        cutpeaks=((ycv_data(peakspos)-ymin)/(ymax-ymin))<htcut;
    else
        cutpeaks=ycv_data(peakspos)<(htcut+ymin);
    end
    npeaks=npeaks-sum(cutpeaks);
    peakspos=peakspos(cutpeaks==0);
end

if nargout==1
    varargout(1)={xcv_data(peakspos)};
elseif nargout==3
    varargout(3)={peakspos};
    varargout(2)={ycv_data};
    varargout(1)={xcv_data};
else
    error('Incorrect number of output arguments.')
end

return


Contact us