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: % % other parameters
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,'threshold')
tomwalters@0:     threshold=params.threshold;
tomwalters@0: else
tomwalters@0:     threshold=0; % in dB
tomwalters@0: end
tomwalters@0: if isfield(params,'dynamic_range')
tomwalters@0:     dynamic_range=params.dynamic_range;
tomwalters@0: else
tomwalters@0:     dynamic_range=20; % in dB
tomwalters@0: end
tomwalters@0: if isfield(params,'jitter_time')
tomwalters@0:     jitter_time=params.jitter_time;
tomwalters@0: else
tomwalters@0:     jitter_time=0.1;
tomwalters@0: end
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,'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: % spont_rate=0.015;
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=1/getsr(sig)*1000;
tomwalters@0: sig_len=getlength(sig);
tomwalters@0: nr_steps=getnrpoints(sig);
tomwalters@0: 
tomwalters@0: 
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: % test_sig=getvalues(sig);
tomwalters@0: % test_sig=test_sig+60; % shift it upwards
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: %             difft=time_now-last_spike;    % how long ago is the last spike?
tomwalters@0: %             if difft<1, continue,  end
tomwalters@0: %
tomwalters@0: % %             % modify eta by the amplitude:
tomwalters@0: % %             cur_amp=test_sig(j)-threshold; % in dB above threshold
tomwalters@0: % %             if cur_amp>dynamic_range % in saturation
tomwalters@0: % %                 cur_eta=eta;
tomwalters@0: % %             else
tomwalters@0: % %                 cur_eta=(eta-eta_null)/dynamic_range*cur_amp+eta_null;
tomwalters@0: % % %                 cur_eta=f2f(cur_amp,0,dynamic_range,eta_null,eta);
tomwalters@0: % %             end
tomwalters@0: %
tomwalters@0: % %             z = (log(difft) - p1) ./ p2;
tomwalters@0: % %             spike_prob = exp(p3 - exp(p3)) ./ p2;
tomwalters@0: %             spike_prob = gevpdf(log(difft),p3,p2,p1);
tomwalters@0: %
tomwalters@0: % %             weibull function
tomwalters@0: % %             spike_prob=(beta/cur_eta)*power(log(difft)/cur_eta,beta-1);
tomwalters@0: %
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: %                 last_spike=time_now; % 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: %         end
tomwalters@0: % %     end
tomwalters@0: %     ret_spikes{i}=spikes;
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: %% the version with 2 different prob functions
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: sigvalues=getvalues(sig);
tomwalters@0: x2=getxvalues(sig).*1000;
tomwalters@0: x2=x2(x2>0);
tomwalters@0: x2=[x2; x2(end)+binwidth];
tomwalters@0: % figure(43)
tomwalters@0: % cla
tomwalters@0: % hold on
tomwalters@0: for i=1:2
tomwalters@0:     p1=params.recovery_parameter.fit_p(1);
tomwalters@0:     p2=params.recovery_parameter.fit_p(2);
tomwalters@0:     p3=params.recovery_parameter.fit_p(3);
tomwalters@0:     if i==1
tomwalters@0:         p2=p2/2;
tomwalters@0:     end
tomwalters@0:     pdf=gevpdf(log(x2),p1,p2,p3);
tomwalters@0:     cdf=gevcdf(log(x2),p1,p2,p3);
tomwalters@0:     spike_prob_function{i}=ones(size(x2))*inf;
tomwalters@0:     for j=1:length(x2)
tomwalters@0:         if cdf(j)<1
tomwalters@0:             spike_prob_function{i}(j)=pdf(j)/(1-cdf(j))/params.acc_fac;
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: times=[20 250];
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % spike_prob_function2=gevpdf(x2,p3,p2/start_boost_beta_val,p1); % for the first 20 ms
tomwalters@0: 
tomwalters@0: % shift the signal by the latency to simulate the latency:
tomwalters@0: zeros_in_front=zeros(latency/binwidth,1);
tomwalters@0: sigvalues=[zeros_in_front;sigvalues];
tomwalters@0: 
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:     swapc=1;
tomwalters@0:     next_swap=times(swapc);
tomwalters@0:     c_function=spike_prob_function{swapc};
tomwalters@0: 
tomwalters@0:     still_in_latency=1;
tomwalters@0:     for j=1:nr_steps
tomwalters@0:         time_now=j*binwidth; % thats the global time counter
tomwalters@0: 
tomwalters@0:         if still_in_latency && time_now> latency
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; % save and add the latency
tomwalters@0:                 spikecounter=spikecounter+1; %remember the spike, when it happens
tomwalters@0:                 still_in_latency=0;
tomwalters@0:             end
tomwalters@0: 
tomwalters@0: 
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=round(difft/binwidth);    sindx=max(1,sindx);    sindx=min(350,sindx);
tomwalters@0: 
tomwalters@0:             spike_prob=c_function(sindx);
tomwalters@0:             %timefound=find(difft<times,1,'first');
tomwalters@0:             %spike_prob=spike_prob_function{timefound}(sindx);
tomwalters@0: 
tomwalters@0:             spike_prob=spike_prob*sigvalues(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; % 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
tomwalters@0: 
tomwalters@0: return
tomwalters@0: % 
tomwalters@0: % %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: % %% the version with 10 different prob functions
tomwalters@0: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: % 
tomwalters@0: % x2=getxvalues(sig).*1000;
tomwalters@0: % x2=x2(x2>0);
tomwalters@0: % x2=[x2; x2(end)+binwidth];
tomwalters@0: % % figure(43)
tomwalters@0: % % cla
tomwalters@0: % % hold on
tomwalters@0: % for i=10:-1:1
tomwalters@0: %     if isfield(p{i},'k') % otherwise no spikes
tomwalters@0: %         p1=p{i}.k;
tomwalters@0: %         p2=log(p{i}.x);
tomwalters@0: %         p3=p{i}.y;
tomwalters@0: %         %         plot3(p1,p2,p3,'o')
tomwalters@0: %         pdf=gevpdf(log(x2),p1,p2,p3);
tomwalters@0: %         cdf=gevcdf(log(x2),p1,p2,p3);
tomwalters@0: %         spike_prob_function{i}=ones(size(x2))*inf;
tomwalters@0: %         for j=1:length(x2)
tomwalters@0: %             if cdf(j)<1
tomwalters@0: %                 spike_prob_function{i}(j)=pdf(j)/(1-cdf(j))/params.acc_fac;
tomwalters@0: %             end
tomwalters@0: %         end
tomwalters@0: %     else
tomwalters@0: %         % otherwise take the "default" one
tomwalters@0: %         spike_prob_function{i}=spike_prob_function{10};
tomwalters@0: %     end
tomwalters@0: % end
tomwalters@0: % 
tomwalters@0: % times=0:10:60;
tomwalters@0: % times=[times 100 250];
tomwalters@0: % 
tomwalters@0: % 
tomwalters@0: % % spike_prob_function2=gevpdf(x2,p3,p2/start_boost_beta_val,p1); % for the first 20 ms
tomwalters@0: % 
tomwalters@0: % jitter_time=0;
tomwalters@0: % latency=0;
tomwalters@0: % 
tomwalters@0: % start_boost_beta_length=20;
tomwalters@0: % test_sig=zeros(length(x2),1);
tomwalters@0: % sigend=53;
tomwalters@0: % test_sig(1:sigend/binwidth)=1;
tomwalters@0: % % build a ramp...
tomwalters@0: % test_sig(1:2/binwidth)=linspace(0,1,2/binwidth);
tomwalters@0: % test_sig(sigend/binwidth-2/binwidth:sigend/binwidth)=linspace(1,0,2/binwidth+1);
tomwalters@0: % 
tomwalters@0: % spont_rate=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: %     swapc=2;
tomwalters@0: %     next_swap=times(swapc);
tomwalters@0: %     c_function=spike_prob_function{swapc};
tomwalters@0: % 
tomwalters@0: %     for j=1:nr_steps
tomwalters@0: %         time_now=j*binwidth; % thats the global time counter
tomwalters@0: %         if time_now>next_swap
tomwalters@0: %             swapc=swapc+1;
tomwalters@0: %             next_swap=times(swapc);
tomwalters@0: %             c_function=spike_prob_function{swapc};
tomwalters@0: %         end
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: % 
tomwalters@0: %             % follow the first spike prob
tomwalters@0: %             spike_prob=spike_prob_function{1}(j);
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: %                 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: % 
tomwalters@0: % 
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=round(difft/binwidth);    sindx=max(1,sindx);    sindx=min(350,sindx);
tomwalters@0: % 
tomwalters@0: %             spike_prob=c_function(sindx);
tomwalters@0: %             %timefound=find(difft<times,1,'first');
tomwalters@0: %             %spike_prob=spike_prob_function{timefound}(sindx);
tomwalters@0: % 
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
tomwalters@0: %