'Create Poisson train
'Script to create a timeview file that contains an event channel holding a train of events with Poisson
'statistics, to be used as a nerve stimulus trigger. You need to specify the duration of the train, the mean
'frequency, and the minimal and maximal frequencies. The latter are important for two reasons. The intervals
'are created sequentially in a random process, so very large intervals (small frequencies) make it difficult
'to reach the desired mean frequency because not enough events will subsequently fit into the assigned time
'window. Second, a maximal frequency is necessary because of refractory effects that cause stimulus failures
'when stimulation intervals are too short.
'This script runs fairly slowly, because it creates random event times repeatedly, until the total number of
'events exactly matches the desired mean frequency. The first event will be at 5s, the last at specified train
'duration + 5s. The total number of events is calculated from mean frequency * duration + 1.
'If you need a waveform channel with square pulses, you can create those from the even channel by running the
'"SquareWave" script.
'Dirk Bucher, June 2016, completely replacing an earlier script by Wolfgang Stein.
var ok%,dur,meanF,maxF,minF, nEvts%, title$;
DlgCreate("Assign");
DlgReal(1,"train dur [s]:",1,6000);
DlgReal(2,"mean freq [Hz]:",0.01,1000);
DlgReal(3,"max inst freq [Hz]",0,1000);
DlgReal(4,"min inst freq [Hz]",0,1000);
dur:=300;
meanF:=19;
maxF:=75;
minF:=0.1;
ok%:=DlgShow(dur,meanF,maxF,minF);
nEvts%:=meanF*dur;
title$:="poiss_"+str$(meanF)+"Hz";
var t,n%,random, int, times[200000];
repeat
meanF:=meanF+0.05;
t:=5;
n%:=0;
repeat
random:=Rand();
int:=-1*ln(random)/meanF;
if int<1/minF and int>1/maxF then
t:=t+int;
n%:=n%+1;
times[n%]:=t;
endif
until t>dur+5-1/maxF;
until n%>=nEvts%;
repeat
t:=5;
n%:=0;
repeat
random:=Rand();
int:=-1*ln(random)/meanF;
if int<1/minF and int>1/maxF then
t:=t+int;
n%:=n%+1;
times[n%]:=t;
endif
until t>dur+5-1/maxF;
until n%=nEvts%;
times[0]:=5;
times[nEvts%]:=dur+5;
var k%,mem%;
FileNew(7,1,1,1,dur+10);
mem%:=MemChan(2);
for k%:=0 to nEvts% step 1 do;
MemSetItem(mem%,0,times[k%]);
next
var evtCh%;
evtCh%:=ChanSave(mem%,0);
ChanShow(evtCh%);
ChanDelete(mem%);
XRange(0, MaxTime());
ChanTitle$(evtCh%,title$);
Message(n%+1," events created. Mean frequency = ",n%/dur," Hz");