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:

   1 from distutils.core import setup
   2 from distutils.extension import Extension
   3 from Cython.Distutils import build_ext
   4 
   5 import numpy
   6 
   7 ext = Extension("brescount", ["brescount.pyx"],
   8     include_dirs = [numpy.get_include()])
   9 
  10 setup(ext_modules=[ext],
  11       cmdclass = {'build_ext': build_ext})

To build the extension module, you must have Cython installed.

You can build the extension module as follows:

$ python setup.py build_ext --inplace

Cookbook/EyeDiagram (last edited 2012-02-04 22:49:55 by WarrenWeckesser)