tomwalters@0: % generating function for 'aim-mat'
tomwalters@0: % 
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %  
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %
tomwalters@0: % 
tomwalters@0: % (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: 
tomwalters@0: function [allstrobeprocesses,allthresholds]=gen_sf2003(nap,strobeoptions)
tomwalters@0: 
tomwalters@0: % find in each channel each strobe and gives it back as an list of
tomwalters@0: % structures of strobes
tomwalters@0: 
tomwalters@0: % possible values of strobe_criterion:
tomwalters@0: % switch strobeoptions.strobe_criterion
tomwalters@0: % 	case 'threshold'
tomwalters@0: % 	case 'peak'
tomwalters@0: % 	case 'temporal_shadow'
tomwalters@0: % 	case 'local_maximum'
tomwalters@0: % 	case 'delta_gamma'
tomwalters@0: % 	case 'parabola'
tomwalters@0: % 	case 'bunt'
tomwalters@0: % 	case 'adaptive'
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: if strcmp(strobeoptions.strobe_criterion,'interparabola')
tomwalters@0: 	[allstrobeprocesses,allthresholds]=do_interparabola(nap,strobeoptions);
tomwalters@0: 	return
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % non interchannel methods:
tomwalters@0: % buffer for all thrsholds afterwards. Exact copy of the signal:
tomwalters@0: allthresholds=nap;
tomwalters@0: cfs=getcf(nap);
tomwalters@0: 
tomwalters@0: waithand=waitbar(0,'generating strobes');
tomwalters@0: nr_channels=getnrchannels(nap);
tomwalters@0: for ii=1:nr_channels
tomwalters@0: 	waitbar(ii/nr_channels);
tomwalters@0: 	single_channel=getsinglechannel(nap,ii);   
tomwalters@0: 	current_cf=cfs(ii);
tomwalters@0: 	strobeoptions.current_cf=current_cf;
tomwalters@0: 	
tomwalters@0: 	% here they are calculated:
tomwalters@0: 	[strobe_points,threshold]=findstrobes(single_channel,strobeoptions);
tomwalters@0: 	
tomwalters@0: 	% return values: strobes in sec
tomwalters@0: 	% threshold= signal indicating the adaptive threshold
tomwalters@0: 	
tomwalters@0: 	strobe_points=strobe_points(find(strobe_points>0));
tomwalters@0: 	strobe_vals=gettimevalue(single_channel,strobe_points);
tomwalters@0: 	
tomwalters@0: 	thresvals=getvalues(threshold);
tomwalters@0: 	
tomwalters@0: 	% make shure, they are in one column
tomwalters@0: 	if size(strobe_points,1) > size(strobe_points,2)
tomwalters@0: 		strobe_points=strobe_points';
tomwalters@0: 		strobe_vals=strobe_vals';
tomwalters@0: 	end
tomwalters@0: 	allstrobeprocesses{ii}.strobes=strobe_points;
tomwalters@0: 	allstrobeprocesses{ii}.strobe_vals=strobe_vals;
tomwalters@0: 	allthresholds=setsinglechannel(allthresholds,ii,thresvals);
tomwalters@0: end
tomwalters@0: close(waithand);
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % interchannel methods:
tomwalters@0: function [all_processes,thresholds]=do_interparabola(nap,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: % with interchannel interaction: A strobe in one channel reduces the
tomwalters@0: % threshold in the neighbouring channels
tomwalters@0: 
tomwalters@0: nr_channels=getnrchannels(nap);
tomwalters@0: 
tomwalters@0: current_threshold=zeros(nr_channels,1);
tomwalters@0: sr=getsr(nap);
tomwalters@0: last_strobe=ones(nr_channels,1)* -inf;
tomwalters@0: 
tomwalters@0: napvals=getvalues(nap);
tomwalters@0: tresval=zeros(size(napvals));
tomwalters@0: nr_dots=length(napvals);
tomwalters@0: 
tomwalters@0: strobe_points=zeros(nr_channels,1000); % makes things faster
tomwalters@0: strobe_times=zeros(size(strobe_points)); % makes things faster
tomwalters@0: 
tomwalters@0: %when the last strobe occured
tomwalters@0: last_strobe_time=ones(nr_channels,1)* -inf;
tomwalters@0: last_threshold_value=zeros(nr_channels,1);
tomwalters@0: last_val=napvals(:,1);
tomwalters@0: 
tomwalters@0: cfs=getcf(nap);
tomwalters@0: % copy options to save time
tomwalters@0: h=options.parabel_heigth;
tomwalters@0: wnull=options.parabel_width_in_cycles*1./cfs;
tomwalters@0: w_variabel=wnull;
tomwalters@0: 
tomwalters@0: strobe_decay_time=options.strobe_decay_time;
tomwalters@0: times_per_ms=round(sr*0.005); % how often the bar should be updated
tomwalters@0: counter=ones(nr_channels,1);
tomwalters@0: 
tomwalters@0: waithand=waitbar(0,'generating Interchannel Strobes');
tomwalters@0: 
tomwalters@0: for ii=2:nr_dots-1
tomwalters@0: 	current_time=ii/sr;
tomwalters@0: 	if mod(ii,times_per_ms)==0
tomwalters@0: 		waitbar(ii/nr_dots);
tomwalters@0: 	end
tomwalters@0: 
tomwalters@0: 	for jj=1:nr_channels
tomwalters@0: 		current_val=napvals(jj,ii);
tomwalters@0: 		next_val=napvals(jj,ii+1);
tomwalters@0: 		
tomwalters@0: 		if current_val>=current_threshold(jj)  % above threshold -> criterium for strobe
tomwalters@0: 			current_threshold(jj)=current_val;
tomwalters@0: 			if current_val > last_val(jj) && current_val > next_val  % only strobe on local maxima
tomwalters@0: 				% take this one as a new one
tomwalters@0: 				strobe_points(jj,counter(jj))=ii;
tomwalters@0: 				strobe_time(jj,counter(jj))=ii/sr;
tomwalters@0: 				counter(jj)=counter(jj)+1;	% strobecounter
tomwalters@0: 				
tomwalters@0: 				last_strobe_time(jj)=ii/sr; % anyhow, its a candidate
tomwalters@0: 				last_threshold_value(jj)=current_threshold(jj);
tomwalters@0: 				a=4*(1-h)/(wnull(jj)*wnull(jj));
tomwalters@0: 				b=-wnull(jj)/2;
tomwalters@0: 				w_variabel(jj)=wnull(jj)-(current_threshold(jj)-2*a*b)/(2*a);
tomwalters@0: 				
tomwalters@0: 			% the interchannel action: reduce the threshold in adjacent channel
tomwalters@0: 				if jj>1
tomwalters@0: 					current_threshold(jj-1)=current_threshold(jj-1)/1.2;
tomwalters@0: 					last_threshold_value(jj-1)=last_threshold_value(jj-1)/1.2;
tomwalters@0: 				end
tomwalters@0: 				
tomwalters@0: 			end
tomwalters@0: 		end
tomwalters@0: 		tresval(jj,ii)=current_threshold(jj);
tomwalters@0: 		
tomwalters@0: 		time_since_last_strobe(jj)=current_time-last_strobe_time(jj);
tomwalters@0: 		if time_since_last_strobe(jj) > w_variabel(jj) % linear falling threshold
tomwalters@0: 			const_decay=last_threshold_value(jj)/sr/strobe_decay_time;
tomwalters@0: 			current_threshold(jj)=current_threshold(jj)-const_decay;
tomwalters@0: 			current_threshold(jj)=max(0,current_threshold(jj));
tomwalters@0: 		else    % parabel for the first time y=a(x+b)^2+c
tomwalters@0: 			a=4*(1-h)/(wnull(jj)*wnull(jj));
tomwalters@0: 			b=-wnull(jj)/2;
tomwalters@0: 			c=h;
tomwalters@0: 			factor=a*(time_since_last_strobe(jj) + b) ^2+c;
tomwalters@0: 			current_threshold(jj)=last_threshold_value(jj)*factor;
tomwalters@0: 		end    
tomwalters@0: 		
tomwalters@0: 		current_threshold(jj)=max(0,current_threshold(jj));     % cant be smaller then 0
tomwalters@0: 		last_val(jj)=current_val;
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 	% give back only the strobes, that really occured:
tomwalters@0: for jj=1:nr_channels
tomwalters@0: 	if counter(jj)>1
tomwalters@0: 		real_strobe_points=strobe_points(jj,1:counter(jj)-1);
tomwalters@0: 		strobe_vals=napvals(jj,real_strobe_points);
tomwalters@0: 		all_processes{jj}.strobe_vals=strobe_vals;
tomwalters@0: 		real_strobe_points=bin2time(nap,real_strobe_points);
tomwalters@0: 		all_processes{jj}.strobes=real_strobe_points;
tomwalters@0: 	else
tomwalters@0: 		all_processes{jj}.strobe_vals=[];
tomwalters@0: 		all_processes{jj}.strobes=[];
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: thresholds=nap;
tomwalters@0: thresholds=setvalues(thresholds,tresval);
tomwalters@0: 
tomwalters@0: close(waithand);
tomwalters@0: return