Skip to content
Snippets Groups Projects
Commit 39f9e3ee authored by Patrick Godwin's avatar Patrick Godwin
Browse files

gstlal_etg_template_overlap: stripped out waveform functions copied over in...

gstlal_etg_template_overlap: stripped out waveform functions copied over in gstlal_etg in favor of refactored waveform class in idq_utils.py
parent 6f7c0578
No related branches found
No related tags found
No related merge requests found
......@@ -33,7 +33,7 @@ import numpy
import lal
#from gstlal import idq_waveforms
from gstlal import idq_utils
###############################
#
......@@ -54,87 +54,6 @@ def parse_command_line():
return options, args
####################
#
# functions
#
####################
### FIXME: should not copy functions over from gstlal_etg
### instead, factor out functions from gstlal_etg to its own module
#
# construct sine gaussian waveforms that taper to 1e-7 at edges of window
#
def duration(phi, q, tolerance = 1e-7):
# return the duration of the waveform such that its edges will die out to tolerance of the peak.
# phi is the central frequency of the frequency band
return q/(2.*numpy.pi*phi)*numpy.log(1./tolerance)
def sine_gaussian(phi, phi_0, q, time_arr):
# edges should be 1e-7 times peak
dt = time_arr[1]-time_arr[0]
rate = 1./dt
assert phi < rate/2.
# phi is the central frequency of the sine gaussian
dur = duration(phi,q)
tau = q/(2.*numpy.pi*phi)
sg_vals = numpy.cos(2.*numpy.pi*phi*time_arr + phi_0)*numpy.exp(-1.*time_arr**2./tau**2.)
# normalize sine gaussians to have unit length in their vector space
#inner_product = numpy.sum(sg_vals*sg_vals)
#norm_factor = 1./(inner_product**0.5)
return sg_vals
def half_sine_gaussian(sg_arr):
samples = sg_arr.size/2 + 1
# only take first half of sine gaussian + peak
hsg_vals = sg_arr[:samples]
# renormalize
inner_product = numpy.sum(hsg_vals*hsg_vals)
norm_factor = 1./(inner_product**0.5)
return norm_factor*hsg_vals
#
# number of tiles in frequency and Q
#
def N_Q(q_min, q_max, mismatch = 0.2):
"""
Minimum number of distinct Q values to generate based on Q_min, Q_max, and mismatch params.
"""
return numpy.ceil(1./(2.*numpy.sqrt(mismatch/3.))*1./numpy.sqrt(2)*numpy.log(q_max/q_min))
def N_phi(phi_min, phi_max, q, mismatch = 0.2):
"""
Minimum number of distinct frequency (phi) values to generate based on phi_min, phi_max, and mismatch params.
"""
return numpy.ceil(1./(2.*numpy.sqrt(mismatch/3.))*(numpy.sqrt(2.+q**2.)/2.)*numpy.log(phi_max/phi_min))
def Qq(q_min, q_max, mismatch = 0.2):
"""
List of Q values to generate based on Q_min, Q_max, and mismatch params.
"""
N_q = numpy.ceil(N_Q(q_min, q_max, mismatch = mismatch))
return [q_min*(q_max/q_min)**((0.5+q)/N_q) for q in range(int(N_q))]
def phi_ql(phi_min, phi_max, q_min, q_max, mismatch = 0.2, phi_nyquist = None):
"""
Generates Q and frequency (phi) pairs based on mismatch, Q and phi ranges.
"""
for q in Qq(q_min, q_max, mismatch = mismatch):
nphi = N_phi(phi_min, phi_max, q, mismatch = mismatch)
for l in range(int(nphi)):
phi = phi_min*(phi_max/phi_min)**((0.5+l)/nphi)
if phi_nyquist:
if phi < phi_nyquist/(1+numpy.sqrt(11)/q):
yield (phi, q)
else:
yield (phi, q)
####################
#
# main
......@@ -157,10 +76,9 @@ if __name__ == '__main__':
fhigh = rate/2.*0.8
qhigh = options.qhigh
qlow = 4
dur = max([duration(phi, q, 5e-3) for (phi, q) in phi_ql(flow, fhigh, qlow, qhigh, mismatch = options.mismatch)])
t_arr = numpy.linspace(-dur/2., dur/2., int(numpy.floor(dur*rate)+1))
#basis_params[rate] = [(phi, q, duration(phi, q, 5e-3)/2.) for (phi, q) in phi_ql(flow, fhigh, qlow, qhigh, mismatch = options.mismatch, phi_nyquist = fhigh)]
templates[rate] = [half_sine_gaussian(sine_gaussian(phi, 0, q, t_arr)) for (phi, q) in phi_ql(flow, fhigh, qlow, qhigh, mismatch = options.mismatch, phi_nyquist = fhigh)]
hsg_waveforms = idq_utils.HalfSineGaussianGenerator((flow, fhigh), (qlow, qhigh), rate, mismatch = options.mismatch)
#basis_params[(channel, rate)] = hsg_waveforms.parameter_grid
templates[rate] = [waveform for waveform in hsg_waveforms.generate_templates()]
# calculate overlap for each template, find maximum
overlaps = {rate: [] for rate in rates}
......
......@@ -272,4 +272,4 @@ class HalfSineGaussianGenerator(object):
for l in range(num_f):
f = f_min * (f_max/f_min)**( (0.5+l) /num_f)
if f < f_max / (1 + (numpy.sqrt(11)/q)):
yield (f, q)
yield (f, q)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment