calibcant/psd_filter_analyze.py

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