找出两个（非谐波）波之间的相位差

Posted 2023-03-12

【中文标题】找出两个（非谐波）波之间的相位差

【英文标题】：Find phase difference between two (inharmonic) waves

【发布时间】：2011-09-03 16:59:02

【问题描述】：

我有两个数据集，列出了两个神经网络组件在时间 t 的平均电压输出，如下所示：

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 之间的相位差，最好是在整个模拟过程中。由于两个组件不太可能完全同相，因此我想将该相位差与某个阈值进行比较。

这些是非谐波振荡器，我不知道它们的功能，只知道这些值，所以我不知道如何确定相位或各自的相位差。

我在 Python 中做这个项目，使用 numpy 和 scipy（这两个程序集是 numpy 数组）。

任何建议将不胜感激！

编辑：添加地块

Example datafile for assembly 1

Example datafile for assembly 2

这是两个数据集的样子：

【问题讨论】：

请检查您的数据文件。他们充满了“80”！

【参考方案1】：

也许您正在寻找互相关：

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

互相关中峰值的位置将是相位差的估计值。

编辑 3： 现在更新，我已经查看了真实的数据文件。您发现相移为零有两个原因。首先，两个时间序列之间的相移确实为零。如果您在 matplotlib 图形上水平放大，您可以清楚地看到这一点。其次，重要的是首先对数据进行正则化（最重要的是，减去均值），否则数组末端的零填充效应会淹没互相关中的真实信号。在以下示例中，我通过添加人为移位然后检查我是否正确恢复它来验证我是否找到了“真实”峰值。

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

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

编辑 2：添加了示例图表。

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"\usepackagetxfonts")
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) = \mathrmsquare(2\pi t/T)\)", A, 2.5), (r"\(B(t) = \mathrmsawtooth(\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"\textbfCross 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")

所以它工作得很好，即使是非常嘈杂的数据和非常非谐波的波。

【讨论】：

我已经尝试过了，即 numpy.argmax(signal.correlate(listX,listY)))，但它只返回列表中元素的数量...

奇数。如果你绘制它，自相关数组是什么样子的？

顺便说一下，“非谐波”是指波形是严格周期性的而不是正弦的吗？或者它们不是严格周期性的？而且，这两个序列的周期是否相同？

@Dow - 抱歉，我不明白您发布的数据文件。每个包含 20,000 行，每行由一个整数（行号）和 100 个浮点数组成。除了第一行，每行中绝大多数（>95%）的浮点数都是相同的值（-80.0）。我们应该如何从这些文件中提取您的时间序列 A 和 B？

@Dow - 你的两个时间序列绝对是“同相”的，因为 A 中的峰值与 B 中的峰值重合。一个和另一个之间的滞后为零。我还没有检查它们的周期性。如果您对此感兴趣，可以寻找功率谱中的峰值 - 例如numpy.abs(numpy.fft.fft(A))**2。

【参考方案2】：

@deprecated 的 cmets 是对纯代码 python 解决方案的确切答案。 cmets 非常有价值，但我觉得我应该为在神经网络的特定上下文中寻找答案的人们添加一些注释。

当您像我一样采用大型神经元集合的平均膜电位时，相关性会相对较弱。您首先要查看的是脉冲序列之间的相关性、延迟或单个组件的兴奋性（即突触功效）。只需查看潜力超过某个阈值的点，就可以相对容易地找到这一点。当你给它脉冲序列时，Scipy 的脉冲序列相关函数将显示神经元或神经组件之间相互依赖的更详细的图片，而不是实际的电位。你也可以看看 Brian 的统计模块，可以在这里找到：

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

至于相位差，它可能是一个不充分的衡量标准，因为神经元不是谐振子。如果您想对相位进行非常精确的测量，最好查看非谐波振荡器的同步。描述这类振荡器的数学模型在神经元和神经网络的背景下非常有用，是仓本模型。对于 Kuramoto 模型和 Integrate-and-fire-synchronization 有大量可用的文档，所以我将把它留给那些。