calibcant/psd_filter_analyze.py

   1 #!/usr/bin/python
   2
   3 if __name__ == "__main__" :
   4
   5     import sys
   6     import filter, calibcant_vib_analyze
   7     from optparse import OptionParser
   8
   9     usage_string = ('%prog <input-file>\n'
  10                     '2008, W. Trevor King.\n'
  11                     '\n'
  12                     'Apply various filters to data y(x).'
  13                     '<input-file> should be whitespace-delimited, 2 column ASCII\n'
  14                     'without a header line.\n'
  15                     'e.g: "<x>\\t<y>\\n"')
  16     parser = OptionParser(usage=usage_string, version='%prog 0.1')
  17     parser.add_option('-o', '--output-file', dest='ofilename',
  18                       help='write output to FILE (default stdout)',
  19                       type='string', metavar='FILE')
  20     parser.add_option('-w', '--width', dest='width', default=None,
  21                       help='window width (i +/- width, inclusive)',
  22                       type='int', metavar='WIDTH')
  23     parser.add_option('-t', '--type', dest='type', default='mean',
  24                       help='filter type (default %default)',
  25                       type='string', metavar='TYPE')
  26     parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
  27                       help='minimum frequency in Hz (default %default)',
  28                       type='float', metavar='FREQ')
  29     parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
  30                       help='maximum frequency in Hz (default %default)',
  31                       type='float', metavar='FREQ')
  32     parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
  33                       help='Print gnuplot fit check script to stderr',
  34                       default=False)
  35     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
  36                       help='Print lots of debugging information',
  37                       default=False)
  38
  39     options,args = parser.parse_args()
  40     parser.destroy()
  41     assert len(args) >= 1, "Need an input file"
  42
  43     ifilename = args[0]
  44
  45     # mirror filter behavior:
  46     freq,psd = filter.read_data(ifilename)
  47
  48     if options.type == 'mean' :
  49         winfunc = filter.windowed_mean_point
  50     elif options.type == 'median' :
  51         winfunc = filter.windowed_median_point
  52     else :
  53         raise Exception, "unrecognized window type '%s'" % options.type
  54
  55     psd = filter.windowed_filter(psd, winfunc, width=options.width)
  56
  57     # mirror calibcant_vib_analyze behavior:
  58     (A,B,C) = calibcant_vib_analyze.fit_psd(freq, psd,
  59                                             minFreq=options.min_freq,
  60                                             maxFreq=options.max_freq)
  61     if options.verbose :
  62         print >> sys.stderr, "Fit f(x) = C / ((A**2-x**2)**2 + (x*B)**2) with"
  63         print >> sys.stderr, "A = %g, \t B = %g, \t C = %g" % (A, B, C)
  64     if options.gnuplot :
  65         print >> sys.stderr, \
  66             calibcant_vib_analyze.gnuplot_check_fit(ifilename, A, B, C)
  67
  68     area = calibcant_vib_analyze.lorentzian_area(A,B,C)
  69
  70     if options.ofilename != None :
  71         print >> file(options.ofilename, 'w'), area
  72     else :
  73         print area