[Commits] [svn:einsteintoolkit] pyGWAnalysis/trunk/DataAnalysis/ (Rev. 8)

[reisswig at tapir.caltech.edu]

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