Eye Diagram
The code below generates the following plot:
The main script generates num_traces traces, and on a grid of 600x600, it counts the number times a trace crosses a grid point. The grid is then plotted using matplotlib's imshow() function. The counting is performed using Bresenham's line algorithm (http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm), to ensure that the counting is correct, and steep parts of the curve don't result in missed counts.
Bresenham's algorithm is slow in pure Python, so a Cython version is included. If you do not build the Cython version of the Bresenham code, be sure to reduce num_traces before running the program!
Here's the main demo script, eye_demo.py.
1 import numpy as np
2
3 use_fast = True
4 try:
5 from brescount import bres_curve_count
6 except ImportError:
7 print "The cython version of the curve counter is not available."
8 use_fast = False
9
10
11 def bres_segment_count_slow(x0, y0, x1, y1, grid):
12 """Bresenham's algorithm.
13
14 The value of grid[x,y] is incremented for each x,y
15 in the line from (x0,y0) up to but not including (x1, y1).
16 """
17
18 nrows, ncols = grid.shape
19
20 dx = abs(x1 - x0)
21 dy = abs(y1 - y0)
22
23 sx = 0
24 if x0 < x1:
25 sx = 1
26 else:
27 sx = -1
28 sy = 0
29 if y0 < y1:
30 sy = 1
31 else:
32 sy = -1
33
34 err = dx - dy
35
36 while True:
37 # Note: this test is moved before setting
38 # the value, so we don't set the last point.
39 if x0 == x1 and y0 == y1:
40 break
41
42 if 0 <= x0 < nrows and 0 <= y0 < ncols:
43 grid[x0, y0] += 1
44
45 e2 = 2 * err
46 if e2 > -dy:
47 err -= dy
48 x0 += sx
49 if e2 < dx:
50 err += dx
51 y0 += sy
52
53 def bres_curve_count_slow(x, y, grid):
54 for k in range(x.size - 1):
55 x0 = x[k]
56 y0 = y[k]
57 x1 = x[k+1]
58 y1 = y[k+1]
59 bres_segment_count_slow(x0, y0, x1, y1, grid)
60
61
62 def random_trace(t):
63 s = 2*(np.random.randint(0, 5) % 2) - 1
64 r = 0.01 * np.random.randn()
65 s += r
66 a = 2.0 + 0.001 * np.random.randn()
67 q = 2*(np.random.randint(0, 7) % 2) - 1
68 t2 = t + q*(6 + 0.01*np.random.randn())
69 t2 += 0.05*np.random.randn()*t
70 y = a * (np.exp(s*t2) / (1 + np.exp(s*t2)) - 0.5) + 0.07*np.random.randn()
71 return y
72
73
74 if __name__ == "__main__":
75 import matplotlib.pyplot as plt
76 grid_size = 600
77 grid = np.zeros((grid_size, grid_size), dtype=np.int32)
78
79 tmin = -10.0
80 tmax = 10.0
81 n = 81
82 t = np.linspace(tmin, tmax, n)
83 dt = (tmax - tmin) / (n - 1)
84
85 ymin = -1.5
86 ymax = 1.5
87
88 num_traces = 100000
89
90 for k in range(num_traces):
91
92 # Add some noise to the times at which the signal
93 # will be sampled. Without this, all the samples occur
94 # at the same times, and this produces an aliasing
95 # effect in the resulting bin counts.
96 # If n == grid_size, this can be dropped, and t2 = t
97 # can be used instead. (Or, implement an antialiased
98 # version of bres_curve_count.)
99 steps = dt + np.sqrt(0.01 * dt) * np.random.randn(n)
100 steps[0] = 0
101 steps_sum = steps.cumsum()
102 t2 = tmin + (tmax - tmin) * steps_sum / steps_sum[-1]
103
104 td = (((t2 - tmin) / (tmax - tmin)) * grid_size).astype(np.int32)
105
106 y = random_trace(t2)
107
108 # Convert y to integers in the range [0,grid_size).
109 yd = (((y - ymin) / (ymax - ymin)) * grid_size).astype(np.int32)
110
111 if use_fast:
112 bres_curve_count(td, yd, grid)
113 else:
114 bres_curve_count_slow(td, yd, grid)
115
116 plt.figure()
117 # Convert to float32 so we can use nan instead of 0.
118 grid = grid.astype(np.float32)
119 grid[grid==0] = np.nan
120 plt.grid(color='w')
121 plt.imshow(grid.T[::-1,:], extent=[0,1,0,1], cmap=plt.cm.coolwarm,
122 interpolation='gaussian')
123 ax = plt.gca()
124 ax.set_axis_bgcolor('k')
125 ax.set_xticks(np.linspace(0,1,11))
126 ax.set_yticks(np.linspace(0,1,11))
127 ax.set_xticklabels([])
128 ax.set_yticklabels([])
129 plt.colorbar()
130 fig = plt.gcf()
131
132 #plt.savefig("eye-diagram.png", bbox_inches='tight')
133 plt.show()
Here's brescount.pyx, the Cython implementation of Bresenham's line algorithm:
1 import numpy as np
2 cimport numpy as np
3 cimport cython
4
5
6 @cython.boundscheck(False)
7 cdef int bres_segment_count(unsigned x0, unsigned y0,
8 unsigned x1, unsigned y1,
9 np.ndarray[np.int32_t, ndim=2] grid):
10 """Bresenham's algorithm.
11
12 See http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
13 """
14
15 cdef unsigned nrows, ncols
16 cdef int e2, sx, sy, err
17 cdef int dx, dy
18
19 nrows = grid.shape[0]
20 ncols = grid.shape[1]
21
22 if x1 > x0:
23 dx = x1 - x0
24 else:
25 dx = x0 - x1
26 if y1 > y0:
27 dy = y1 - y0
28 else:
29 dy = y0 - y1
30
31 sx = 0
32 if x0 < x1:
33 sx = 1
34 else:
35 sx = -1
36 sy = 0
37 if y0 < y1:
38 sy = 1
39 else:
40 sy = -1
41
42 err = dx - dy
43
44 while True:
45 # Note: this test occurs before increment the
46 # grid value, so we don't count the last point.
47 if x0 == x1 and y0 == y1:
48 break
49
50 if (x0 < nrows) and (y0 < ncols):
51 grid[x0, y0] += 1
52
53 e2 = 2 * err
54 if e2 > -dy:
55 err -= dy
56 x0 += sx
57 if e2 < dx:
58 err += dx
59 y0 += sy
60
61 return 0
62
63
64 def bres_curve_count(np.ndarray[np.int32_t, ndim=1] x,
65 np.ndarray[np.int32_t, ndim=1] y,
66 np.ndarray[np.int32_t, ndim=2] grid):
67 cdef unsigned k
68 cdef int x0, y0, x1, y1
69
70 for k in range(len(x)-1):
71 x0 = x[k]
72 y0 = y[k]
73 x1 = x[k+1]
74 y1 = y[k+1]
75 bres_segment_count(x0, y0, x1, y1, grid)
76 if 0 <= x1 < grid.shape[0] and 0 <= y1 < grid.shape[1]:
77 grid[x1, y1] += 1
This file, setup.py, can be used to build the Cython extension module:
To build the extension module, you must have Cython installed.
You can build the extension module as follows:
$ python setup.py build_ext --inplace