Trouver la différence de phase entre les deux (inharmoniques) des ondes

Peut-être vous êtes à la recherche pour le cross-corrélation:

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

La position du pic de corrélation croisée est une estimation de la différence de phase.

EDIT 3: mettre à Jour maintenant que j'ai regardé le réel fichiers de données. Il y a deux raisons que vous trouverez un décalage de phase de zéro. Tout d'abord, le décalage de phase est vraiment zéro entre les deux temps de la série. Vous pouvez voir clairement si vous effectuez un zoom avant à l'horizontale sur votre matplotlib graphique. Deuxièmement, il est important de régulariser les données de la première (le plus important, de soustraire au large de la moyenne), sinon l'effet de zéro-padding à la fin des tableaux marais le vrai signal de la corrélation croisée. Dans l'exemple suivant, je vérifie que je suis de trouver le "vrai" sommet par l'ajout d'une artificielle shift et ensuite vérifier que je le récupérer correctement.

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: Voici un exemple de la façon dont il fonctionne avec des données fausses.

EDIT 2: Ajouté un graphique de l'exemple.

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

De sorte qu'il fonctionne assez bien, même très bruyant données et de très aharmonic vagues.