bleeck@4: %  
bleeck@4: % function out = PeakPicker(sig_in, params)
bleeck@4: %
bleeck@4: %   Find the peaks of a signal
bleeck@4: %
bleeck@4: %   INPUT VALUES:
bleeck@4: %       sig_in                  Input signal
bleeck@4: %       params.dyn_thresh       dynamic threshold. Off if not used
bleeck@4: %       params.smooth_sigma     sigma for smoothing
bleeck@4: %                               
bleeck@4: %       
bleeck@4: % 		
bleeck@4: %
bleeck@4: %   RETURN VALUE:
bleeck@4: %       out is an array of a struct
bleeck@4: %       out.x                    x position of the Peak
bleeck@4: %       out.t                    according time value
bleeck@4: %		out.y                    y value of the peak
bleeck@4: %       out.left.[x,t,y]         left Minumum
bleeck@4: %       out.right.[x,t,y]        right Minumum
bleeck@4: % 
bleeck@4: % (c) 2003, University of Cambridge, Medical Research Council 
bleeck@4: % Christoph Lindner 
bleeck@4: % (c) 2011, University of Southampton
bleeck@4: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@4: % download of current version is on the soundsoftware site: 
bleeck@4: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@4: % documentation and everything is on http://www.acousticscale.org
bleeck@4: 
bleeck@4: function out = PeakPicker(sig_in, params)
bleeck@4: 
bleeck@4: if nargin<2
bleeck@4:     params=[];
bleeck@4: end
bleeck@4: 
bleeck@4: % % % % ----- Other Parameters -----
bleeck@4: % % % % Lowpass param for the higpassfilter
bleeck@4: % % % LP_sigma_for_HP_filter = getnrpoints(sig_in)/7;
bleeck@4: % % % LP_sigma_for_smooth = 3;
bleeck@4: % % % % min width of a peak: upper_threshold+lower_thresh
bleeck@4: % % % upper_thresh =  0.02*getnrpoints(sig_in);
bleeck@4: % % % lower_thresh =  0.03*getnrpoints(sig_in);
bleeck@4: 
bleeck@4: % x is index or position in vector
bleeck@4: % y is value of the
bleeck@4: % t is the time domain wich is assigned to the x dimension
bleeck@4: 
bleeck@4: % the original values befor filtering
bleeck@4: orig_values = getdata(sig_in)';
bleeck@4: 
bleeck@4: if isfield(params,'smooth_sigma') 
bleeck@4:     if (params.smooth_sigma~=0)
bleeck@4:         % smooth the curve to kill small side peaks 
bleeck@4:         sig_in = smooth(sig_in, params.smooth_sigma);
bleeck@4:     end
bleeck@4: end
bleeck@4: 
bleeck@4: values=getdata(sig_in)';
bleeck@4: 
bleeck@4: % ------------------- Find the local maxima ------------------------------
bleeck@4: % find x positions of ALL local maxima, incl. zero!!
bleeck@4: max_x = find((values >= [0 values(1:end-1)]) & (values > [values(2:end) 0]));
bleeck@4: if isfield(params,'smooth_sigma') 
bleeck@4:     if (params.smooth_sigma~=0)
bleeck@4:         % smoothing might have shifted the positions of maxima.
bleeck@4:         % Therefore the maximum is the highest of the neighbours in
bleeck@4:         % distance +- smooth_sigma
bleeck@4:         for i=1:length(max_x)
bleeck@4:             start = max_x(i)-params.smooth_sigma;
bleeck@4:             if start<1
bleeck@4:                 start=1;
bleeck@4:             end
bleeck@4:             stop = max_x(i)+params.smooth_sigma;
bleeck@4:             if stop>length(orig_values)
bleeck@4:                 stop=length(orig_values);
bleeck@4:             end
bleeck@4:             m = find(orig_values(start:stop) == max(orig_values(start:stop)));
bleeck@4:             m=m(1);
bleeck@4:             max_x(i)=start-1+m;
bleeck@4:         end
bleeck@4:     end
bleeck@4: end
bleeck@4: max_y = orig_values(max_x);
bleeck@4: orig_max_y = orig_values(max_x);
bleeck@4: 
bleeck@4: % ------------------- Find the local minima -----------------------------
bleeck@4: min_x = find((values < [inf values(1:end-1)]) & (values <= [values(2:end) inf]));
bleeck@4: min_y = values(min_x);
bleeck@4: 
bleeck@4: 
bleeck@4: % peakpos_x=zeros(1,length(max_x));
bleeck@4: peakpos_x=[];
bleeck@4: for i=1:length(max_x),
bleeck@4:     % only take the highest peak 
bleeck@4:     my = [max_y==max(max_y)];  % find the highest peak
bleeck@4: %     peakpos_x(i) = max_x(my);   % x pos of highest peak
bleeck@4:     peakpos_x = [peakpos_x max_x(my)]; 
bleeck@4:     max_y = max_y([max_y<max(max_y)]);   % del max value in y domain
bleeck@4:     max_x = max_x([max_x ~= peakpos_x(end)]); % and in x domain
bleeck@4: end
bleeck@4: peakpos_y = orig_values(peakpos_x);  % extract the y vector 
bleeck@4: 
bleeck@4: % --------------- Dynamic Threshold ---------------------
bleeck@4: % works relativ to the mean
bleeck@4: if isfield(params,'dyn_thresh') 
bleeck@4:     if (params.dyn_thresh~=0)
bleeck@4:         % dynamic thresholding 
bleeck@4:         m = mean(orig_values);  
bleeck@4:         dthr = params.dyn_thresh.*m;
bleeck@4:         peakpos_x = peakpos_x([peakpos_y>=dthr]);
bleeck@4:         peakpos_y = peakpos_y([peakpos_y>=dthr]);
bleeck@4:     end
bleeck@4: end
bleeck@4: 
bleeck@4: maxima = cell(1, length(peakpos_x));
bleeck@4: % find the left end right minima that belong to a maximum
bleeck@4: for i=1:length(peakpos_x)
bleeck@4:     maxima{i}.x = peakpos_x(i);
bleeck@4:     maxima{i}.t = bin2time(sig_in, maxima{i}.x);
bleeck@4:     maxima{i}.y = orig_values(peakpos_x(i));
bleeck@4:     
bleeck@4:     % find left and right minimum for this maximum
bleeck@4:     maxima{i}.left.x = max(min_x([min_x < maxima{i}.x]));
bleeck@4:     if isempty(maxima{i}.left.x)
bleeck@4:         maxima{i}.left.x = 1;
bleeck@4:         maxima{i}.left.t = 0;
bleeck@4:         maxima{i}.left.y = orig_values(maxima{i}.left.x);
bleeck@4:     else
bleeck@4:         maxima{i}.left.y = orig_values(maxima{i}.left.x);
bleeck@4:         maxima{i}.left.t = bin2time(sig_in, maxima{i}.left.x);
bleeck@4:     end
bleeck@4:     maxima{i}.right.x = min(min_x([min_x > maxima{i}.x]));
bleeck@4:     if isempty(maxima{i}.right.x)
bleeck@4:         maxima{i}.right.x = length(orig_values);
bleeck@4:         maxima{i}.right.t = 0;
bleeck@4:         maxima{i}.right.y = orig_values(maxima{i}.right.x);
bleeck@4:     else
bleeck@4:         maxima{i}.right.y = orig_values(maxima{i}.right.x);
bleeck@4:         maxima{i}.right.t = bin2time(sig_in, maxima{i}.right.x);
bleeck@4:     end
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: out = maxima;
bleeck@4: 
bleeck@4: