找出两个(非谐波)波之间的相位差 [英] Find phase difference between two (inharmonic) waves

问题描述

我有两个数据集，列出了时间t处两个神经网络组件的平均电压输出，看起来像这样:

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]

当两个神经组件在一定程度上同相"时，这意味着它们是相互关联的.我想做的是计算A和B之间的相位差，最好是在整个仿真过程中.由于两个程序集不可能完全同相，因此我想将该相位差与某个阈值进行比较.

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.

我正在使用numpy和scipy(这两个程序集都是numpy数组)在Python中进行此项目.

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

任何建议将不胜感激！

添加了图

程序集1的示例数据文件

程序集2的示例数据文件

这是两个数据集的外观图:

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

推荐答案

也许您正在寻找互相关:

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.

现在更新，因为我已经查看了真实的数据文件.您发现相移为零有两个原因.首先，两个时间序列之间的相移实际上为零.如果在matplotlib图上水平放大，则可以清楚地看到这一点.其次，重要的是首先对数据进行正则化(最重要的是，减去均值)，否则在阵列末端的零填充效应会使互相关中的实际信号陷入沼泽.在下面的示例中，我通过添加人为偏移并检查是否正确恢复来验证自己是否找到了真实"峰值.

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

编辑:这是一个如何处理假数据的示例.

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

添加了示例图.

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.

这篇关于找出两个(非谐波)波之间的相位差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！