bleeck@3: % This external file is included as part of the 'aim-mat' distribution package
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
tomwalters@0: function ret_spikes=sim_spikes(sig,nr_sweeps,params)
tomwalters@0: % simulates the response to a signal
tomwalters@0: 
tomwalters@0: % parameters of the recovery function
tomwalters@0: mu=params.mu;
tomwalters@0: beta=params.beta;
tomwalters@0: 
tomwalters@0: % other parameters
tomwalters@0: if isfield(params,'spont_rate')
tomwalters@0:     spont_rate=params.spont_rate;
tomwalters@0: else
tomwalters@0:     spont_rate=0;
tomwalters@0: end
tomwalters@0: if isfield(params,'latency')
tomwalters@0:     latency=params.latency;
tomwalters@0: else
tomwalters@0:     latency=4;
tomwalters@0: end
tomwalters@0: % if isfield(params,'first_spike_boost')
tomwalters@0: %     first_spike_boost=params.first_spike_boost;
tomwalters@0: % else
tomwalters@0: %     first_spike_boost=2;
tomwalters@0: % end
tomwalters@0: % if isfield(params,'start_boost_beta_length')
tomwalters@0: %     start_boost_beta_length=params.start_boost_beta_length;
tomwalters@0: % else
tomwalters@0: %     start_boost_beta_length=15;
tomwalters@0: % end
tomwalters@0: % if isfield(params,'start_boost_beta_val')
tomwalters@0: %     start_boost_beta_val=params.start_boost_beta_val;
tomwalters@0: % else
tomwalters@0: %     start_boost_beta_val=2;
tomwalters@0: % end
tomwalters@0: if isfield(params,'jitter')
tomwalters@0:     jitter_time=params.jitter;
tomwalters@0: else
tomwalters@0:     jitter_time=0.2;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % calcualates an artificail spiketrain with so many sweeps
tomwalters@0: % the important parameter are mu (the location) and beta (the scale) of the
tomwalters@0: % distribution
tomwalters@0: 
tomwalters@0: % changing some parameters so that the result looks good
tomwalters@0: %jitter to reduce some effect that the odd (or even?) multiples
tomwalters@0: %of the sr have significant and repeatable higher nr intervals
tomwalters@0: %then the others. Cant work out why. solution: let the interval
tomwalters@0: %jitter around the peak a bit
tomwalters@0: 
tomwalters@0: % calculate the x values of the distribution
tomwalters@0: % take a length of the distribution up to the point where it it 99.9 %
tomwalters@0: % that saves a lot of time
tomwalters@0: maxpoint=0.999;
tomwalters@0: tlength=mu-beta*log(-log(maxpoint));
tomwalters@0: 
tomwalters@0: % set the random number generator to a new value
tomwalters@0: % rand('state', sum(100*clock));
tomwalters@0: % dont do that because the values are too close to each other and actually
tomwalters@0: % not random at all!
tomwalters@0: 
tomwalters@0: binwidth=getsr(sig);
tomwalters@0: sig_len=getlength(sig);
tomwalters@0: nr_steps=sig_len/binwidth;
tomwalters@0: x2=binwidth:binwidth:tlength;
tomwalters@0: 
tomwalters@0: spike_prob_function=hazard_function(x2,mu,beta); % for the full signal length
tomwalters@0: spike_prob_function2=hazard_function(x2,mu,beta/start_boost_beta_val); % for the first 20 ms
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % build a "test signal" out of zeros and ones
tomwalters@0: test_sig=zeros(nr_steps,1);
tomwalters@0: test_sig(1:50/sig_sr)=1;  % the first 50 ms are ones
tomwalters@0: % build a ramp...
tomwalters@0: test_sig(1:2/sig_sr)=linspace(0,1,length(5/sig_sr));
tomwalters@0: 
tomwalters@0: 
tomwalters@0: for i=1:nr_sweeps
tomwalters@0:     spikes=[];
tomwalters@0:     last_spike=-1; % inital condition: how long ago was the last spike
tomwalters@0:     spikecounter=1; % spikecounter    
tomwalters@0:     for j=1:nr_steps
tomwalters@0:         time_now=j*binwidth; % thats the global time counter
tomwalters@0:         % implementation of a simple solution for the first spike problem: if the spike is the first then assume a very high probability
tomwalters@0:         if spikecounter==1 % yes, its the first
tomwalters@0:             if rand<binwidth*first_spike_boost % if a random number is smaller, then ...
tomwalters@0:                 jitter=randfloat(-jitter_time,jitter_time);
tomwalters@0:                 last_spike=time_now+jitter; % yes, a spike has happend now!
tomwalters@0:                 spikes(spikecounter)=last_spike+latency; % save and add the latency
tomwalters@0:                 spikecounter=spikecounter+1; %remember the spike, when it happens
tomwalters@0:             end
tomwalters@0:         else % its a follow up spike
tomwalters@0:             difft=time_now-last_spike;    % how long ago is the last spike?
tomwalters@0:             sindx=ceil(difft/binwidth);    sindx=max(1,sindx);    sindx=min(length(spike_prob_function),sindx);
tomwalters@0:             if time_now<start_boost_beta_length
tomwalters@0:                 spike_prob=spike_prob_function2(sindx);
tomwalters@0:             else
tomwalters@0:                 spike_prob=spike_prob_function(sindx);
tomwalters@0:             end
tomwalters@0:             spike_prob=spike_prob*test_sig(j)+spont_rate; %    % modulate the probability with the height of the signal
tomwalters@0:             spike_prob=spike_prob/(1.2+binwidth/2); %correction factor
tomwalters@0:             if rand<spike_prob % if a random number is smaller, then ...
tomwalters@0:                 jitter=randfloat(-jitter_time,jitter_time);
tomwalters@0:                 last_spike=time_now+jitter; % yes, a spike has happend now!
tomwalters@0:                 % make sure that it is not too close to the last one (as a result of the jitter)
tomwalters@0:                 if last_spike<spikes(spikecounter-1)+0.1;
tomwalters@0:                     last_spike=time_now;
tomwalters@0:                 end
tomwalters@0:                 spikes(spikecounter)=last_spike+latency; % save and add the latency
tomwalters@0:                 spikecounter=spikecounter+1; %remember the spike, when it happens
tomwalters@0:             end
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0:     ret_spikes{i}=spikes;
tomwalters@0: end