- if plotVerbose or config.PYLAB_VERBOSE :
- common._import_pylab()
- common._pylab.figure(config.BASE_FIGNUM+2)
-
- if deflection_bits != None:
- # plot time series
- common._pylab.subplot(311)
- common._pylab.hold(False)
- common._pylab.plot(deflection_bits, 'r.')
- common._pylab.title("free oscillation")
-
- # plot histogram distribution and gaussian fit
- common._pylab.subplot(312)
- common._pylab.hold(False)
- n, bins, patches = \
- common._pylab.hist(deflection_bits, bins=30,
- normed=1, align='center')
- gaus = numpy.zeros((len(bins),), dtype=numpy.float)
- mean = deflection_bits.mean()
- std = deflection_bits.std()
- pi = numpy.pi
- exp = numpy.exp
- for i in range(len(bins)) :
- gaus[i] = (2*pi)**(-.5)/std * exp(-0.5*((bins[i]-mean)/std)**2)
- common._pylab.hold(True)
- common._pylab.plot(bins, gaus, 'r-');
- common._pylab.hold(False)
-
- # plot FFTed data
- axes = common._pylab.subplot(313)
- else:
- # use a nice big subplot just for the FFTed data
- axes = common._pylab.subplot(111)
- common._pylab.hold(False)
- common._pylab.semilogy(freq_axis, power, 'r.-')
- common._pylab.hold(True)
- xmin,xmax = axes.get_xbound()
- ymin,ymax = axes.get_ybound()
-
- if minFreq is not None and maxFreq is not None:
- # highlight the region we're fitting in
- common._pylab.axvspan(minFreq, maxFreq, facecolor='g', alpha=0.1,
- zorder = -1)
-
- fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2) + D
- common._pylab.plot(freq_axis, fitdata, 'b-')
- noisefloor = D + 0*freq_axis;
- common._pylab.plot(freq_axis, noisefloor)
- common._pylab.hold(False)
- axes.axis([xmin,xmax,ymin,ymax])
-
- common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE) and False: # TODO: cleanup and test
- # write all the ft data now
- fd = file(datafilename, 'w')
- for i in range(len(freq_axis)) :
- fd.write("%g\t%g\n" % (freq_axis[i], power[i]))
- fd.write("\n") # blank line to separate fit data.
- # write the fit data
- for i in range(imin, imax) :
- x = freq_axis[i]
- fd.write("%g\t%g\n" % (freq_axis[i],
- C / ((A**2-x**2)**2 + (x*B)**2) ) )
- fd.close()
- print Vphoto_t.var(), A, B, C, ((numpy.pi * C)/(2 * B * A**2) * 1e-12)
- g("set terminal epslatex color solid")
- g("set output 'calibration.tex'")
- g("set size 2,2") # multipliers on default 5"x3.5"
- #g("set title 'Thermal calibration'")
- g("set logscale y")
- g("set xlabel 'Frequency (Hz)'")
- g("set ylabel 'Power ($\mu$V$^2$/Hz)'")
- g("set xrange [0:20000]")
- g("g(x) = x < %g || x > %g ? 1/0 : f(x)" % (minFreq, maxFreq)) # only plot on fitted region
- g("plot '%s' title 'Vibration data', g(x) title 'Lorentzian fit'" % datafilename)
-
- g("set mouse")
- g("pause mouse 'click with mouse'")
- g.getvar("MOUSE_BUTTON")
- del(g)
-
-# Make vib_analyze splittable (was postponed until the child functions
-# were defined).
-make_splittable_kwargs_function(vib_analyze,
- (fit_psd, 'freq_axis', 'psd_data'),
- (vib_plot, 'deflection_bits', 'freq_axis', 'power', 'A', 'B', 'C', 'D',
- 'minFreq', 'maxFreq'))
-
-@splittableKwargsFunction((vib_analyze, 'deflection_bits', 'freq'))
-def vib_load_analyze_tweaked(tweak_file, analyze=vib_analyze, **kwargs) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- analyze_kwargs, = \
- vib_load_analyze_tweaked._splitargs(vib_load_analyze_tweaked, kwargs)
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- path = path.strip()
- if path[0] == '#': # a comment
- continue
- # read the data
- data = vib_load(path)
- deflection_bits = data['Deflection input']
- freq = float(path.split('_')[-1].split('Hz')[0])
- if config.TEXT_VERBOSE :
- print "Analyzing '%s' at %g Hz" % (path, freq)
- # analyze
- if 'freq' in analyze._kwargs(analyze):
- var = analyze(deflection_bits, freq, **analyze_kwargs)
- else:
- var = analyze(deflection_bits, **analyze_kwargs)
- Vphoto_var.append(var)
- return numpy.array(Vphoto_var, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
- """
- ifile can be a filename string or open (seekable) file object.
- returns an array of data values (1 column)
- """
- #returns (column 1 array, column 2 array)
- #"""
- if ifile == None : ifile = sys.stdin
- unlabeled_data=scipy.io.read_array(ifile)
- return unlabeled_data #(unlabeled_data[:,0], unlabeled_data[:,1])
-
-if __name__ == '__main__' :
- # command line interface
- from optparse import OptionParser
-
- usage_string = ('%prog <input-file>\n'
- '2007-2009 W. Trevor King.\n'
- '\n'
- 'There are several operation modes, each handling a different <input-file>\n'
- 'type. In each case, the output is the fitted variance (in square Volts).\n'
- '\n'
- 'Single file mode without datalogger mode (the default) :\n'
- '<input-file> should be a 1 column ASCII without a header line of cantilever\n'
- 'deflection (in bits)\n'
- 'e.g: "<deflection bits>\\n..."\n'
- '\n'
- 'Single file mode with datalogger mode :\n'
- 'Same as without datalogger mode, except the input should be a datalogger file\n'
- 'with data stored in bits (e.g. as saved by the unfold module).\n'
- '\n'
- 'Tweak file mode:\n'
- 'Runs the same analysis as in single file mode for each vib file in a tweak\n'
- 'file. Each line in the tweak file specifies a single vib file. Blank lines\n'
- 'and those beginning with a pound sign (#) are ignored. Vib files are valid\n'
- 'datalogger files as per single-file-with-datalogger-mode.\n'
- )
- parser = OptionParser(usage=usage_string, version='%prog 0.1')
- parser.add_option('-f', '--sample-freq', dest='freq', default=100e3,
- help='sample frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-m', '--min-freq', dest='min_freq', default=500,
- help='minimum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-M', '--max-freq', dest='max_freq', default=7000,
- help='maximum frequency in Hz (default %default)',
- type='float', metavar='FREQ')
- parser.add_option('-o', '--output-file', dest='ofilename',
- help='write output to FILE (default stdout)',
- type='string', metavar='FILE')
- parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
- help='Output comma-seperated values (default %default)',
- default=False)
- parser.add_option('-g', '--gnuplot', dest='gnuplot', action='store_true',
- help='Print gnuplot fit check script to stderr',
- default=False)
- parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
- help='Produce pylab fit checks during execution',
- default=False)
- parser.add_option('-G', '--plot-guess', dest='plot_guess', action='store_true',
- help='Produce plots of the pre-fit Lorentzian',
- default=False)
- parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
- help='Read input files with datalogger.read_dict_of_arrays(). This is useful, for example, to test a single line from a tweakfile.',
- default=False)
- parser.add_option('-s', '--disable-spectrum-fitting', dest='spectrum_fitting',
- action='store_false',
- help='Disable PSD fitting, just use the raw variance',
- default=True)
- parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
- help='Print lots of debugging information',
- default=False)
-
- options,args = parser.parse_args()
- parser.destroy()
- assert len(args) >= 1, "Need an input file"
-
- ifilename = args[0]
-
- if options.ofilename != None :
- ofile = file(options.ofilename, 'w')
- else :
- ofile = sys.stdout
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = options.gnuplot
- PLOT_GUESSED_LORENTZIAN = options.plot_guess
- if options.spectrum_fitting == True:
- analyze = vib_analyze
- analyze_args = {'freq':options.freq,
- 'minFreq':options.min_freq,
- 'maxFreq':options.max_freq}