tomwalters@0: % tool
tomwalters@0: % 
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %  
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %
tomwalters@0: % 
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
bleeck@3: 
tomwalters@0: 
tomwalters@0: function [newsig,strobe_points,threshold]=adaptivthreshold(sig,options)
tomwalters@0: 
tomwalters@0: if nargin < 2
tomwalters@0: 	options=[];
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % the way, in which the strobes are calculated
tomwalters@0: if isfield(options,'strobe_decay_time')
tomwalters@0: 	strobe_decay_time=options.strobe_decay_time;
tomwalters@0: else
tomwalters@0: 	strobe_decay_time=0.04;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % the time constant of the strobe decay
tomwalters@0: if isfield(options,'strobe_decay_time')
tomwalters@0: 	strobe_decay_time=options.strobe_decay_time;
tomwalters@0: else
tomwalters@0: 	strobe_decay_time=0.04;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % the timeout, after that time a strobe must occure
tomwalters@0: if isfield(options,'timeout')
tomwalters@0: 	timeout=options.timeout;
tomwalters@0: else
tomwalters@0: 	strobe_decay_time=0.01;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % the strobe lag, ie the time that is waited for the next strobe to occure
tomwalters@0: if isfield(options,'strobe_lag')
tomwalters@0: 	strobe_lag=options.strobe_lag;
tomwalters@0: else
tomwalters@0: 	strobe_lag=0.005;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % the threshold of a strobe
tomwalters@0: if isfield(options,'threshold')
tomwalters@0: 	threshold=options.threshold;
tomwalters@0: else
tomwalters@0: 	threshold=0.04;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % if we want to start with a different threshold
tomwalters@0: if isfield(options,'startingthreshold')
tomwalters@0: 	current_threshold=options.startingthreshold;
tomwalters@0: else
tomwalters@0: 	current_threshold=0;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % if some graphic should be displayed during run
tomwalters@0: if isfield(options,'grafix')
tomwalters@0: 	grafix=options.grafix;
tomwalters@0: else
tomwalters@0: 	grafix=0; % 
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % % if the result shell be filtered by this lowpass filter:
tomwalters@0: % if isfield(options,'lowpass') 
tomwalters@0: %     lowpassvalue=options.lowpass;
tomwalters@0: % else
tomwalters@0: %     lowpassvalue=-1; % 
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: sr=getsr(sig);
tomwalters@0: newsig=sig;
tomwalters@0: newsig=setname(newsig,sprintf('adaptive threshold of %s',getname(sig)));
tomwalters@0: % newvals=getvalues(sig);
tomwalters@0: threshold=sig;
tomwalters@0: tresval=getvalues(sig);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: sigvals=getvalues(sig);
tomwalters@0: options.current_threshold=current_threshold;
tomwalters@0: options.sr=sr;
tomwalters@0: 
tomwalters@0: switch options.strobe_criterion
tomwalters@0: 	case 'parabola'
tomwalters@0: 		[strobe_points,tresval,newvals]=doadaptivethresholdparabola(sigvals,options);
tomwalters@0: 	case 'threshold'
tomwalters@0: 	case 'peak'
tomwalters@0: 		[strobe_points,tresval,newvals]=doadaptivethresholdpeak(sigvals,options);
tomwalters@0: 	case 'peak_shadow+'
tomwalters@0: 		[strobe_points,tresval,newvals]=doadaptivethresholdpeakshadowplus(sigvals,options);
tomwalters@0: 	case 'peak_shadow-'
tomwalters@0: 	case 'delta_gamma'
tomwalters@0: 	case 'adaptive'
tomwalters@0: 		[strobe_points,tresval,newvals]=doadaptivethresholdadaptive(sigvals,options);
tomwalters@0: 	case 'bunt'
tomwalters@0: 		[strobe_points,tresval,newvals]=doadaptivethresholdbunt(sigvals,options);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: strobe_points=bin2time(sig,strobe_points);
tomwalters@0: 
tomwalters@0: newsig=setvalues(newsig,newvals);
tomwalters@0: threshold=setvalues(threshold,tresval);
tomwalters@0: 
tomwalters@0: % % do a low pass filtering of the result to smooth it
tomwalters@0: % if lowpassvalue > -1
tomwalters@0: %     if lowpassvalue > 2000
tomwalters@0: %         lowpassvalue=2000;
tomwalters@0: %     end  
tomwalters@0: %     newsig=lowpass(newsig,lowpassvalue);
tomwalters@0: %     newsig=halfwayrectify(newsig);
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0: 	clf
tomwalters@0: 	plot(sig); hold on
tomwalters@0: 	ax=axis;
tomwalters@0: 	plot(newsig,'r');
tomwalters@0: 	plot(threshold,'g');
tomwalters@0: 	axis(ax);
tomwalters@0: end
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: function [strobe_points,tresval,newvals]=doadaptivethresholdparabola(sigvals,options)
tomwalters@0: % the threshold is calculated relativ to the hight of the last strobe
tomwalters@0: % the sample rate is needed for the calculation of the next thresholds
tomwalters@0: % for time_constant_alpha this is a linear decreasing function that goes
tomwalters@0: % from the maximum value to 0 in the time_constant
tomwalters@0: 
tomwalters@0: current_threshold=options.current_threshold;
tomwalters@0: sr=options.sr;
tomwalters@0: last_strobe=-inf;
tomwalters@0: 
tomwalters@0: tresval=zeros(size(sigvals));
tomwalters@0: newvals=zeros(size(sigvals));
tomwalters@0: nr=length(sigvals);
tomwalters@0: strobe_points=zeros(100,1);
tomwalters@0: 
tomwalters@0: %when the last strobe occured
tomwalters@0: % last_strobe=-inf;
tomwalters@0: last_threshold_value=0;
tomwalters@0: 
tomwalters@0: % copy options to save time
tomwalters@0: h=options.parabel_heigth;
tomwalters@0: w=options.parabel_width;
tomwalters@0: %     w=options.parabel_width_in_cycles*1/options.current_cf;
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: 
tomwalters@0: counter=1;
tomwalters@0: for ii=1:nr
tomwalters@0: 	current_val=sigvals(ii);
tomwalters@0: 	current_time=ii/sr;
tomwalters@0: 	
tomwalters@0: 	if current_val>current_threshold
tomwalters@0: 		new_val=current_val-current_threshold;
tomwalters@0: 		current_threshold=current_val;
tomwalters@0: 		strobe_points(counter)=ii;
tomwalters@0: 		counter=counter+1;
tomwalters@0: 		last_strobe=ii;
tomwalters@0: 		last_threshold_value=current_threshold;
tomwalters@0: 	else
tomwalters@0: 		new_val=0;
tomwalters@0: 	end
tomwalters@0: 	tresval(ii)=current_threshold;
tomwalters@0: 	
tomwalters@0: 	time_since_last_strobe=(ii-last_strobe)/sr;
tomwalters@0: 	if time_since_last_strobe > w % linear falling threshold
tomwalters@0: 		const_decay=last_threshold_value/sr/strobe_decay_time;
tomwalters@0: 		current_threshold=current_threshold-const_decay;
tomwalters@0: 		current_threshold=max(0,current_threshold);
tomwalters@0: 	else    % parabel for the first time y=a(x+b)^2+c
tomwalters@0: 		a=4*(1-h)/(w*w);
tomwalters@0: 		b=-w/2;
tomwalters@0: 		c=h;
tomwalters@0: 		factor=a*(time_since_last_strobe + b) ^2+c;
tomwalters@0: 		current_threshold=last_threshold_value*factor;
tomwalters@0: 	end    
tomwalters@0: 	
tomwalters@0: 	current_threshold=max(0,current_threshold);     % cant be smaller then 0
tomwalters@0: 	
tomwalters@0: 	newvals(ii)=new_val;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function [strobe_points,tresval,newvals]=doadaptivethresholdpeak(sigvals,options)
tomwalters@0: % finds every single local maximum
tomwalters@0: sr=options.sr;
tomwalters@0: tresval=zeros(size(sigvals));
tomwalters@0: newvals=zeros(size(sigvals));
tomwalters@0: strobe_points=zeros(100,1);
tomwalters@0: sig=signal(sigvals);
tomwalters@0: sig=setsr(sig,sr);
tomwalters@0: [t,h]=getlocalmaxima(sig);
tomwalters@0: strobe_points=t*sr;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: function [strobe_points,tresval,newvals]=doadaptivethresholdpeakshadowplus(sigvals,options)
tomwalters@0: % finds every single peak and starts from that an falling threshold
tomwalters@0: 
tomwalters@0: current_threshold=options.current_threshold;
tomwalters@0: sr=options.sr;
tomwalters@0: tresval=zeros(size(sigvals));
tomwalters@0: newvals=zeros(size(sigvals));
tomwalters@0: nr=length(sigvals);
tomwalters@0: strobe_points=zeros(100,1);
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: last_strobe=-inf;
tomwalters@0: counter=1;
tomwalters@0: last_threshold_value=0;
tomwalters@0: 
tomwalters@0: for ii=2:nr-1
tomwalters@0: 	current_val=sigvals(ii);
tomwalters@0: 	current_time=ii/sr;
tomwalters@0: 	
tomwalters@0: 	if current_val>current_threshold
tomwalters@0: 		if sigvals(ii) > sigvals(ii-1) && sigvals(ii) > sigvals(ii+1)
tomwalters@0: 			new_val=current_val-current_threshold;
tomwalters@0: 			current_threshold=current_val;
tomwalters@0: 			strobe_points(counter)=ii;
tomwalters@0: 			counter=counter+1;
tomwalters@0: 			last_strobe=ii;
tomwalters@0: 			last_threshold_value=current_threshold;
tomwalters@0: 		end
tomwalters@0: 	end
tomwalters@0: 	const_decay=last_threshold_value/sr/strobe_decay_time;
tomwalters@0: 	current_threshold=current_threshold-const_decay;
tomwalters@0: 	current_threshold=max(0,current_threshold);
tomwalters@0: 	tresval(ii)=current_threshold;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: function [strobe_points,tresval,newvals]=doadaptivethresholdadaptive(sigvals,options)
tomwalters@0: % the threshold is calculated relativ to the hight of the last strobe
tomwalters@0: % the sample rate is needed for the calculation of the next thresholds
tomwalters@0: % for time_constant_alpha this is a linear decreasing function that goes
tomwalters@0: % from the maximum value to 0 in the time_constant
tomwalters@0: 
tomwalters@0: current_threshold=options.current_threshold;
tomwalters@0: sr=options.sr;
tomwalters@0: last_strobe=-inf;
tomwalters@0: 
tomwalters@0: tresval=zeros(size(sigvals));
tomwalters@0: newvals=zeros(size(sigvals));
tomwalters@0: nr=length(sigvals);
tomwalters@0: strobe_points=zeros(100,1);
tomwalters@0: 
tomwalters@0: %when the last strobe occured
tomwalters@0: % last_strobe=-inf;
tomwalters@0: last_threshold_value=0;
tomwalters@0: 
tomwalters@0: % copy options to save time
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: 
tomwalters@0: bunt=0.5;
tomwalters@0: 
tomwalters@0: % decay_time=options.strobe_decay_time;
tomwalters@0: % threshold_decay_constant=power(0.5,1./(options.strobe_decay_time*sr));
tomwalters@0: 
tomwalters@0: slope_coefficient=options.slope_coefficient;
tomwalters@0: slope_coefficient=0.0005;
tomwalters@0: minoffset=0.2;
tomwalters@0: 
tomwalters@0: threshold_decay_offset=-1/(options.strobe_decay_time*sr);
tomwalters@0: default_threshold_decay_offset=threshold_decay_offset;
tomwalters@0: 
tomwalters@0: counter=1;
tomwalters@0: for ii=1:nr
tomwalters@0: 	current_val=sigvals(ii);
tomwalters@0: 	current_time=ii/sr;
tomwalters@0: 	
tomwalters@0: 	if current_val>current_threshold
tomwalters@0: 		new_val=current_val-current_threshold;
tomwalters@0: 		current_threshold=current_val;
tomwalters@0: 		strobe_points(counter)=ii;
tomwalters@0: 		counter=counter+1;
tomwalters@0: 		time_offset=ii-last_strobe;	% time since last strobe
tomwalters@0: 		last_strobe=ii;
tomwalters@0: 		
tomwalters@0: 		amplitude_offset=current_threshold-last_threshold_value;
tomwalters@0: 		
tomwalters@0: 		last_threshold_value=current_threshold;
tomwalters@0: 		
tomwalters@0: 		current_bunt=0;
tomwalters@0: 		% 		if amplitude_offset>0
tomwalters@0: 		% 			current_bunt=amplitude_offset/1;
tomwalters@0: 		% 		else
tomwalters@0: 		% 			current_bunt=0;
tomwalters@0: 		% 		end			
tomwalters@0: 		current_threshold=current_threshold+current_bunt+minoffset;
tomwalters@0: 		
tomwalters@0: 		offset=amplitude_offset/time_offset*slope_coefficient;
tomwalters@0: 		
tomwalters@0: 		threshold_decay_offset=threshold_decay_offset+offset;
tomwalters@0: 		% 		threshold_decay_constant=power(0.5,1./(decay_time*sr));
tomwalters@0: 	else
tomwalters@0: 		new_val=0;
tomwalters@0: 	end
tomwalters@0: 	tresval(ii)=current_threshold;
tomwalters@0: 	time_since_last_strobe=(ii-last_strobe)/sr;
tomwalters@0: 	
tomwalters@0: 	
tomwalters@0: 	% 	current_threshold=current_threshold*threshold_decay_constant;
tomwalters@0: 	current_threshold=current_threshold+threshold_decay_offset;
tomwalters@0: 	current_threshold=max(current_threshold,0);
tomwalters@0: 	
tomwalters@0: 	if time_since_last_strobe>0.035
tomwalters@0: 		current_threshold=0;
tomwalters@0: 		threshold_decay_offset=default_threshold_decay_offset;
tomwalters@0: 	end
tomwalters@0: 	
tomwalters@0: 	newvals(ii)=new_val;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%% BUNT
tomwalters@0: function [strobe_points,tresval,newvals]=doadaptivethresholdbunt(sigvals,options)
tomwalters@0: % the bunt is relative to the last peak hight
tomwalters@0: 
tomwalters@0: current_threshold=options.current_threshold;
tomwalters@0: sr=options.sr;
tomwalters@0: last_strobe=-inf;
tomwalters@0: 
tomwalters@0: tresval=zeros(size(sigvals));
tomwalters@0: newvals=zeros(size(sigvals));
tomwalters@0: nr=length(sigvals);
tomwalters@0: strobe_points=zeros(100,1);
tomwalters@0: 
tomwalters@0: %when the last strobe occured
tomwalters@0: last_strobe_time=-inf;
tomwalters@0: last_threshold_value=0;
tomwalters@0: 
tomwalters@0: % copy options to save time
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: 
tomwalters@0: % wait that many cycles to confirm a strobe
tomwalters@0: wait_time=options.wait_cycles/options.cf;
tomwalters@0: 
tomwalters@0: % if waited for too long, then strobe anyhow after some passed candidates:
tomwalters@0: wait_candidate_max=options.nr_strobe_candidates;
tomwalters@0: current_jumped_candidates=1;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: bunt=options.bunt;
tomwalters@0: 
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: last_val=sigvals(1);
tomwalters@0: 
tomwalters@0: counter=1;
tomwalters@0: for ii=2:nr-1
tomwalters@0: 	current_val=sigvals(ii);
tomwalters@0: 	next_val=sigvals(ii+1);
tomwalters@0: 	current_time=ii/sr;
tomwalters@0: 	if current_val>=current_threshold  % above threshold -> criterium for strobe
tomwalters@0: 		current_threshold=current_val;
tomwalters@0: 		if current_val > last_val && current_val > next_val  % only strobe on local maxima
tomwalters@0: 
tomwalters@0: 			% so far its a candidate, but is it a real strobe?
tomwalters@0: 			% look if there was a strobe in the past, that is deleted
tomwalters@0: 			if (counter>1 && current_time-strobe_time(counter-1)<wait_time) %&& current_threshold >last_threshold_value
tomwalters@0: 				% if its too long in the past, fire anyway
tomwalters@0: 				if current_jumped_candidates > wait_candidate_max	
tomwalters@0: 					current_jumped_candidates=1;	% reset counter
tomwalters@0: 				else
tomwalters@0: 					current_jumped_candidates=current_jumped_candidates+1;
tomwalters@0: 					counter=counter-1;	% delete the last one
tomwalters@0: 				end
tomwalters@0: 			else
tomwalters@0: 				current_jumped_candidates=current_jumped_candidates+1;
tomwalters@0: 			end				
tomwalters@0: 			
tomwalters@0: 				
tomwalters@0: 			% take this one as a new one
tomwalters@0: 			strobe_points(counter)=ii;
tomwalters@0: 			strobe_time(counter)=ii/sr;
tomwalters@0: 			counter=counter+1;	% strobecounter
tomwalters@0: 			current_threshold=current_threshold*options.bunt;	%increase threshold
tomwalters@0: 			
tomwalters@0: 			last_strobe_time=ii/sr; % anyhow, its a candidate
tomwalters@0: 			last_threshold_value=current_threshold;
tomwalters@0: 			
tomwalters@0: 		end
tomwalters@0: 	end
tomwalters@0: 	tresval(ii)=current_threshold;
tomwalters@0: 	
tomwalters@0: 	const_decay=last_threshold_value/sr/strobe_decay_time;
tomwalters@0: 	current_threshold=current_threshold-const_decay;
tomwalters@0: 	
tomwalters@0: 	current_threshold=max(current_threshold,0);
tomwalters@0: 	last_val=current_val;
tomwalters@0: end
tomwalters@0: