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

问题描述

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

A = [-80.0,-80.0,-79.58,-79.55,-79.08,-78.95,-78.77,-78.45,-77.75,-77.18,-77.08,-77.16,-76.6,-76.34,-76.35]

B = [-80.0,-78.74,-78.65,-78.08,-77.31,-76.55,-75.55,-75.18,-75.34,-75.32,-75.43,-74.94,-74.7,-74.68]

当两个神经组件在一定程度上“同相”时，这意味着它们是相互关联的。我想做的是计算A和B之间的相位差，最好是在整个仿真过程中。由于两个程序集不可能完全同相，因此我想将该相位差与某个阈值进行比较。 这些是非谐波振荡器，我不知道它们的功能，只有这些值，所以我不知道如何确定相位或各自的相位差。 我正在使用numpy和scipy（两个程序集是numpy数组）在Python中进行此项目。 任何建议将不胜感激！ 编辑：添加图 程序集1的示例数据文件 程序集2的示例数据文件 这是两个数据集的外观图：     

解决方法

        也许您正在寻找互相关：

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

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

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\"\\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\")

因此，即使对于非常嘈杂的数据和非谐波波，它也能很好地工作。     ,        当涉及纯代码python解决方案时，@ deprecated \的注释是该问题的确切答案。这些评论非常有价值，但是我觉得我应该为在神经网络特定上下文中寻找答案的人们添加一些注释。 当像我一样获取大型神经元组件的平均膜电位时，相关性将相对较弱。您要查看的主要是各个序列的峰值序列，潜伏时间或兴奋性（即突触功效）之间的相关性。只需查看电势超过某个阈值的点，就可以相对容易地找到它。当您给尖峰序列提供Scipy时，Scipy的相关函数将显示神经元或神经组件之间相互依存的详细信息，而不是实际的电势。您还可以查看Brian的统计信息模块，可以在这里找到： http://neuralensemble.org/trac/brian/browser/trunk/brian/tools/statistics.py 至于相位差，这可能是一个不足的度量，因为神经元不是谐波振荡器。如果要非常精确地测量相位，则最好查看非谐波振荡器的同步。仓本模型是描述这类振荡器的数学模型，它在神经元和神经网络的环境中非常有用。有大量有关Kuramoto模型和“集成并发射同步”的文档，因此我将其保留。