From 8b9c7258d41e3767170a7fe7a7c86f5744f47434 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Mon, 15 Dec 2008 09:20:40 -0500 Subject: [PATCH] Added comparison to vib_analyze_naive() in vib_analyze(). This will protect against really hideous fits :p. It also lead to fixes for outdated code in vib_analyze_naive() (it hadn't been used in a while) and some outdated comments in vib_analyze(). Finally, added abs()ing for the always-positive fit parameters A and B when using the scipy.optimize.leastsq() method. --- calibcant/vib_analyze.py | 37 ++++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py index b8e8c7c..72c2140 100755 --- a/calibcant/vib_analyze.py +++ b/calibcant/vib_analyze.py @@ -51,10 +51,9 @@ def vib_analyze_naive(deflection_bits, Vphoto_in2V=config.Vphoto_in2V, Vphoto_in2V is a function converting Vphoto values from bits to Volts. This function is assumed linear, see bump_analyze(). - zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts. """ Vphoto_std_bits = deflection_bits.std() - Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zpiezo_zeroV_in) + Vphoto_std = Vphoto_in2V(Vphoto_std_bits + zeroVphoto_bits) return Vphoto_std**2 def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000, @@ -64,18 +63,16 @@ def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000, Calculate the variance in the raw data due to the cantilever's thermal oscillation and convert it to Volts**2. - Improves on vib_analyze_naive() by first converting Vphoto(t) to - frequency space, and fitting a Lorentzian in the relevant frequency - range (see cantilever_calib.pdf for derivation). + Improves on vib_analyze_naive() by first converting Vphoto(t) to + frequency space, and fitting a Lorentzian in the relevant + frequency range (see cantilever_calib.pdf for derivation). + However, there may be cases where the fit is thrown off by noise + spikes in frequency space. To protect from errors, the fitted + variance is compared to the naive variance (where all noise is + included), and the minimum variance is returned. - The conversion to frequency space generates an average power spectrum - by breaking the data into windowed chunks and averaging the power - spectrums for the chunks together. - See FFT_tools.avg_power_spectrum - Vphoto_in2V is a function converting Vphoto values from bits to Volts. This function is assumed linear, see bump_analyze(). - zeroVphoto_bits is the bit value corresponding to Vphoto = zero Volts. """ Vphoto_t = numpy.zeros((len(deflection_bits),), dtype=numpy.float) @@ -102,7 +99,15 @@ def vib_analyze(deflection_bits, freq, minFreq=500, maxFreq=7000, minFreq, maxFreq, plotVerbose=plotVerbose) # Our A is in uV**2, so convert back to Volts**2 - return lorentzian_area(A,B,C) * 1e-12 + fitted_variance = lorentzian_area(A,B,C) * 1e-12 + + naive_variance = vib_analyze_naive(deflection_bits, + Vphoto_in2V=Vphoto_in2V) + if config.TEXT_VERBOSE: + print "Fitted variance: %g V**2", fitted_variance + print "Naive variance : %g V**2", naive_variance + + return min(fitted_variance, naive_variance) def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : """ @@ -112,14 +117,10 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : Calculate the variance in the raw data due to the cantilever's thermal oscillation and convert it to Volts**2. - Improves on vib_analyze_naive() by working on frequency space data - and fitting a Lorentzian in the relevant frequency range (see - cantilever_calib.pdf for derivation). - The conversion to frequency space generates an average power spectrum by breaking the data into windowed chunks and averaging the power spectrums for the chunks together. - See FFT_tools.unitary_avg_power_spectrum(). + See FFT_tools.unitary_avg_power_spectrum(). """ # cut out the relevent frequency range if config.TEXT_VERBOSE : @@ -226,6 +227,8 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : else : print "Solution did not converge" A,B,C = p + A=abs(A) # A and B only show up as squares in f(x) + B=abs(B) # so ensure we get positive values return (A, B, C) def lorentzian_area(A, B, C) : -- 2.26.2