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.

diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py
index b8e8c7c3397e50ebab2af581a660c64e14a0efb9..72c2140fecf6c9183408065506082e402c6b943b 100755 (executable)
--- 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) :