var spikecount,totalampl, logfilename$,vh%,n,BINSZ,datapoinz,STaDEV,STaDEV$,above,below,above2, below2,above3, below3, Spikes, delay,eventzahlabove,eventzahlbelow,x;

' this script find action potentials in an extracellular recording according to Farzan Nadims idea that one can
' use the standard deviation of a trace to determine signal / noise ratio

' it calculates the standard deviation of a channel and then uses every minimum and maximum between 1 x stddev
' and 10 x stddev to define spikes. If a maximum is followed by a minimum within 5ms, a spike is counted.
' the script writes spike times and spike minimum and maximum values as well as total amplitude (min -> max)
' to the logfile for further processing.


vh%:=View(); 'get view of data file

'clear log window
View(LogHandle());
		EditSelectAll();
		EditClear();

'bring data file to front again

View(vh%);


' make new cursors
Cursordelete(-1);
CursorNew();
Cursornew();

' set region to analyze
interact("set Cursors around region to analyze",1023);

' make toolbar
ToolbarClear();
ToolbarSet(1,"find spikes",findpeak%);
ToolbarSet(2, "save Logfile", dothat%);
ToolbarSet(3, "Stop", stop%);

repeat
  Toolbar("",1023);
until 0;
Halt;

Halt;


' find the spikes 
func findpeak%();

' get channel to look for spikes
n:=input("What channel?",2);

' get minimum delay between spikes
delay:=input("minimum time between spikes in ms?",5);
delay:=delay/1000;

' calculate and state standard deviation
STaDEV:=Chanmeasure(n,12,Cursor(1),Cursor(2));
STaDEV$:=("the standarddeviation is "+Str$(STaDEV));
message(STaDEV$);


'using 1 x STDDEV as minimum and 10xSTDDEV as maximum
' define and create necessary memory channels
below:=memchan(3,0);
below2:=memchan(3,0);
below3:=memchan(3,0);
above:=memchan(3,0);
above2:=memchan(3,0);
above3:=memchan(3,0);
Chanshow(below);
Chanshow(below2);
Chanshow(below3);
Chanshow(above);
Chanshow(above2);
Chanshow(above3);

spikes:=MemChan(3,0);'holds spikes
Chanshow(spikes);

'these are the channels for rising and falling through min and max values
Memimport(below,n,Cursor(1),Cursor(2),3,delay,(-1)*STaDEV);
'Memimport(below2,n,Cursor(1),Cursor(2),3,delay,(-10)*STaDEV);
Memimport(above,n,Cursor(1),Cursor(2),2,delay,STaDEV);
'Memimport(above2,n,Cursor(1),Cursor(2),2,delay,10*STaDEV);

'these are the channels for FALLING /Rising through treshold 1 
Memimport(below3,n,Cursor(1),Cursor(2),2,delay,(-1)*STaDEV);
Memimport(above3,n,Cursor(1),Cursor(2),3,delay,STaDEV);

'window discriminator: All spikes below 1*STDDEV and all above 10* STDEV are omitted. Same for negative values

var punkt1a; 'variable defining the rising time of the spike
var punkt1b; 'variable defining the falling time of the spike

var punkt1c; 'variable defining the falling time of the spike
var punkt1d; 'variable defining the rising time of the spike

var peakamp1; ' variable defining the peak voltage of the spike
var peakampl2;


punkt1a:=Cursor(1);
punkt1c:=Cursor(1);

while punkt1a>0 do 
	punkt1a:=nexttime(above,punkt1a); 'rising time of spike
	punkt1b:=nexttime(above3,punkt1a); 'falling time of spike
if punkt1a>0 and punkt1b>0 then 
peakamp1:=ChanMeasure(n,8,punkt1a,punkt1b);  'amplitude of spike
endif;

if peakamp1<10*STaDEV then 'look for time of AP peak and write to channel

MemSetItem(above2,0,punkt1a);
'MemImport(eventCh1%, inputCh1%, punkt1a, punkt1b, 0, 0, 0); 'finds the time for the peak voltage
endif;

wend;  



while punkt1c>0 do 
	punkt1c:=nexttime(below,punkt1c); 'rising time of spike
	punkt1d:=nexttime(below3,punkt1c); 'falling time of spike
if punkt1c>0 and punkt1d>0 then 
peakampl2:=ChanMeasure(n,8,punkt1c,punkt1d);  'amplitude of spike
endif;

if peakampl2>-10*STaDEV then 'look for time of AP peak and write to channel

MemSetItem(below2,0,punkt1c);
'MemImport(eventCh1%, inputCh1%, punkt1a, punkt1b, 0, 0, 0); 'finds the time for the peak voltage
endif;

wend;  

'that got data that is above STD DEV and below -1*STDDEV, above and below zero. above2 and below2 are the channels that hold these spikes
'now deleting unnecessary channels:

Chandelete(above);
Chandelete(below);
Chandelete(above3);
Chandelete(below3);


'Now: find spikes that go up and down!
spikecount:=0;
punkt1a:=Cursor(1);
while punkt1a>0 do 
	punkt1a:=nexttime(above2,punkt1a); 'next rising spike
	punkt1b:=nexttime(below2,punkt1a); 'next falling spike
if punkt1a>0 and punkt1b>0 and punkt1b-punkt1a<0.005 then  'assumes 5 msec spike duration!!!!!!!!!!!!!!!!!!!!!
 
peakamp1:=ChanMeasure(n,8,punkt1a,punkt1b);  'positive amplitude of spike
peakampl2:=Chanmeasure(n,9,punkt1b,punkt1b+delay); 'negative amplitude of spike
totalampl:=ChanMeasure(n,8,punkt1a,punkt1b)-Chanmeasure(n,9,punkt1b,punkt1b+delay);
spikecount:=spikecount+1;

'write to new channel

memsetitem(spikes, 0, punkt1a);
'CHECK SIGN OF AMPLITUDES!!!!!!!!!!!!!!!!!!!!!!
Printlog("Spike_number ", spikecount, " Time: ",punkt1a, " pos.amplitude: ",peakamp1," total_amplitude ",totalampl);
Printlog("Spike_number ", spikecount, " Time: ",punkt1b, " neg.amplitude: ",peakampl2," total_amplitude ",totalampl);

endif;


wend;  


eventzahlbelow:=Count(above,0,maxtime());
eventzahlbelow:=Count(below,0,maxtime());


'BinSZ:=Binsize(n);
'datapoinz:=(Cursor(2)-Cursor(1))/BinSZ;
'datapoinz:=round(datapoinz);
'var all[datapoinz];
'all[]:=Chanmeasure(n,14,Cursor(1), Cursor(2));
'Chandata(n,all[],Cursor(1), Cursor(2));

'Printlog(all[]);

return;
end;

func dothat%();

logfilename$:=input$("give logfile name","here.txt");
' and save the file
View(LogHandle()); ' show it
FileSaveAs(logfilename$,1); ' save it

return;
end;

func stop%()


halt;
end;