StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>Perhaps you are looking for the cross-correlation:</p> <pre><code>scipy.signal.signaltools.correlate(A, B) </code></pre> <p>The position of the peak in the cross-correlation will be an estimate of the phase difference. </p> <p><strong>EDIT 3:</strong> 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.</p> <pre><code>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 </code></pre> <p><strong>EDIT:</strong> Here is an example of how it works with fake data.</p> <p><strong>EDIT 2:</strong> Added a graph of the example.</p> <p><img src="https://i.stack.imgur.com/V0a7V.png" alt="Cross correlation of noisy anharmonic signals"></p> <pre><code>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") </code></pre> <p>So it works pretty well, even for very noisy data and very aharmonic waves. </p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload