Find phase difference between two (inharmonic) waves

Perhaps you are looking for the cross-correlation:

scipy.​signal.​signaltools.correlate(A, B)

The position of the peak in the cross-correlation will be an estimate of the phase difference.

EDIT 3: Update now that I have looked at the real data files. There are two reasons that you find a phase shift of zero. First, the phase shift really is zero between your two time series. You can see this clearly if you zoom in horizontally on your matplotlib graph. Second, it is important to regularize the data first (most importantly, subtract off the mean), otherwise the effect of zero-padding at the ends of the arrays swamps the real signal in the cross-correlation. In the following example, I verify that I am finding the "true" peak by adding an artificial shift and then checking that I recover it correctly.

import numpy, scipyfrom scipy.signal import correlate# Load datasets, taking mean of 100 values in each table rowA = numpy.loadtxt("vb-sync-XReport.txt")[:,1:].mean(axis=1)B = numpy.loadtxt("vb-sync-YReport.txt")[:,1:].mean(axis=1)nsamples = A.size# regularize datasets by subtracting mean and dividing by s.d.A -= A.mean(); A /= A.std()B -= B.mean(); B /= B.std()# Put in an artificial time shift between the two datasetstime_shift = 20A = numpy.roll(A, time_shift)# Find cross-correlationxcorr = correlate(A, B)# delta time array to match xcorrdt = numpy.arange(1-nsamples, nsamples)recovered_time_shift = dt[xcorr.argmax()]print "Added time shift: %d" % (time_shift)print "Recovered time shift: %d" % (recovered_time_shift)# SAMPLE OUTPUT:# Added time shift: 20# Recovered time shift: 20

EDIT: Here is an example of how it works with fake data.

EDIT 2: Added a graph of the example.

import numpy, scipyfrom scipy.signal import square, sawtooth, correlatefrom numpy import pi, randomperiod = 1.0                            # period of oscillations (seconds)tmax = 10.0                             # length of time series (seconds)nsamples = 1000noise_amplitude = 0.6phase_shift = 0.6*pi                   # in radians# construct time arrayt = numpy.linspace(0.0, tmax, nsamples, endpoint=False)# Signal A is a square wave (plus some noise)A = square(2.0*pi*t/period) + noise_amplitude*random.normal(size=(nsamples,))# Signal B is a phase-shifted saw wave with the same periodB = -sawtooth(phase_shift + 2.0*pi*t/period) + noise_amplitude*random.normal(size=(nsamples,))# calculate cross correlation of the two signalsxcorr = correlate(A, B)# The peak of the cross-correlation gives the shift between the two signals# The xcorr array goes from -nsamples to nsamplesdt = numpy.linspace(-t[-1], t[-1], 2*nsamples-1)recovered_time_shift = dt[xcorr.argmax()]# force the phase shift to be in [-pi:pi]recovered_phase_shift = 2*pi*(((0.5 + recovered_time_shift/period) % 1.0) - 0.5)relative_error = (recovered_phase_shift - phase_shift)/(2*pi)print "Original phase shift: %.2f pi" % (phase_shift/pi)print "Recovered phase shift: %.2f pi" % (recovered_phase_shift/pi)print "Relative error: %.4f" % (relative_error)# OUTPUT:# Original phase shift: 0.25 pi# Recovered phase shift: 0.24 pi# Relative error: -0.0050# Now graph the signals and the cross-correlationfrom pyx import canvas, graph, text, color, style, trafo, unitfrom pyx.graph import axis, keytext.set(mode="latex")text.preamble(r"\usepackage{txfonts}")figwidth = 12gkey = key.key(pos=None, hpos=0.05, vpos=0.8)xaxis = axis.linear(title=r"Time, \(t\)")yaxis = axis.linear(title="Signal", min=-5, max=17)g = graph.graphxy(width=figwidth, x=xaxis, y=yaxis, key=gkey)plotdata = [graph.data.values(x=t, y=signal+offset, title=label) for label, signal, offset in (r"\(A(t) = \mathrm{square}(2\pi t/T)\)", A, 2.5), (r"\(B(t) = \mathrm{sawtooth}(\phi + 2 \pi t/T)\)", B, -2.5)]linestyles = [style.linestyle.solid, style.linejoin.round, style.linewidth.Thick, color.gradient.Rainbow, color.transparency(0.5)]plotstyles = [graph.style.line(linestyles)]g.plot(plotdata, plotstyles)g.text(10*unit.x_pt, 0.56*figwidth, r"\textbf{Cross correlation of noisy anharmonic signals}")g.text(10*unit.x_pt, 0.33*figwidth, "Phase shift: input \(\phi = %.2f \,\pi\), recovered \(\phi = %.2f \,\pi\)" % (phase_shift/pi, recovered_phase_shift/pi))xxaxis = axis.linear(title=r"Time Lag, \(\Delta t\)", min=-1.5, max=1.5)yyaxis = axis.linear(title=r"\(A(t) \star B(t)\)")gg = graph.graphxy(width=0.2*figwidth, x=xxaxis, y=yyaxis)plotstyles = [graph.style.line(linestyles + [color.rgb(0.2,0.5,0.2)])]gg.plot(graph.data.values(x=dt, y=xcorr), plotstyles)gg.stroke(gg.xgridpath(recovered_time_shift), [style.linewidth.THIck, color.gray(0.5), color.transparency(0.7)])ggtrafos = [trafo.translate(0.75*figwidth, 0.45*figwidth)]g.insert(gg, ggtrafos)g.writePDFfile("so-xcorr-pyx")

So it works pretty well, even for very noisy data and very aharmonic waves.

@deprecated's comments are the exact answer to the question, when it comes to the pure-code python solution. The comments were very valuable, but I feel like I should add some notes for people searching for an answer in the specific context of neural networks.

When you take the average membrane potential of large assemblies of neurons, like I did, the correlation will be relatively weak. What you want to look at, primarily, is either the correlation between spike trains, the latency or the excitability (i.e. synaptic efficacy) of the individual assemblies. This can be found relatively easily by just looking at points where the potential exceeds a certain threshold. Scipy's correlation function on spike trains will show a much more detailed picture of interdependence between neurons or neural assemblies when you give it spike trains, as opposed to the actual potentials. You can also take a look at Brian's statistics module, which can be found here:

http://neuralensemble.org/trac/brian/browser/trunk/brian/tools/statistics.py

As for phase difference, it probably is an inadequate measure, because neurons are not harmonic oscillators. If you want to take very precise measurements of phase, it's best to look at the synchronization of inharmonic oscillators. The mathematical model that describes these kinds of oscillators, which is very useful in the context of neurons and neural networks, is the Kuramoto model. There is extensive documentation available for the Kuramoto model and Integrate-and-fire-synchronization, so I'll leave it to that.