Phase between two 2kHz

phase.py

#!/usr/bin/python

import warnings
warnings.filterwarnings("ignore")

import sys
import numpy as np
import scipy as sy 
from numpy import genfromtxt
import matplotlib.pyplot as plt

# fetch all data from csv file
alldata = genfromtxt(sys.argv[1], delimiter=',')
alldata = alldata[2:] # remove header

# pilot frequency 
pilotfreq = 2000.0
pilotperiod = 500 # in microseconds

# time 
t = alldata[:, 0]
timestep = t[2]-t[1]

sTime = t[-1]-t[0]

print "total time: %f seconds" % sTime

# sampling freq
Fs = 1.0/timestep

print "Sampling freq: %f" % Fs

# original data and response data
ori = alldata[:, 1] 
rsp = alldata[:, 2]

# apply fft
f_ori = np.fft.fft(ori)
f_rsp = np.fft.fft(rsp)

# frequency spectrum (x axis)
freq = np.fft.fftfreq(ori.size, timestep)

# remove phase noise because frequency bins not properly aligned
threshold = max(np.absolute(f_ori))/1000
f_ori_flt = np.array(map(lambda x:complex(0,0) if np.absolute(x) < threshold else x, f_ori))
threshold = max(np.absolute(f_rsp))/1000
f_rsp_flt = np.array(map(lambda x:complex(0,0) if np.absolute(x) < threshold else x, f_rsp))

# calculate angle from fft result (complex number)
f_ori_ang = np.arctan2(f_ori_flt.imag, f_ori_flt.real)*180/np.pi
f_rsp_ang = np.arctan2(f_rsp_flt.imag, f_rsp_flt.real)*180/np.pi

# angle difference
ang_dif =  f_rsp_ang - f_ori_ang

# frequency bins 
bins = Fs/f_ori.size
at = round(bins/pilotfreq)
print "At: %f" % at
print "Angle difference : %f degrees" % ang_dif[at] # angle difference at pilot freq

# from angle difference to time lag at pilot frequency
print "Time delay (approx) : %f us" % ((ang_dif[at] * pilotperiod) / 360.0)

# plot the sine waves
x = np.arange(0, sTime, timestep)
plt.figure(1)
ori, = plt.plot(x, ori[1:], label="original")
rsp, = plt.plot(x, rsp[1:], label="response")
plt.legend(handles=[ori, rsp])

# plot phase angle difference and frequency distribution
plt.figure(2)
fre, = plt.plot(freq, f_ori, 'bo', label="frequency distribution")
ang, = plt.plot(freq, ang_dif, 'ro', label="angle difference")
plt.legend(handles=[fre, ang])
plt.show()

or just use np.angle, so you can np.unwrap it if it's wrapped.