bleeck@4: % function out = PeakPicker(sig_in, params)
bleeck@4: %
bleeck@4: %   Find the peaks of a signal and their neighbours
bleeck@4: %
bleeck@4: %   INPUT VALUES:
bleeck@4: %       sig_in                  Input signal
bleeck@4: %         threshold     dynamic threshold. Off if not used
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) 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,threshold)
bleeck@4: % threshold is the percentage of the maximum value of the signal,
bleeck@4: % under which a peak is not counted
bleeck@4: 
bleeck@4: if nargin <2
bleeck@4: 	threshold=0;
bleeck@4: end
bleeck@4: 
bleeck@4: plot_switch = 0;
bleeck@4: 
bleeck@4: % the original values befor filtering
bleeck@4: orig_values = getdata(sig_in)';
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: 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=[];
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: % maxima = cell(1, length(peakpos_x));
bleeck@4: maxima = [];
bleeck@4: 
bleeck@4: if plot_switch
bleeck@4: 	figure(123);
bleeck@4: 	clf
bleeck@4: 	plot(sig_in,'k');
bleeck@4: 	hold on;
bleeck@4: end
bleeck@4: 
bleeck@4: threshold_val=threshold*max(sig_in);
bleeck@4: counter=1;
bleeck@4: 
bleeck@4: % find the left end right minima that belong to a maximum
bleeck@4: for i=1:length(peakpos_x)
bleeck@4: 	y_val= orig_values(peakpos_x(i));
bleeck@4: 	
bleeck@4: 	if y_val< threshold_val
bleeck@4: 		continue
bleeck@4: 	end
bleeck@4: 
bleeck@4: 	maxima{counter}.y =y_val;
bleeck@4: 	maxima{counter}.x = peakpos_x(i);
bleeck@4: 	maxima{counter}.t = bin2time(sig_in, maxima{i}.x);
bleeck@4: 	maxima{counter}.fre = 1/maxima{counter}.t;
bleeck@4: 	
bleeck@4: 	if plot_switch
bleeck@4: 		plot(maxima{counter}.x,maxima{counter}.y,'go');
bleeck@4: 	end    
bleeck@4: 	% find left and right minimum for this maximum
bleeck@4: 	maxima{counter}.left.x = max(min_x([min_x < maxima{counter}.x]));
bleeck@4: 	if isempty(maxima{counter}.left.x)
bleeck@4: 		maxima{counter}.left.x = 1;
bleeck@4: 		maxima{counter}.left.t = 0;
bleeck@4: 		maxima{counter}.left.y = orig_values(maxima{counter}.left.x);
bleeck@4: 	else
bleeck@4: 		maxima{counter}.left.y = orig_values(maxima{counter}.left.x);
bleeck@4: 		maxima{counter}.left.t = bin2time(sig_in, maxima{counter}.left.x);
bleeck@4: 	end
bleeck@4: 	maxima{counter}.right.x = min(min_x([min_x > maxima{counter}.x]));
bleeck@4: 	if isempty(maxima{counter}.right.x)
bleeck@4: 		maxima{counter}.right.x = length(orig_values);
bleeck@4: 		maxima{counter}.right.t = 0;
bleeck@4: 		maxima{counter}.right.y = orig_values(maxima{counter}.right.x);
bleeck@4: 	else
bleeck@4: 		maxima{counter}.right.y = orig_values(maxima{counter}.right.x);
bleeck@4: 		maxima{counter}.right.t = bin2time(sig_in, maxima{counter}.right.x);
bleeck@4: 	end
bleeck@4: 	
bleeck@4: 	if plot_switch
bleeck@4: 		plot(maxima{counter}.right.x,maxima{counter}.right.y,'ro');
bleeck@4: 		plot(maxima{counter}.left.x,maxima{counter}.left.y,'ro');
bleeck@4: 	end
bleeck@4: 	counter=counter+1;
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: % umrechnung in die Darstellung, die wir brauchen:
bleeck@4: for i=1:counter-1
bleeck@4: 	logtime=maxima{i}.t;
bleeck@4: 	time=logtime2time(logtime);
bleeck@4: 	maxima{i}.t=time;
bleeck@4: 	maxima{i}.fre=1/time;
bleeck@4: 	maxima{i}.y=gettimevalue(sig_in,logtime);
bleeck@4: 
bleeck@4: 	left=maxima{i}.left;
bleeck@4: 	lefttime=logtime2time(left.t);
bleeck@4: 	maxima{i}.left.t=lefttime;
bleeck@4: 	maxima{i}.left.fre=1/lefttime;
bleeck@4: 	maxima{i}.left.y=gettimevalue(sig_in,left.t);
bleeck@4: 
bleeck@4: 	right=maxima{i}.right;
bleeck@4: 	righttime=logtime2time(right.t);
bleeck@4: 	maxima{i}.right.t=righttime;
bleeck@4: 	maxima{i}.right.fre=1/righttime;
bleeck@4: 	maxima{i}.right.y=gettimevalue(sig_in,right.t);
bleeck@4: end
bleeck@4: 
bleeck@4: 
bleeck@4: out = maxima;
bleeck@4: 
bleeck@4: 
bleeck@4: 
bleeck@4: function time=logtime2time(logtime)
bleeck@4: 	time=f2f(logtime,0,0.035,0.001,0.035,'linlog');
bleeck@4: