123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081 |
- function [T,tvec]=PechePourPoisson(rate,dt)
- % Rate-Varying Poisson-Process generator.
- %
- % T=PechePourPoisson(rate);
- % T=PechePourPoisson(rate,dt);
- % [T tvec]=PechePourPoisson(rate,dt);
- %
- % Example:
- % rate=conv(randn(1,10000),ones(1,1000),'same');
- % [T t]=PechePourPoisson(rate,.001);
- % plot(t,rate);
- % addline(T); addline(0,'h'); % You may not have the addline function
- %
- % This function spits out a poisson process with a rate that can vary with
- % time. Each sample of the "rate" vector is considered to be the beginning
- % of a time bin of width dt. The number of events occurring in each bin is
- % sampled from a Poisson Distribution. The timing of each event within the
- % bin is then randomized.
- %
- % Note: You can easily specify a constant-rate process running for a
- % certain time by just making "rate" a scalar rate and setting dt to your
- % desired time span.
- %
- % Inputs Description
- % "rate" Vector representing the rate as a function of time. Negative
- % rates are interpreted as zero rates.
- % ("dt") Optional scalar defining the time resolution of the rate
- % vector. If neglected, dt=1, and the rate is interpreted as
- % meaning "events per time-sample". If dt is set, rate will be
- % scaled by the value of dt. eg, say you want 90events/s, where
- % the rate vector represents the spiking rate as a function of
- % time with millisecond resolution. Then dt=0.001, and
- % rate(i)=90, where i is the particular instant with a 90event/s
- % rate.
- %
- % Outputs Description
- % "T" Vector of event times. The first sample of "rate" is
- % considered to be time 0. If dt is specified, times will be
- % scaled accordingly.
- % "tvec" If you're too lazy to make your own time vector for your rate
- % vector, this function can optionally do it for you.
- %
- % Enjoy
- % Peter
- % oconnorp _at_ ethz _dot_ ch
- % dt defaults to 1
- if nargin<2, dt=1; end
- % Determine number of events in each time bin
- p=poissrnd(max(rate(:),0)'*dt);
- % Find the times of non-empty bins
- ix=find(p);
- adds=diff([0 ix-1]);
- cum=cumsum([1 p(ix)]);
- % Term is same as sum(p(ix)), just avoids repeating calculation
- T=zeros(cum(end)-1,1);
- T(cum(1:end-1))=adds;
- six=find(~T);
- T=cumsum(T);
- % Sort the randoms (skipping bins with just 1 spike obviously to save time,
- % thought the time taken by the interpreter to parse out this comment is
- % probably even more than that, but I guess its the priniple that counts.)
- T=T+rand(size(T));
- sortable=sort([six(diff([0; six])>1)-1; six]);
- T(sortable)=sort(T(sortable));
- % Efficiency of this last chunk could definitely be improved but really I
- % don't care enough to do it now.
- % Scale times
- if dt~=1, T=T*dt; end
- % For the lazy...
- if nargout>1
- tvec=0:dt:(length(rate)-1)*dt;
- end
- end
|