function HP = MalvestioPRE2017Hpointprocess(spiketimes,Q,q,s,W,W_ex,dchoice) %% MalvestioPRE2017Hpointprocess v 1.0 02.10.2018 % Adaptation of AndrzejakKreuzDpointprocess with the implemention of spike % train distance made by Eero Satuvuori (http://wwwold.fi.isc.cnr.it/users/thomas.kreuz/Source-Code/cSPIKE.html ). % It calculates the distances between every pairs of windows of the same spike train, and from this the corresponding % matrix of the ranks. % In particular, HP(i,j) = g_{i,j}, introduced in the paper [1] before Eq. 1. % % spiketimes: 1D vector containing times of spikes % Q: interval of time during which two dynamics are simultaneously measured. % q: breakdown of Q into smaller overlaping intervals (q < Q). % s: step size of the overlaping intervals q (s < q). % s and q need to be choose in order to have (Q-q)/s as an integer(otherwise the very end of the signal may be ignored). % W: Theiler Window, to avoid comparison between neighboring windows % W_ex: number of windows excluded at the very beginning and at the very end. % dchoice allow to select the distance: % 1 : ISI distance % 2 : SPIKE distance % 3 : Adaptive ISI distance % 4 : Adaptive SPIKE distance % % For further information, please refer to % MalvestioPRE2017Readme.pdf % and to [1]: % Malvestio I, Kreuz T, Andrzejak RG: % Robustness and versatility of a nonlinear interdependence method for directional coupling detection from spike trains. % PRE, 96.2 (2017) 022203. InitializecSPIKE; N=fix((Q-q)/(s)+1); DP=zeros(N,N); if dchoice == 3 || 4 duration_total = Q; spikes_total1 = spiketimes(spiketimes>0 & spiketimes< duration_total); spikemat_total{1} = spikes_total1; spikemat_total{2} = spikes_total1; STS_total = SpikeTrainSet(spikemat_total, 0, duration_total); [threshold, ~] = giveTHR(STS_total); end for wc1=W+1:N-1 % calculate distance profile winspikes1=spiketimes(spiketimes>0 & spiketimes0+wc1*s & spiketimes