User: reisswig Date: 2011/06/14 06:19 PM
Added: /trunk/DataAnalysis/ WaveDiff.py
Modified: /trunk/DataAnalysis/ DataAnalysis.py /trunk/DiscreteFunction/ DiscreteFunction.py
Log: Introduced a new class for dealing with waveform difference computations more elegantly.
Fixed a bug in FFI regarding negative m-modes.
Optimizations.
Directory Changes:
Directory: /svn:executable/ ===========================
+ *
File Changes:
Directory: /trunk/DataAnalysis/ ===============================
File [modified]: DataAnalysis.py Delta lines: +3 -3 =================================================================== --- trunk/DataAnalysis/DataAnalysis.py 2011-04-18 22:27:52 UTC (rev 2) +++ trunk/DataAnalysis/DataAnalysis.py 2011-06-14 23:19:15 UTC (rev 3) @@ -266,9 +266,9 @@
for ii in range(0, fr.Length()): if (fr.x[ii] != 0): - if fr.x[ii] < 0: fr.f[ii] = 0 + #if fr.x[ii] < 0: fr.f[ii] = 0 div = 4*pi**2*fr.x[ii]**2 - if fr.x[ii] < frange[0]: div = 4*pi**2*frange[0]**2 + if abs(fr.x[ii]) < abs(frange[0]): div = 4*pi**2*frange[0]**2 #if fr.x[ii] > frange[1]: div = 4*pi**2*frange[1]**2 #if DF1fr.x[ii] > frange[0] and DF1fr.x[ii] < frange[1]: #div = (1.0-GaussBlend(DF1fr.x[ii], frange[0], (frange[1]-frange[0])/100)) * div + GaussBlend(DF1fr.x[ii], frange[0], (frange[1]-frange[0])/100) * 4*pi**2*frange[0]**2 # div = (1 - (DF1fr.x[ii]-frange[0])/(frange[1]-frange[0])) * 4*pi**2*frange[0]**2 + (DF1fr.x[ii]-frange[0])/(frange[1]-frange[0]) * div @@ -300,7 +300,7 @@ if (fr.x[ii] != 0): if fr.x[ii] < 0: fr.f[ii] = 0 div = 2*pi*fr.x[ii] - if fr.x[ii] < frange[0]: div = 2*pi*frange[0] + if abs(fr.x[ii]) < abs(frange[0]): div = 2*pi*frange[0]*sign(fr.x[ii]) #if fr.x[ii] > frange[1]: div = 4*pi**2*frange[1]**2 #if DF1fr.x[ii] > frange[0] and DF1fr.x[ii] < frange[1]: #div = (1.0-GaussBlend(DF1fr.x[ii], frange[0], (frange[1]-frange[0])/100)) * div + GaussBlend(DF1fr.x[ii], frange[0], (frange[1]-frange[0])/100) * 4*pi**2*frange[0]**2 # div = (1 - (DF1fr.x[ii]-frange[0])/(frange[1]-frange[0])) * 4*pi**2*frange[0]**2 + (DF1fr.x[ii]-frange[0])/(frange[1]-frange[0]) * div
File [added]: WaveDiff.py Delta lines: +215 -0 =================================================================== --- trunk/DataAnalysis/WaveDiff.py (rev 0) +++ trunk/DataAnalysis/WaveDiff.py 2011-06-14 23:19:15 UTC (rev 3) @@ -0,0 +1,215 @@ +#!/usr/bin/python + +#---------------------------------------------------------- +# +# This module computes differences in two given waveforms. +# +# +# This module depends on DataAnalysis and DiscreteFunction +# +# Released under the MIT License. +# (C) Christian Reisswig 2010-2011 +# +#---------------------------------------------------------- + + +import sys +import os +import glob +from DiscreteFunction import * +from WaveFunction import * +from harmonics import * +from DataAnalysis import * + + + +class WaveDiff: + """A class holding two wavefunction objects + and their differences""" + + def __init__(self, WF1, WF2, dx, x0, x1): + """Initializes this class by loading two wavefunction. + The wavefunctions are assumed to be in complex form (_not_ decomposed into amplitude and phase!). + The waveforms will be resampled to new dx at common points and restricted to interval defined by [x0,x1].""" + self.WF1 = WaveFunction(WF1.x, WF1.f) + self.WF2 = WaveFunction(WF2.x, WF2.f) + self.align = "" + self.phi = 0 + self.omega = 0 + self.dx = dx + self.x0 = x0 + self.x1 = x1 + self.isAligned = False + # we cache the aligned waveforms + self.aWF1 = WaveFunction([],[]) + self.aWF2 = WaveFunction([],[]) + + + def AlignInPhase(self, phi): + """Align given waveforms at same phase value phi""" + self.align = "phase" + self.phi = phi + self.isAligned = False + + def AlignInFreq(self, omega): + """Align given waveforms at same frequency omega""" + self.align = "freq" + self.omega = omega + self.isAligned = False + + def AlignAtAmpMax(self): + """Align two waveforms at their amplitude maximum""" + self.align = "max" + self.isAligned = False + + def NoAlign(self): + self.align = "" + self.isAligned = False + + def Align(self): + """Returns aligned waveforms according to settings""" + + if self.isAligned: + # return copy of cached waveforms + return WaveFunction(self.aWF1.x, self.aWF1.f), WaveFunction(self.aWF2.x, self.aWF2.f) + + if (self.align == ""): + aWF1 = WaveFunction(self.WF1.x, self.WF1.f) + aWF2 = WaveFunction(self.WF2.x, self.WF2.f) + if (self.align == "max"): + # align at amplitude maximum + tmax1 = self.WF1.Amplitude().FindAbsMaxInterpolated(1e-6) + tmax2 = self.WF2.Amplitude().FindAbsMaxInterpolated(1e-6) + aWF1 = WaveFunction(self.WF1.x-tmax1, self.WF1.f) + aWF2 = WaveFunction(self.WF2.x-tmax2, self.WF2.f) + if (self.align == "freq"): + # get omega(t) + freq1 = self.WF1.Frequency() + freq2 = self.WF2.Frequency() + # invert omega(t) -> t(omega) + t1_vs_omega1 = DiscreteFunction(abs(freq1.f), freq1.x) + t2_vs_omega2 = DiscreteFunction(abs(freq2.f), freq2.x) + t1_vs_omega1 = t1_vs_omega1.MakeMonotonicCoordinates() + t2_vs_omega2 = t2_vs_omega2.MakeMonotonicCoordinates() + # find time t at which WF has frequency omega + t1 = t1_vs_omega1.Interpolate(self.omega) + t2 = t2_vs_omega2.Interpolate(self.omega) + #print t1, t2 + # align waveforms at that frequency + aWF1 = WaveFunction(self.WF1.x-t1, self.WF1.f) + aWF2 = WaveFunction(self.WF2.x-t2, self.WF2.f) + + #self.aWF1 = WaveFunction(freq1.x-t1, freq1.f) + #self.aWF2 = WaveFunction(freq2.x-t2, freq2.f) + + # interpolate waveforms to same points in same intervall + self.aWF1, self.aWF2 = MakeConsistent(aWF1, aWF2, self.x0, self.x1, self.dx) + + self.isAligned = True + + return WaveFunction(self.aWF1.x, self.aWF1.f), WaveFunction(self.aWF2.x, self.aWF2.f) + + + def AlignedPhases(self, reparam=""): + """Compute aligned/reparametrized phase as function of time, phase or frequency""" + + aWF1, aWF2 = self.Align() + + phase1 = aWF1.Phase() + phase2 = aWF2.Phase() + + if (self.align == "max"): + # shift phase to zero at t=0 (where the maximum is) + shift1 = phase1.Interpolate(0.0) + shift2 = phase2.Interpolate(0.0) + + phase1.f -= shift1 + phase2.f -= shift2 + + return phase1, phase2 + + + def AlignedAmps(self, reparam=""): + """Computes aligned/reparametrized amplitude as function of time, phase or frequency""" + + aWF1, aWF2 = self.Align() + + amp1 = aWF1.Amplitude() + amp2 = aWF2.Amplitude() + + if (reparam == "phase"): + # Reparametrize in terms of phase... + # ...and align, if necessary + phase1 = aWF1.Phase() + phase2 = aWF2.Phase() + if (self.align == "max"): + # shift phase to zero at t=0 (where the maximum is) + shift1 = phase1.Interpolate(0.0) + shift2 = phase2.Interpolate(0.0) + + phase1.f -= shift1 + phase2.f -= shift2 + + # we take the absolute value of the phase since the phase can be negative and then our (positive) monotonicity does not hold + amp1_vs_phase1 = DiscreteFunction(abs(phase1.f), amp1.f) + amp1_vs_phase1 = amp1_vs_phase1.MakeMonotonicCoordinates() + + amp2_vs_phase2 = DiscreteFunction(abs(phase2.f), amp2.f) + amp2_vs_phase2 = amp2_vs_phase2.MakeMonotonicCoordinates() + + amp1_vs_phase1, amp2_vs_phase2 = MakeConsistent(amp1_vs_phase1, amp2_vs_phase2, self.x0, self.x1, self.dx) + + return amp1_vs_phase1, amp2_vs_phase2 + + # if no reparametrization is requested, just return difference in time + return amp1, amp2 + + + def DiffAmp(self, reparam=""): + """Computes difference in amplitude as function of time, phase or frequency""" + + a, b = self.AlignedAmps(reparam) + + return a-b + + def DiffPhase(self, reparam=""): + """Computes difference in phase as function of time, phase or frequency""" + + a, b = self.AlignedPhases(reparam) + + return a-b + + + def Diff(self, reparam=""): + """Computes the complex wave difference as function of time, phase or frequency""" + + aWF1, aWF2 = self.Align() + + return diffWF + + + def Mismatch(self, mass, ftype="h", Sn=SensitivityAdLIGO, frange=[10,8000]): + """Computes the mismatch for a given mass (in Msol), assuming the wave-function to be any of ftype=(h,news,psi4). + Sn is the detector sensitivity function, frange the (detector) frequency range in Hz.""" + + aWF1, aWF2 = self.Align() + + MM = -1 + + if ftype == "h": + aWF1SI = StrainToSI(aWF1, mass) + aWF2SI = StrainToSI(aWF2, mass) + MM = MisMatch(aWF1SI, aWF2SI, Sn, frange[0], frange[1]) + if ftype == "news": + aWF1SI = NewsToSI(aWF1, mass) + aWF2SI = NewsToSI(aWF2, mass) + MM = MisMatch(aWF1SI, aWF2SI, Sn, frange[0], frange[1], 'news') + if ftype == "psi4": + aWF1SI = Psi4ToSI(aWF1, mass) + aWF2SI = Psi4ToSI(aWF2, mass) + MM = MisMatch(aWF1SI, aWF2SI, Sn, frange[0], frange[1], 'psi4') + + return MM + + + \ No newline at end of file
Property changes on: trunk/DataAnalysis/WaveDiff.py ___________________________________________________________________
Directory: /trunk/DiscreteFunction/ ===================================
File [modified]: DiscreteFunction.py Delta lines: +12 -5 =================================================================== --- trunk/DiscreteFunction/DiscreteFunction.py 2011-04-18 22:27:52 UTC (rev 2) +++ trunk/DiscreteFunction/DiscreteFunction.py 2011-06-14 23:19:15 UTC (rev 3) @@ -70,15 +70,18 @@ else: self.dx = 0.0
- def MakeMonotonicCoordinates(self): - """Given a function whose coordinates are not monotonic, we remove all entries that violate monotonicity""" + def MakeMonotonicCoordinates(self, maxvariance=1e10): + """Given a function whose coordinates are not monotonic, we remove all entries that violate monotonicity. + maxvariance is the maximal change that may occur from one point to the next""" x_ = [] f_ = [] - + jj = 0 for ii in range(0, len(self.x)-1): if (self.x[jj] >= self.x[ii]): continue + if (self.x[ii] > self.x[jj]+maxvariance): + continue jj = ii x_.append(self.x[ii]) f_.append(self.f[ii]) @@ -509,7 +512,7 @@ res = self.__class__(zeros(npoints, float), zeros(npoints, type(self.f[0])))
dt = self.dx - res.f = dt * fft.fft(self.f) #/2 # divide by two becuase we're going to cut the negative frequencies! + res.f = dt * fft.fft(self.f) #/2 # divide by two because we're going to cut the negative frequencies!
# get the freq
@@ -638,7 +641,11 @@ if (x0 < max(DF1.x[0], DF2.x[0])): x0 = max(DF1.x[0], DF2.x[0]) if (x1 > min(DF1.x[-1], DF2.x[-1])): x1 = min(DF1.x[-1], DF2.x[-1])
- npoints = int ((x1-x0)/dx)+1 + if ((x1-x0) < 0): + print "Warning: MakeConsistent: Empty overlap." + npoints = 1 + else: + npoints = int ((x1-x0)/dx)+1
newDF1 = DF2.__class__(zeros(npoints, float), zeros(npoints, type(DF1.f[0]))) newDF2 = DF2.__class__(zeros(npoints, float), zeros(npoints, type(DF2.f[0])))