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: 
tomwalters@0: 
tomwalters@0: function ret=gen_twoDat2003(bmm,options)
tomwalters@0: % two D adaptive thresholding
tomwalters@0: 
tomwalters@0: % first do the normal nap:
tomwalters@0: nap=gen_hcl(bmm,options.nap);
tomwalters@0: 
tomwalters@0: if nargin < 2
tomwalters@0:     options=[];
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: vals=getvalues(nap);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%normalise to 0-60dB range
tomwalters@0: %max_nap=max(max(vals));
tomwalters@0: %nap_norm=60./max_nap;
tomwalters@0: %vals=vals*nap_norm;
tomwalters@0: %save vals.mat vals;
tomwalters@0: 
tomwalters@0: new_vals=zeros(size(vals));
tomwalters@0: nr_chan=size(vals,1);
tomwalters@0: nr_dots=size(vals,2);
tomwalters@0: sr=getsr(nap);
tomwalters@0: cfs=getcf(nap);
tomwalters@0: %load bmm.mat;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % the delays can be set by an option
tomwalters@0: if isfield(options,'time_constant_factor') % this is multiplied to the threshold_time_const
tomwalters@0:     time_constant_factor=options.time_constant_factor;
tomwalters@0: else
tomwalters@0:     time_constant_factor=0.8;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if isfield(options,'frequency_constant_factor') % the influence of the left and the right channel
tomwalters@0:     frequency_constant_factor=options.frequency_constant_factor;
tomwalters@0: else
tomwalters@0:     frequency_constant_factor=0.9;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if isfield(options,'threshold_rise_constant') % the rate at which the threshold may rise
tomwalters@0:     threshold_rise_constant=options.threshold_rise_constant;
tomwalters@0: else
tomwalters@0:     threshold_rise_constant=1;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %the following calculates the envelope of nap response
tomwalters@0: %out_li=zeros(nr_chan,nr_dots);
tomwalters@0: %for chan=1:nr_chan
tomwalters@0: %out_li(chan,:)=standard_leaky_integrator(vals(chan,:),0.003,1,16000);
tomwalters@0: %out_li_scale=max(vals(chan,:))./max(out_li(chan,:));
tomwalters@0: %out_li(chan,:)=out_li(chan,:).*out_li_scale;
tomwalters@0: %end
tomwalters@0: %[ERB_no, dummy] = Freq2ERB(cfs);
tomwalters@0: %for chan=1:nr_chan
tomwalters@0: %out_li(chan,:)=standard_leaky_integrator(vals(chan,:),0.003,1,16000);
tomwalters@0: %out_li_scale=max(vals(chan,:))./max(out_li(chan,:));
tomwalters@0: %out_li(chan,:)=out_li(chan,:).*out_li_scale;
tomwalters@0: %end
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %save cfs.mat cfs;
tomwalters@0: order = 4;
tomwalters@0: [rate ERB] = Freq2ERB(cfs);
tomwalters@0: B=options.b *ERB;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %note: decay of gammatone is described by: exp(-2*pi*B*t)
tomwalters@0: % time to drop to 1/2 (6dB) --in seconds
tomwalters@0: %i.e. exp(-2*pi*B*t)=0.5
tomwalters@0: drop6dbtime=-log(0.5)./(2*pi.*B);
tomwalters@0: 
tomwalters@0: % this means so many per samplepoint
tomwalters@0: threshold_decay_constant=-log10(2)*20./(drop6dbtime*sr);
tomwalters@0: threshold_decay_constant=threshold_decay_constant./time_constant_factor;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %threshold_rise_constant=.5;
tomwalters@0: 
tomwalters@0: threshold_rise=(100.*B.^0.9)./sr; 
tomwalters@0: threshold_rise=threshold_rise.*threshold_rise_constant;
tomwalters@0: %threshold_rise=zeros(nr_chan);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: %calc frequency slope values
tomwalters@0: if nr_chan>1
tomwalters@0:     %frequency_slope_c=calcfreqslope(cfs(2),cfs(1));
tomwalters@0:     %frequency_slope_a=calcfreqslope(cfs(1),cfs(2));
tomwalters@0:     
tomwalters@0:     %c is effect from channel above but we must input
tomwalters@0:     %the higher frequency (2) as the adjacent band and
tomwalters@0:     %the channel in which we want the level (1)as 
tomwalters@0:     %the centre freq
tomwalters@0:     frequency_slope_c=calcfreqslope(cfs(1),cfs(2),options.b);
tomwalters@0:     
tomwalters@0:     %a is effect from channel below but we must input
tomwalters@0:     %the lower frequency (1) as the adjacent band and
tomwalters@0:     %the channel in which we want the level (2) as 
tomwalters@0:     %the centre freq
tomwalters@0:     frequency_slope_a=calcfreqslope(cfs(2),cfs(1),options.b);
tomwalters@0:     
tomwalters@0: else
tomwalters@0:     frequency_slope_c=0;
tomwalters@0:     frequency_slope_a=0;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: %we divide by the frequency_constant_factor because we want a
tomwalters@0: %frequency_constant_factor of 1 to mean we use the actual
tomwalters@0: %frequency slope and a frequency_constant_factor of 0.5
tomwalters@0: %to mean that the frequency slope falls away at twice the rate
tomwalters@0: %i.e. less effect
tomwalters@0: frequency_slope_c=frequency_slope_c./frequency_constant_factor;
tomwalters@0: frequency_slope_a=frequency_slope_a./frequency_constant_factor;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: oldfigure=gcf;
tomwalters@0: if nr_chan==1
tomwalters@0: 	onechannelfigure=figure(3);	% assuming, the other one is 1
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % channel_select=1:nr_chan;
tomwalters@0: times_per_ms=round(sr*0.01);
tomwalters@0: if nr_chan>1
tomwalters@0: 	waithand=waitbar(0,'generating 2dat');
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %val_e=zeros(nr_chan,nr_dots);
tomwalters@0: %val_inp=zeros(nr_chan,nr_dots);
tomwalters@0: %val_thres=zeros(nr_chan,nr_dots);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % a is the decayed working variable from the channel above
tomwalters@0: % b is the range limit
tomwalters@0: % c is the decayed working variable from the channel below
tomwalters@0: % d is the delayed and decayed nap signal from on-channel
tomwalters@0: % e is the maximum of [a b c d]
tomwalters@0: % threshold is the working variable
tomwalters@0: % thresholdlast is the working variable in the last round
tomwalters@0: % inp is the current input
tomwalters@0: 
tomwalters@0: 
tomwalters@0: thresholdlast=zeros(nr_chan,1);
tomwalters@0: threshold=zeros(nr_chan,1);
tomwalters@0: a=0;b=1;c=0;d=0;e=0;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: for ii=1:nr_dots % through the time
tomwalters@0: 	if nr_chan>1
tomwalters@0: 		if mod(ii,times_per_ms)==0
tomwalters@0: 			waitbar(ii/nr_dots);
tomwalters@0: 		end
tomwalters@0: 	end
tomwalters@0:     
tomwalters@0:     for jj=1:nr_chan  % through all channels: prepare working variable
tomwalters@0:         inp=vals(jj,ii);% current input
tomwalters@0:         
tomwalters@0:         if jj< nr_chan          
tomwalters@0:             chan_above=thresholdlast(jj+1);
tomwalters@0:             a=chan_above+frequency_slope_a;
tomwalters@0:         else
tomwalters@0:             a=0;
tomwalters@0:         end
tomwalters@0:         if jj>1
tomwalters@0:             chan_below=thresholdlast(jj-1);
tomwalters@0:             c=chan_below+frequency_slope_c;
tomwalters@0:         else
tomwalters@0:             c=0;
tomwalters@0:         end
tomwalters@0:              
tomwalters@0:         %a=0;
tomwalters@0:         %c=0;
tomwalters@0:         
tomwalters@0:        
tomwalters@0:         
tomwalters@0:         
tomwalters@0:         d=max(thresholdlast(jj)+threshold_decay_constant(jj),0);
tomwalters@0:         
tomwalters@0:         
tomwalters@0:         e1=max(a,b);   
tomwalters@0:         e2=max(c,d);
tomwalters@0:         e=max(e1,e2);   
tomwalters@0:         
tomwalters@0:         
tomwalters@0:                
tomwalters@0:         
tomwalters@0:         if inp>e   %threshold is exceeded
tomwalters@0:             
tomwalters@0:            threshold(jj)=min(thresholdlast(jj)+threshold_rise(jj),inp);
tomwalters@0:            %threshold(jj)=inp; 
tomwalters@0:            
tomwalters@0:             %new_vals(jj,ii)=threshold(jj)-e;
tomwalters@0:            new_vals(jj,ii)=inp-threshold(jj);
tomwalters@0:             %    % else threshold decays
tomwalters@0:             else
tomwalters@0:               
tomwalters@0:             threshold(jj)=e; 
tomwalters@0:             new_vals(jj,ii)=0;  
tomwalters@0:         end
tomwalters@0:     
tomwalters@0:              
tomwalters@0:         
tomwalters@0:          %val_e(jj,ii)=e;
tomwalters@0:          %val_inp(jj,ii)=inp;
tomwalters@0:          %val_thres(jj,ii)=threshold(jj);
tomwalters@0:         
tomwalters@0:                                
tomwalters@0:     end
tomwalters@0:     
tomwalters@0:     
tomwalters@0:     
tomwalters@0:     
tomwalters@0:     
tomwalters@0:     % prepare next round 
tomwalters@0:     if ii< nr_dots                   
tomwalters@0:     thresholdlast=threshold;
tomwalters@0:     end
tomwalters@0:     
tomwalters@0:     if nr_chan==1
tomwalters@0:         tvals(ii)=thresholdlast;
tomwalters@0:     end
tomwalters@0:     
tomwalters@0:     
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if nr_chan==1
tomwalters@0:     figure(onechannelfigure);
tomwalters@0:     plot(nap);
tomwalters@0:     hold on
tomwalters@0:     a=plot(tvals,'r');
tomwalters@0:     set(a,'LineWidth',1.5);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: nap=setvalues(nap,new_vals);
tomwalters@0: %save nap.mat nap;
tomwalters@0: 
tomwalters@0: vals=getvalues(nap);
tomwalters@0: nap=setallmaxvalue(nap,max(max(vals))*2);
tomwalters@0: nap=setallminvalue(nap,min(min(vals)));
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %save val_e.mat val_e;
tomwalters@0: %save new_vals.mat new_vals;
tomwalters@0: %save val_inp.mat val_inp;
tomwalters@0: %save val_thres.mat val_thres;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: ret=nap;
tomwalters@0: if nr_chan>1
tomwalters@0: 	close(waithand);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: figure(oldfigure);