SPI Data Reading and Linear Fit

2014

This is a program that reads in the file media:spi_indistinguishable.txt. A plot media:spi_data.png is created. Then a linear fit is performed. The results of the fit are printed at the end of the program (including errors on parameters), and then the model is plotted with the data in media:spi_data_linear_model.png.

import numpy as np
import pylab as plt

_volts, _A, _B, _C, _AB, _AC, _BC, _ABC = np.loadtxt('spi_indistinguishable.txt',\
                    unpack=True, skiprows=9, usecols=(2,3,4,5,6,7,8,9), \
                    delimiter=',')

print '_volts =', _volts
print '_ABC =', _ABC
print '_AC =', _AC

volts = np.arange(30,71)
print 'volts =', volts
N = len(volts)
AC = np.zeros(N)
AC_err = np.zeros(N)
nsamp = np.zeros(N)

print 'range(N) =', range(N)
for i in range(N):
    print "i =", i
    print "voltage =", volts[i]
    indices = np.where(_volts == volts[i])
    print "indices =", indices
    print "_volts at indices =", _volts[indices]
    print "_AC at indices =", _AC[indices]
    AC[i] = np.mean(_AC[indices])
    AC_err[i] = np.std(_AC[indices])
    nsamp[i] = len(_AC[indices])

plt.clf()
plt.errorbar(volts, AC, AC_err/np.sqrt(nsamp), fmt='r.', label='AC')
plt.xlabel('Piezoelectric Voltage (V)')
plt.ylabel('Count (per 0.01 s)')
plt.title('Indistiguishable Paths')
plt.legend(loc='lower left')
plt.savefig('spi_data.png')

#
# Try linear fit a*cos(2*pi*V/period)+b*sin(2*pi*V/period)+c
#

period=9 #Volts

cos = np.cos(2*np.pi*volts/period)
sin = np.sin(2*np.pi*volts/period)
one = np.ones(len(volts))

MT = np.array([cos, sin, one])  #Transpose of projection matrix
M = np.transpose(MT) #Projection matrix
Ninv = np.diag(1/(AC_err**2/nsamp)) #Inverse noise matrix

# np.dot function produces
D = np.dot(MT,np.dot(Ninv,M))
b = np.dot(MT,np.dot(Ninv,AC))

cov = np.linalg.inv(D)  #Parameter covariance matrix
param = np.dot(cov,b)   #Best fit parameters: a=(D^{-1}M^TN^-1)(data)

print "a =", param[0], '+/-', np.sqrt(cov[0][0])
print "b =", param[1], '+/-', np.sqrt(cov[1][1])
print "c =", param[2], '+/-', np.sqrt(cov[2][2])

plt.plot(volts, param[0]*cos+param[1]*sin+param[2]*one, label='Model')
plt.legend(loc='lower left')
plt.savefig("spi_data_linear_model.png")