Matlab function
All functions in a single file:
Listing
1 function [phasedrift] = phasefunc(A, B, tau)
2 %
3 % Computes phase function at the specified lag (tau) and returns it.
4 %
5 % Usage:
6 % phasefunc(A, B, tau)
7 %
8 % where A is row matrix holding spike/discrete time events,
9 % B is another such row matrix (A itself for finding auto phase function),
10 % tau is the lag at which the phase function is sought.
11 %
12 % Reference:
13 % A Phase Function to Quantify Serial Dependence between Discrete Samples.
14 % Ramana Dodla and Charles J. Wilson.
15 % Biophys J. 2010 February 17; 98(4): L5-L7.
16 % http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2820633/
17 %
18 % Example: The following code plots a curve in Fig. 2b of the above reference.
19 %
20 % A = [0, 0.18595, 0.33955, 0.4997, 0.61495, 0.7413, 0.9024]
21 % p = [];
22 % for i = 0:400
23 % p=[p phasefunc(A, A, i * 0.001)];
24 % end
25 % x=[0:1:400] * 0.001;
26 % plot(x, p, 'o')
27 %
28 % Author: Ramana Dodla, 2009
29 %
30 [c] = phasedisp_c(A, B + tau);
31 thetadrift = mean(abs(c(:, 1)));
32 norm = sqrt(2.0);
33 phasedrift = thetadrift / norm;
34 end
35 %------- Compute c_i: -----------
36 function [c] = phasedisp_c(A, B)
37 pairs = phasedisp_pairs(A, B);
38 mg = mean(pairs(:, 1));
39 md = mean(pairs(:, 2));
40 f1 = mg * mg + md * md;
41 f2 = sqrt(f1);
42 c = (f1 - pairs(:, 1) * mg - pairs(:, 2) * md) / f2;
43 end
44 %------ Compute (gamma_i, delta_i): --
45 function [pairs] = phasedisp_pairs(A, B)
46 i = 2;
47 j = 2;
48 pairs = [];
49 T = max(A(1), B(1));
50 while i <= length(A) && j <= length(B)
51 if A(i) >= B(j)
52 if A(i-1) >= B(j)
53 T = A(i-1);
54 j = j + 1;
55 else
56 g = (B(j) - T) / (A(i) - A(i-1));
57 d = (B(j) - T) / (B(j) - B(j-1));
58 pairs = [pairs; [g d]];
59 T = B(j);
60 j = j + 1;
61 end;
62 else
63 if B(j-1) >= A(i)
64 T = B(j-1);
65 i = i+1;
66 else
67 g = (A(i) - T)/(A(i) - A(i-1));
68 d = (A(i) - T)/(B(j) - B(j-1));
69 pairs = [pairs; [g d]];
70 T = A(i);
71 i = i + 1;
72 end;
73 end
74 end
75 end