diff --git a/gstlal-ugly/bin/gstlal_feature_extractor_template_overlap b/gstlal-ugly/bin/gstlal_feature_extractor_template_overlap
index 4b0a2534cb12e38f53e7c2510a98fe940cbfa620..2ae23886b1b62ed6aa7b57498537ccdc5ca0d774 100755
--- a/gstlal-ugly/bin/gstlal_feature_extractor_template_overlap
+++ b/gstlal-ugly/bin/gstlal_feature_extractor_template_overlap
@@ -233,7 +233,6 @@ if __name__ == '__main__':
aggregator.makedir(options.output_dir)
# common parameters we will use throughout
- #max_samp_rate = 2048
max_samp_rate = 8192
min_samp_rate = 32
n_rates = int(numpy.log2(max_samp_rate/min_samp_rate) + 1)
@@ -243,22 +242,20 @@ if __name__ == '__main__':
# generate templates for each rate considered
rates = [min_samp_rate*2**i for i in range(n_rates)]
- waveforms = {}
- templates = {}
- basis_params = {}
- for rate in rates:
- flow = rate/4.*0.8
- fhigh = rate/2.*0.8
- qhigh = options.qhigh
- qlow = 3.3166
- if options.waveform == 'sine_gaussian':
- waveforms[rate] = idq_utils.SineGaussianGenerator((flow, fhigh), (qlow, qhigh), rate, mismatch = options.mismatch)
- elif options.waveform == 'half_sine_gaussian':
- waveforms[rate] = idq_utils.HalfSineGaussianGenerator((flow, fhigh), (qlow, qhigh), rate, mismatch = options.mismatch)
- else:
- raise NotImplementedError
- basis_params[rate] = waveforms[rate].parameter_grid
- templates[rate] = [waveform for waveform in waveforms[rate].generate_templates(quadrature=False)]
+ downsample_factor = 0.8
+ qhigh = options.qhigh
+ qlow = 3.3166
+ fhigh = max_samp_rate / 2.
+ flow = min_samp_rate / 4.
+
+ if options.waveform == 'sine_gaussian':
+ waveforms = idq_utils.SineGaussianGenerator((flow, fhigh), (qlow, qhigh), rates, mismatch = options.mismatch, downsample_factor=downsample_factor)
+ elif options.waveform == 'half_sine_gaussian':
+ waveforms = idq_utils.HalfSineGaussianGenerator((flow, fhigh), (qlow, qhigh), rates, mismatch = options.mismatch, downsample_factor=downsample_factor)
+ else:
+ raise NotImplementedError
+ basis_params = waveforms.parameter_grid
+ templates = {rate: [waveform for waveform in waveforms.generate_templates(rate, quadrature=False)] for rate in rates}
# get all templates and params
all_templates = list(itertools.chain(*templates.values()))
@@ -267,11 +264,16 @@ if __name__ == '__main__':
if options.verbose:
print >>sys.stderr, "Creating template overlaps..."
+ # zero pad templates to make them the same length
+ max_sample_pts = max(len(template) for template in all_templates)
+ all_templates = [numpy.pad(template, ((max_sample_pts - len(template)) // 2, (max_sample_pts - len(template)) // 2), 'constant') for template in all_templates]
+
# calculate overlap for each template and find maximum
overlaps = []
for this_template in all_templates:
overlaps.append(max([numpy.dot(this_template,template) for template in all_templates if not numpy.array_equal(template, this_template)]))
+ print >>sys.stderr, "total number of templates: %d" % len(all_templates)
print >>sys.stderr, "min overlap specified: %f" % (1 - options.mismatch)
print >>sys.stderr, "max template overlap: %f" % max(overlaps)
print >>sys.stderr, "min template overlap: %f" % min(overlaps)
@@ -291,13 +293,15 @@ if __name__ == '__main__':
print >>sys.stderr, "Creating waveform plots..."
for rate in rates:
- for template_id, template in enumerate(random.sample(templates[rate], num_samples)):
+ #for template_id, template in enumerate(random.sample(templates[rate], num_samples)):
+ for template_id in random.sample(numpy.arange(len(templates[rate])), num_samples):
waveform_params = ["%s: %.3f" % (name, param) for param, name in zip(basis_params[rate][template_id], param_names)]
+ template = templates[rate][template_id]
if options.verbose:
print >>sys.stderr, "\tCreating waveform plot with parameters: %s" % repr(waveform_params)
- series_fig = plot_waveform(waveforms[rate].times, template, waveform_type, waveform_params)
+ series_fig = plot_waveform(waveforms.times[rate], template, waveform_type, waveform_params)
fname = 'plot_%s_%s-timeseries.png' % (str(rate).zfill(4), str(template_id).zfill(4))
plot_paths.append((template_id*int(rate),fname))
series_fig.savefig(os.path.join(options.output_dir, fname))
diff --git a/gstlal-ugly/python/idq_utils.py b/gstlal-ugly/python/idq_utils.py
index 424b934be49aa6780bdf630f0af0fa840de4135f..a8257fa1d079b61bc57a34b2379d2d43e1f85511 100644
--- a/gstlal-ugly/python/idq_utils.py
+++ b/gstlal-ugly/python/idq_utils.py
@@ -25,6 +25,8 @@
####################
+import bisect
+from collections import defaultdict
import glob
import logging
import os
@@ -296,32 +298,59 @@ class HalfSineGaussianGenerator(object):
"""
Generates half sine gaussian templates based on a f, Q range and a sampling frequency.
"""
- def __init__(self, f_range, q_range, f_samp, mismatch=0.2, tolerance=5e-3):
+ def __init__(self, f_range, q_range, rates, mismatch=0.2, tolerance=5e-3, downsample_factor=0.8):
+
### define parameter range
- self.f_low, self.f_high = f_range
+ self.f_low, self.f_high = (downsample_factor * f_range[0], downsample_factor * f_range[1])
self.q_low, self.q_high = q_range
- self.f_samp = f_samp
+
+ ### define grid spacing and edge rolloff for templates
self.mismatch = mismatch
self.tol = tolerance
+ ### determine how to scale down templates considered
+ ### in frequency band to deal with low pass rolloff
+ self.downsample_factor = downsample_factor
+
+ ### determine frequency bands considered
+ self.rates = sorted(set(rates))
+
+ ### determine lower frequency limits for rate conversion based on downsampling factor
+ self.freq_breakpoints = [self.downsample_factor * rate / 2. for rate in self.rates[:-1]]
+
### define grid and template duration
- self.parameter_grid = [(f, q, self.duration(f, q)) for f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high)]
- self.max_duration = max(duration for f, q, duration in self.parameter_grid)
+ self.parameter_grid = defaultdict(list)
+ for rate, f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high):
+ self.parameter_grid[rate].append((f, q, self.duration(f, q)))
+ self.max_duration = {rate: max(duration for f, q, duration in parameters) for rate, parameters in self.parameter_grid.items()}
+
+ ### determine rest of waveform properties
self.phases = [0., numpy.pi/2.]
- self.times = numpy.linspace(-self.max_duration, 0, int(numpy.ceil(self.max_duration*self.f_samp)), endpoint=True)
- self.latency = 0
+ self.times = {rate: numpy.linspace(-self.max_duration[rate], 0, self.round_to_next_odd(self.max_duration[rate]*rate), endpoint=True) for rate in self.rates}
+ self.latency = {rate: 0 for rate in self.rates}
+
+ def frequency_to_rate(self, frequency):
+ """
+ Maps a frequency to the correct sampling rate considered.
+ """
+ idx = bisect.bisect(self.freq_breakpoints, frequency)
+ return self.rates[idx]
+
+ def round_to_next_odd(self, n):
+ return int(numpy.ceil(n) // 2 * 2 + 1)
- def generate_templates(self, quadrature = True):
+ def generate_templates(self, rate, quadrature = True):
"""
generate all half sine gaussian templates corresponding to a parameter range and template duration
+ for a given sampling rate.
"""
- for f, q in self.generate_f_q_grid(self.f_low, self.f_high, self.q_low, self.q_high):
+ for f, q, _ in self.parameter_grid[rate]:
if quadrature:
for phase in self.phases:
yield self.waveform(f, q, phase)
else:
- yield self.waveform(f, q, self.phases[0])
+ yield self.waveform(f, q, self.phases[0], rate)
def duration(self, f, q):
"""
@@ -329,18 +358,16 @@ class HalfSineGaussianGenerator(object):
"""
return 0.5 * (q/(2.*numpy.pi*f)) * numpy.log(1./self.tol)
- def waveform(self, f, q, phase):
+ def waveform(self, f, q, phase, rate):
"""
construct half sine gaussian waveforms that taper to tolerance at edges of window
f is the central frequency of the waveform
"""
- dt = self.times[1] - self.times[0]
- rate = 1./dt
assert f < rate/2.
# phi is the central frequency of the sine gaussian
tau = q/(2.*numpy.pi*f)
- template = numpy.cos(2.*numpy.pi*f*self.times + phase)*numpy.exp(-1.*self.times**2./tau**2.)
+ template = numpy.cos(2.*numpy.pi*f*self.times[rate] + phase)*numpy.exp(-1.*self.times[rate]**2./tau**2.)
# normalize sine gaussians to have unit length in their vector space
inner_product = numpy.sum(template*template)
@@ -369,23 +396,26 @@ class HalfSineGaussianGenerator(object):
def generate_f_q_grid(self, f_min, f_max, q_min, q_max):
"""
- Generates (f, Q) pairs based on f, Q ranges
+ Generates (f_samp, f, Q) pairs based on f, Q ranges
"""
for q in self.generate_q_values(q_min, q_max):
num_f = self.num_f_templates(f_min, f_max, q)
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)
+ rate = self.frequency_to_rate(f)
+ if f < (rate/2) / (1 + (numpy.sqrt(11)/q)):
+ yield rate, f, q
+ elif rate != max(self.rates):
+ yield (2 * rate), f, q
class SineGaussianGenerator(HalfSineGaussianGenerator):
"""
Generates sine gaussian templates based on a f, Q range and a sampling frequency.
"""
- def __init__(self, f_range, q_range, f_samp, mismatch=0.2, tolerance=5e-3):
- super(SineGaussianGenerator, self).__init__(f_range, q_range, f_samp, mismatch=mismatch, tolerance=tolerance)
- self.times = numpy.linspace(-self.max_duration/2., self.max_duration/2., int(numpy.ceil(self.max_duration*self.f_samp)), endpoint=True)
- self.latency = self.max_duration / 2.
+ def __init__(self, f_range, q_range, rates, mismatch=0.2, tolerance=5e-3, downsample_factor=0.8):
+ super(SineGaussianGenerator, self).__init__(f_range, q_range, rates, mismatch=mismatch, tolerance=tolerance, downsample_factor=0.8)
+ self.times = {rate: numpy.linspace(-self.max_duration[rate]/2., self.max_duration[rate]/2., self.round_to_next_odd(self.max_duration[rate]*rate), endpoint=True) for rate in self.rates}
+ self.latency = {rate: self.max_duration[rate] / 2. for rate in self.rates}
def duration(self, f, q):
"""