pythonnumpyneural-networkphysicsscipy# Find phase difference between two (inharmonic) waves

I have two datasets listing the average voltage outputs of two assemblies of neural networks at times t, that look something like this:

```
A = [-80.0, -80.0, -80.0, -80.0, -80.0, -80.0, -79.58, -79.55, -79.08, -78.95, -78.77, -78.45,-77.75, -77.18, -77.08, -77.18, -77.16, -76.6, -76.34, -76.35]
B = [-80.0, -80.0, -80.0, -80.0, -80.0, -80.0, -78.74, -78.65, -78.08, -77.75, -77.31, -76.55, -75.55, -75.18, -75.34, -75.32, -75.43, -74.94, -74.7, -74.68]
```

When two neural assemblies are "in phase" to a reasonable extent, that means that they are interrelated. What I want to do is calculate the phase difference between A and B, preferably over the whole time of the simulation. Since two assemblies are unlikely to be totally in phase, I want to compare that phase difference to a certain threshold.

These are inharmonic oscillators and I don't know their functions, only these values, so I have no idea how to determine the phase or the respective phase difference.

I am doing this project in Python, using numpy and scipy (the two assemblies are numpy arrays).

Any suggestions would be greatly appreciated!

**EDIT: Added plots**

Example datafile for assembly 1

Example datafile for assembly 2

Here is a plot of what the two datasets look like:

Solution

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, scipy
from scipy.signal import correlate
# Load datasets, taking mean of 100 values in each table row
A = 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 datasets
time_shift = 20
A = numpy.roll(A, time_shift)
# Find cross-correlation
xcorr = correlate(A, B)
# delta time array to match xcorr
dt = 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, scipy
from scipy.signal import square, sawtooth, correlate
from numpy import pi, random
period = 1.0 # period of oscillations (seconds)
tmax = 10.0 # length of time series (seconds)
nsamples = 1000
noise_amplitude = 0.6
phase_shift = 0.6*pi # in radians
# construct time array
t = 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 period
B = -sawtooth(phase_shift + 2.0*pi*t/period) + noise_amplitude*random.normal(size=(nsamples,))
# calculate cross correlation of the two signals
xcorr = correlate(A, B)
# The peak of the cross-correlation gives the shift between the two signals
# The xcorr array goes from -nsamples to nsamples
dt = 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-correlation
from pyx import canvas, graph, text, color, style, trafo, unit
from pyx.graph import axis, key
text.set(mode="latex")
text.preamble(r"\usepackage{txfonts}")
figwidth = 12
gkey = 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.

- Python Jinja2 LaTeX Table
- Getting attributes of a class
- How can I print many significant figures in Python?
- How to allow list append() method to return the new list
- Calculate Last Friday of Month in Pandas
- Python type hint for Iterable[str] that isn't str
- How to iterate over a list in chunks
- How to exit the entire application from a Python thread?
- Running shell command and capturing the output
- How do I pass a variable by reference?
- Convert range(r) to list of strings of length 2 in python
- How can I get the start and end dates for each week?
- how to use send_message() in python-telegram-bot
- Python conditional replacement based on element type
- How can I count the number of items in an arbitrary iterable (such as a generator)?
- Find longest consecutive range of numbers in list
- Insert text in braces with asyncpg
- How does one put a link / url to the web-site's home page in Django?
- How to determine if a path is a subdirectory of another?
- Custom Keybindings for Ipython terminal
- FastAPI asynchronous background tasks blocks other requests?
- How to make sure that information from one file is duplicated into several text documents, without specific lines
- Installing a Python environment with Anaconda
- sklearn pipeline model predicting same results for all input
- Brew command not found after installing Anaconda Python
- How to get an XPath from selenium webelement or from lxml?
- Pipe PuTTY console to Python script
- How to align the axes of a figure in matplotlib?
- Persist ParentDocumentRetriever of langchain
- How to reset index in a pandas dataframe?