- common._pylab.subplot(312)
- 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
- common._pylab.subplot(313)
- common._pylab.semilogy(freq_axis, power, 'r.-')
- fitdata = C / ((A**2-freq_axis**2)**2 + (B*freq_axis)**2)
- common._pylab.hold(True)
- common._pylab.plot(freq_axis, fitdata, 'b-');
- common._pylab.hold(False)
-
- common._flush_plot()
- if (plotVerbose or config.GNUPLOT_VERBOSE): # 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)
-
-def vib_load_analyze_tweaked(tweak_file, minFreq=1000, maxFreq=7000,
- chunk_size=2048, overlap=True,
- window=FFT_tools.window_hann,
- Vphoto_in2V=config.Vphoto_in2V,
- plotVerbose=False) :
- """
- Read the vib files listed in tweak_file.
- Return an array of Vphoto variances (due to the cantilever) in Volts**2
- """
- dl = data_logger.data_load()
- Vphoto_var = []
- for path in file(tweak_file, 'r') :
- if path[-1] == '\n':
- path = path.split('\n')[0]
- # 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
- Vphoto_var.append(vib_analyze(deflection_bits, freq, minFreq, maxFreq,
- chunk_size, overlap, window,
- Vphoto_in2V, plotVerbose))
- 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'
- '2008, W. Trevor King.\n'
- '\n'
- 'Deflection power spectral density (PSD) data (X^2/Hz)\n'
- 'and returns the variance in X units (X^2)'
- '<input-file> should be whitespace-delimited, 2 column ASCII\n'
- 'without a header line.\n'
- 'e.g: "<frequency_Hz>\\t<deflection_psd X^2/Hz>\\n"\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('-t', '--tweak-mode', dest='tweakmode', action='store_true',
- help='Run in tweak-file mode',
- default=False)
- 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
- if options.verbose == True :
- vfile = sys.stderr
- else :
- vfile = None
- config.TEXT_VERBOSE = options.verbose
- config.PYLAB_INTERACTIVE = False
- config.PYLAB_VERBOSE = options.pylab
- config.GNUPLOT_VERBOSE = options.gnuplot
-
- if options.tweakmode == False :
- data = read_data(ifilename)
- deflection_bits = data['Deflection input']
- Vphoto_var = vib_analyze(deflection_bits, freq=options.freq,
- minFreq=options.min_freq,
- maxFreq=options.max_freq)
- print >> ofile, Vphoto_var
- else : # tweak mode
- Vphoto_vars = vib_load_analyze_tweaked(ifilename,
- minFreq=options.min_freq,
- maxFreq=options.max_freq)
- if options.comma_out :
- sep = ','
- else :
- sep = '\n'
- common.write_array(ofile, Vphoto_vars, sep)
-
- if common._final_flush_plot != None:
- common._final_flush_plot()
-
- if options.ofilename != None :
- ofile.close()