User: reisswig Date: 2012/03/21 10:30 AM

Modified: /trunk/DataAnalysis/ DataAnalysis.py

Log: Add phase alignement procedure by Boyle et al.

File Changes:

Directory: /trunk/DataAnalysis/ ===============================

File [modified]: DataAnalysis.py Delta lines: +92 -0 =================================================================== --- trunk/DataAnalysis/DataAnalysis.py 2012-02-02 01:08:28 UTC (rev 7) +++ trunk/DataAnalysis/DataAnalysis.py 2012-03-21 15:30:39 UTC (rev 8) @@ -697,3 +697,95 @@ return err

+def PhaseAlignOverRange(DF1, DF2, x0, x1, start_xoff, delta, dx, accuracy=1e-2, monotonic=False): + """Phase-aligns two waveforms by minimizing their difference in the L2 norm and + returns the time-shift (and phase shift) that needs to be applied. The procedure is that + from Boyle et al, arxiv:0804.4184v2. + DF1: phase of waveform 1 + DF2: phase of waveform 2 + [x0,x1]: The minimization interval. + start_xoff: An initial guess to speed up computation. + delta: the shift will be applied in the range [start_xoff-delta, start_xoff+delta]. + The smaller delta, the faster (but keep in mind that you need to be sure that the minimum is located in that range! + dx: The delta spacing of interpolated function values on which the pointwise difference is taken. + accuracy: The accuracy in finding the shift that minimizes the L2 norm of the error ('the delta spacing of the shift'). + monotonic: Set to 'True' if function is monotonic (this will greatly speed up computation)""" + + # Note: this function is essentially identical to AlignOverRange in DiscreteFunction.py. + # But this also return a phase shift! + + # Keep DF1 fixed. Find offset x_off that needs to be applied to DF2.x + # in order to make the L2-norm of the difference over a range [x0,x1] + # as small as possible. + + x_off = start_xoff-delta + minErr = 10e10 + minX = x_off + trueMinX = x_off + trueDeltaPhi = 0 + + # If function is monotonic, we can use a faster algorithm + if (monotonic==True): + factor = 2.0*delta/accuracy*0.1 + maxX = start_xoff+delta + while (factor >= 1.0): + x_off = minX + minErr = 10e10 + while (x_off < maxX): + DF2shift = DiscreteFunction(DF2.x+x_off, DF2.f) + subDF1, subDF2 = MakeConsistent(DF1, DF2shift, x0, x1, dx) + # first, compute DeltaPhi analytically according to Boyle et al. + DeltaPhi = 0 + for ii in range (0, subDF1.Length()): + DeltaPhi += subDF1.f[ii]-subDF2.f[ii] + DeltaPhi /= (x1-x0) + # Now compute error + err = 0 + for ii in range (0, subDF1.Length()): + err += (subDF1.f[ii]-subDF2.f[ii]-DeltaPhi)**2 + err = sqrt(err) + #print minErr, err + if (err < minErr): + minErr = err + minX = x_off-accuracy*factor + trueMinX = x_off + trueDeltaPhi = DeltaPhi + else: + maxX = x_off + break + x_off += accuracy*factor + factor *= 0.1 + return trueMinX, trueDeltaPhi + + print "ERROR! ...using brute force method!" + + # Proceed with brute force algorithm... + + x_off = start_xoff-delta + + minErr = 10e10 + minX = x_off + minDeltaPhi = 0 + + + # We shift DF2 and compute the error over range [x0,x1] + while (x_off < start_xoff+delta): + DF2shift = DiscreteFunction(DF2.x+x_off, DF2.f) + subDF1, subDF2 = MakeConsistent(DF1, DF2shift, x0, x1, dx) + # first, compute DeltaPhi analytically according to Boyle et al. + DeltaPhi = 0 + for ii in range (0, subDF1.Length()): + DeltaPhi += subDF1.f[ii]-subDF2.f[ii] + DeltaPhi /= (x1-x0) + err = 0 + for ii in range (0, subDF1.Length()): + err += (subDF1.f[ii]-subDF2.f[ii]-DeltaPhi)**2 + err = sqrt(err) + #print err + if (err < minErr): + minErr = err + minX = x_off + minDeltaPhi = DeltaPhi + x_off += accuracy + + return minX, minDeltaPhi