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

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

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

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