From fea5ed8b55c43b68c258494a0a989a7088df6620 Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Tue, 2 Dec 2008 15:27:44 -0500 Subject: [PATCH] Moved fitting from GnuplotBiDir to scipy.optimize.leastsq in vib_analyze. Much better fit quality. See bug f3ba76cd-eddd-4520-bb2a-17ca371ac6b6 --- .../7bb7514d-f094-461c-8b96-04590415282c/body | 18 +++++ .../values | 28 +++++++ .../d239edab-0769-496e-8c49-fdd25677afb7/body | 28 +++++++ .../values | 21 ++++++ .../d6896169-8e43-4e31-a2c6-53a580ba5ec6/body | 6 ++ .../values | 28 +++++++ .../values | 35 +++++++++ calibcant/vib_analyze.py | 75 ++++++++++++------- 8 files changed, 213 insertions(+), 26 deletions(-) create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/body create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/values create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/body create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/values create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/body create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/values create mode 100644 calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/values diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/body b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/body new file mode 100644 index 0000000..c1c765f --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/body @@ -0,0 +1,18 @@ +I tried resetting and fitting again: + + for i in range(5): + if PLOT_GUESSED_LORENTZIAN: + vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess, + minFreq, maxFreq, plotVerbose=True) + g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess)) + g("fit f(x) '%s' via A,B,C" % datafilename) + A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x) + B=abs(float(g.getvar('B'))) # so ensure we get positive values + C=float(g.getvar('C')) + A_guess = A + B_guess = B + C_guess = C + +But the extra fit cycles had no effect ("After 1 iterations the fit +converged." with no real parameter change.). + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/values b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/values new file mode 100644 index 0000000..93be287 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/7bb7514d-f094-461c-8b96-04590415282c/values @@ -0,0 +1,28 @@ + + + +Content-type=text/plain + + + + + + +Date=Tue, 02 Dec 2008 17:26:04 +0000 + + + + + + +From=W. Trevor King + + + + + + +In-reply-to=d239edab-0769-496e-8c49-fdd25677afb7 + + + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/body b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/body new file mode 100644 index 0000000..189e4e6 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/body @@ -0,0 +1,28 @@ +Consider + +$ cc_vib_analyze.py -m 300 -M 7000 -f50000 -p -G -d 20080919/20080919121041_zp_time_50000Hz + +The guess looks reasonable. It is a bit off due to some popcorn +noise, but it looks like it should fit. The fit, however, looks bad. +It is definately too flat. Looking at the git.log gnuplot output, we see + + After 10 iterations the fit converged. + final sum of squares of residuals : 2.58573e+09 + rel. change during last iteration : -6.04806e-06 + ... + Final set of parameters Asymptotic Standard Error + ======================= ========================== + + A = 5115 +/- 184.1 (3.599%) + B = 5476.03 +/- 433.1 (7.909%) + C = 2.73854e+18 +/- 4.604e+17 (16.81%) + ... + correlation matrix of the fit parameters: + + A B C + A 1.000 + B 0.905 1.000 + C 0.970 0.963 1.000 + +The high correlation matrix can't be helping... + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/values b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/values new file mode 100644 index 0000000..9de4f63 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d239edab-0769-496e-8c49-fdd25677afb7/values @@ -0,0 +1,21 @@ + + + +Content-type=text/plain + + + + + + +Date=Tue, 02 Dec 2008 17:22:52 +0000 + + + + + + +From=W. Trevor King + + + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/body b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/body new file mode 100644 index 0000000..74c0d70 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/body @@ -0,0 +1,6 @@ +Moving from GnuplotBiDir fitting to scipy.optimize.leastsq fitting. + +Major improvement in fit quality on our problem case: + +$ cc_vib_analyze.py -m 300 -M 7000 -f50000 -p -G -d 20080919/20080919121041_zp_time_50000Hz + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/values b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/values new file mode 100644 index 0000000..32f0604 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/comments/d6896169-8e43-4e31-a2c6-53a580ba5ec6/values @@ -0,0 +1,28 @@ + + + +Content-type=text/plain + + + + + + +Date=Tue, 02 Dec 2008 20:26:45 +0000 + + + + + + +From=W. Trevor King + + + + + + +In-reply-to=7bb7514d-f094-461c-8b96-04590415282c + + + diff --git a/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/values b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/values new file mode 100644 index 0000000..3f2da80 --- /dev/null +++ b/calibcant/.be/bugs/f3ba76cd-eddd-4520-bb2a-17ca371ac6b6/values @@ -0,0 +1,35 @@ + + + +creator=W. Trevor King + + + + + + +severity=minor + + + + + + +status=fixed + + + + + + +summary=Poor fit for 20080919121041_zp_time_50000Hz + + + + + + +time=Tue, 02 Dec 2008 17:08:41 +0000 + + + diff --git a/calibcant/vib_analyze.py b/calibcant/vib_analyze.py index 283f1ae..b8e8c7c 100755 --- a/calibcant/vib_analyze.py +++ b/calibcant/vib_analyze.py @@ -33,7 +33,8 @@ The relevent physical quantities are : """ import os, time, numpy -import GnuplotBiDir # used for fitting lorentzian +import GnuplotBiDir # used for fitting lorentzian +import scipy.optimize # alternative for fitting lorentzian import common # common module for the calibcant package import config # config module for the calibcant package import data_logger @@ -130,13 +131,6 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : assert imax >= imin + 10 , \ "Less than 10 points in freq range (%g,%g)" % (minFreq, maxFreq) - # save about-to-be-fitted data to a temporary file - datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time()) - fd = file(datafilename, 'w') - for i in range(imin, imax) : - fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i])) - fd.close() - # generate guesses for Lorentzian parameters A,B,C max_psd_index = numpy.argmax(psd_data[imin:imax]) + imin max_psd = psd_data[max_psd_index] @@ -178,31 +172,60 @@ def fit_psd(freq_axis, psd_data, minFreq=500, maxFreq=7000) : # (a-w)**2 = (2A-B)/2 * [1 +/- sqrt(1 - 4A(A-2B)/(B-2A)**2] # so w is a disaster ;) # For some values of A and B (non-underdamped), W is imaginary or negative. + # + # The mass m is given by m = k_B T / (pi A) + # The spring constant k is given by k = m (omega_0)**2 + # The quality factor Q is given by Q = omega_0 m / gamma if config.TEXT_VERBOSE : print "guesses : %g, %g, %g" % (A_guess, B_guess, C_guess) if PLOT_GUESSED_LORENTZIAN: vib_plot(None, freq_axis, psd_data, A_guess, B_guess, C_guess, minFreq, maxFreq, plotVerbose=True) - # fit Lorentzian using Gnuplot's 'fit' command - g = GnuplotBiDir.Gnuplot() - # The Lorentzian is the predicted one-sided power spectral density - # For an oscillator whose position x obeys - # m x'' + gamma x' + kx = F_thermal(t) - # A is the ?driving noise? - # B is gamma, the damping term - # C is omega_0, the resonant frequency # Fitting the PSD of V = photoSensitivity*x just rescales the parameters - g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)") - ## The mass m is given by m = k_B T / (pi A) - ## The spring constant k is given by k = m (omega_0)**2 - ## The quality factor Q is given by Q = omega_0 m / gamma - g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess)) - g("fit f(x) '%s' via A,B,C" % datafilename) - A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x) - B=abs(float(g.getvar('B'))) # so ensure we get positive values - C=float(g.getvar('C')) - os.remove(datafilename) + + if False: # Gnuplot fit worse than scipy.optimize.leastsq fit. + # fit Lorentzian using Gnuplot's 'fit' command + + # save about-to-be-fitted data to a temporary file + datafilename = "%s_%d" % (config.GNUFIT_DATA_BASE, time.time()) + fd = file(datafilename, 'w') + for i in range(imin, imax) : + fd.write("%g\t%g\n" % (freq_axis[i], psd_data[i])) + fd.close() + + g = GnuplotBiDir.Gnuplot() + g("f(x) = C / ((A**2-x**2)**2 + (B*x)**2)") + g("A=%g; B=%g; C=%g" % (A_guess, B_guess, C_guess)) + g("fit f(x) '%s' via A,B,C" % datafilename) + A=abs(float(g.getvar('A'))) # A and B only show up as squares in f(x) + B=abs(float(g.getvar('B'))) # so ensure we get positive values + C=float(g.getvar('C')) + + os.remove(datafilename) + else: + # fit Lorenzian using scipy.optimize.leastsq + def residual(p, y, x): + A = p[0] + B = p[1] + C = p[2] + Y = C / ((A**2-x**2)**2 + (B*x)**2) + return Y-y + leastsq = scipy.optimize.leastsq + p,cov,info,mesg,ier = leastsq(residual, (A_guess, B_guess, C_guess), + args=(psd_data[imin:imax], + freq_axis[imin:imax]), + full_output=True, maxfev=10000) + if config.TEXT_VERBOSE : + print "Fitted params:",p + print "Covariance mx:",cov + print "Info:", info + print "mesg:", mesg + if ier == 1 : + print "Solution converged" + else : + print "Solution did not converge" + A,B,C = p return (A, B, C) def lorentzian_area(A, B, C) : -- 2.26.2