calibcant/filter.py

   1 #!/usr/bin/python
   2
   3 """
   4 windowed_filter(array, winfunc, s)
   5 with winfuncs :
   6   windowed_mean_point
   7   windowed_median_point
   8
   9 Command line bindings for filtering 2 column ACSII data.
  10 """
  11
  12 import numpy
  13
  14 def windowed_mean_point(array, i, width=50) :
  15     """
  16     Filter data with a windowed average (mean).
  17     The ith point is replaced by the average of the points with indices in
  18     the range i +/- width (inclusive).
  19
  20     This function returns the ith point in the filtered array,
  21     given the raw array and window half-width.
  22
  23     The weights on the average are uniform at the moment.
  24     """
  25     x=0
  26     n=0
  27     for j in range(i-width, i+width) :
  28         if j >= 0 and j < len(array) :
  29             x+=array[j]
  30             n+=1
  31     x /= float(n)
  32     return x
  33
  34 def windowed_median_point(array, i, width=9) :
  35     """
  36     Filter data with a windowed median.
  37     The ith point is replaced by the median of the points with indices in
  38     the range i +/- width (inclusive).
  39
  40     This function returns the ith point in the filtered array,
  41     given the raw array and window half-width.
  42     """
  43     imin = i-width
  44     if imin < 0 : imin = 0
  45     imax = i+width
  46     if imax >= len(array) : imax = len(array-1)
  47     slice = array[imin:imax+1].copy()
  48     slice.sort()
  49     imid = numpy.floor((imax-imin)/2)
  50     return slice[imid]
  51
  52 def windowed_filter(array, winfunc, width=None) :
  53     """
  54     Filter data with a windowing function winfunc.
  55     The ith point is replaced by the winfunc(array, i, s).
  56
  57     See the windowed_* functions for possible winfunc options.
  58     """
  59     out = array.copy()
  60     if width == None :
  61         win = lambda i : winfunc(array, i) # user winfunc's default s
  62     else :
  63         win = lambda i : winfunc(array, i, width)
  64     for i in range(len(out)) :
  65         out[i] = win(i)
  66     return out
  67
  68
  69 # commandline interface functions
  70 import scipy.io, sys
  71
  72 def read_data(ifile):
  73     """
  74     ifile can be a filename string or open (seekable) file object.
  75     returns (column 1 array, column 2 array)
  76     """
  77     if ifile == None :  ifile = sys.stdin
  78     data=scipy.io.read_array(ifile)
  79     return (data[:,0], data[:,1])
  80
  81 def write_data(ofile, x, y):
  82     """
  83     ofile can be a filename string or open (seekable) file object.
  84     """
  85     data = numpy.zeros((len(x),2))
  86     data[:,0] = x
  87     data[:,1] = y
  88     if ofile == None :  ofile = sys.stdout
  89     scipy.io.write_array(ofile, data, separator='\t')
  90
  91 if __name__ == '__main__' :
  92     # command line interface
  93     from optparse import OptionParser
  94
  95     usage_string = ('%prog <input-file>\n'
  96                     '2008, W. Trevor King.\n'
  97                     '\n'
  98                     'Apply various filters to data y(x).'
  99                     '<input-file> should be whitespace-delimited, 2 column ASCII\n'
 100                     'without a header line.\n'
 101                     'e.g: "<x>\\t<y>\\n"')
 102     parser = OptionParser(usage=usage_string, version='%prog 0.1')
 103     parser.add_option('-o', '--output-file', dest='ofilename',
 104                       help='write output to FILE (default stdout)',
 105                       type='string', metavar='FILE')
 106     parser.add_option('-w', '--width', dest='width', default=None,
 107                       help='window width (i +/- width, inclusive)',
 108                       type='int', metavar='WIDTH')
 109     parser.add_option('-t', '--type', dest='type', default='mean',
 110                       help='filter type (default %default)',
 111                       type='string', metavar='TYPE')
 112     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
 113                       help='Print lots of debugging information',
 114                       default=False)
 115
 116     options,args = parser.parse_args()
 117     parser.destroy()
 118     assert len(args) >= 1, "Need an input file"
 119
 120     ifilename = args[0]
 121
 122     x,y = read_data(ifilename)
 123
 124     if options.type == 'mean' :
 125         winfunc = windowed_mean_point
 126     elif options.type == 'median' :
 127         winfunc = windowed_median_point
 128     else :
 129         raise Exception, "unrecognized window type '%s'" % options.type
 130
 131     y = windowed_filter(y, winfunc, width=options.width)
 132
 133     write_data(options.ofilename, x, y)