Skip to content
Snippets Groups Projects
Commit 088e2692 authored by ChiWai Chan's avatar ChiWai Chan
Browse files

gstlal_inspiral_calc_snr: clean up the code, minor changes for command line options 

svd_bank_snr: In mode 1, use mkwhitened_multirate_src instead of mkwhiten
parent ad044820
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
gstlal-inspiral/bin/gstlal_inspiral_calc_snr
+ 20
− 39
View file @ 088e2692
......@@ -18,7 +18,7 @@ Typical Usages:
2. SVD bank options:
--svd-bank (require)
--bank-id (require)
--row-number
--row-number (optional, calculate all SNRs if not given)
3. PSD options:
--reference-psd (optional)
......@@ -28,7 +28,7 @@ Typical Usages:
4. Output options:
--output-width (default = 32bits)
--instrument (require)
--output (require)
--outdir (require)
--verbose (default = False)
--mode 1 (calculate SNR using Finite Impulse Response):
......@@ -43,10 +43,10 @@ Typical Usages:
2. Template options:
--table (require)
--approximant (require)
--template-duration (default = 600s)
--sample-rate (default = 16384Hz)
--template-duration (require)
--sample-rate (default = 2048Hz)
--f-low (default = 10)
--f-high (default = 1000)
--f-high (optional)
3. PSD / Whiten options:
--reference-psd (optional)
......@@ -57,26 +57,21 @@ Typical Usages:
--zero-pad (default = 0s)
4. Output options:
--output-rate (default = 2048Hz)
--output-width (default = 32bits)
--instrument (require)
--output (require)
--outdir (require)
--verbose (default = False)
"""
import numpy
from optparse import OptionParser, OptionGroup, IndentedHelpFormatter
import os
import sys
from gstlal import cbc_template_fir
from gstlal import chirptime
from gstlal import datasource
from gstlal import pipeparts
from gstlal import reference_psd
from gstlal import spawaveform
from gstlal import svd_bank
from gstlal import svd_bank_snr
from gstlal import templates
import lal
import lal.series
......@@ -125,8 +120,8 @@ def parse_command_line():
group.add_option("--row-number", type = "int", help = "The row number of the template (optional). All the SNRs will be outputed if it is not given.")
group.add_option("--table", metavar = "filename", help = "A LIGO light-weight xml.gz file containing SnglInspiral Table. Expecting one template for each instrument only.")
group.add_option("--approximant", metavar = "name", type = "str", help = "Name of the Waveform model (require).")
group.add_option("--template-duration", metavar = "seconds", default = 600, type = "int", help = "Duration of the template")
group.add_option("--sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sampling rate of the template")
group.add_option("--template-duration", metavar = "seconds", type = "int", help = "Duration of the template")
group.add_option("--sample-rate", metavar = "Hz", default = 2048, type = "int", help = "Sampling rate of the template and SNR for mode 1")
group.add_option("--f-low", metavar = "Hz", default = 10, type = "float", help = "The minimum frequency of GW signal")
group.add_option("--f-high", metavar = "Hz", type = "float", help = "The maximum frequency of GW signal")
parser.add_option_group(group)
......@@ -143,7 +138,6 @@ def parse_command_line():
group.add_option("--drop-first", metavar = "seconds", type = "int", default = 0, help = "Dropping the initital '--drop-first' seconds of SNR data (default = 0).")
group.add_option("--drop-last", metavar = "seconds", type = "int", default = 0, help = "Dropping the last '--drop-last' seconds of SNR data (default = 0).")
group.add_option("--output-width", metavar = "bits", type = "int", default = 32, help = "The size of the output data, can only be 32 or 64 bits (default = 32 bits).")
group.add_option("--output-rate", metavar = "Hz", type = "int", default = 2048, help = "The output sampling rate (default 2048).")
group.add_option("--instrument", metavar = "name", help = "The detector from which the --reference-psd and --frame-cache are loaded (require).")
parser.add_option_group(group)
......@@ -229,37 +223,22 @@ def parse_command_line():
missing_required_options.append("--outdir")
if options.instrument is None:
missing_required_options.append("--instrument")
if options.template_duration is None:
missing_required_options.append("--template-duration")
# Raise VauleError is required option(s) is/are missing
if missing_required_options:
raise ValueError("Missing required option(s) %s" % ", ".join(sorted(missing_required_options)))
#
# Temporary work for comparison
#
banks = svd_bank.read_banks(options.svd_bank, svd_bank.DefaultContentHandler)
bank = banks[options.sub_bank_id]
template_table = bank.sngl_inspiral_table
row = [template_table[options.row_number]]
time_slice = []
for frag in bank.bank_fragments[::-1]:
time_slice.append((frag.rate, frag.start, frag.end))
time_slice = numpy.array(time_slice, dtype = [("rate", "int"),("begin", "float"),("end", "float")])
workspace = cbc_template_fir.templates_workspace(template_table, options.approximant, psd[options.instrument], options.f_low, time_slice, autocorrelation_length = 351, fhigh= options.f_high)
#
# delele me afterward
#
#xmldoc = ligolw_utils.load_url(options.table, contenthandler = ContentHandler, verbose = options.verbose)
#template_table = table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName)
#row = [row for row in template_table if row.ifo == options.instrument]
#if len(row) != 1 :
# raise ValueError("Expecting only one template for --instrument=%s or cannot find template for --instrument=%s" %(options.instrument, options.instrument))
xmldoc = ligolw_utils.load_url(options.table, contenthandler = ContentHandler, verbose = options.verbose)
template_table = table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName)
row = [row for row in template_table if row.ifo == options.instrument]
if len(row) != 1 :
raise ValueError("Expecting only one template for --instrument=%s or cannot find template for --instrument=%s" %(options.instrument, options.instrument))
# work around for building a single whitened template
#time_slice = numpy.array([(options.sample_rate, 0, options.template_duration)], dtype = [("rate", "int"),("begin", "float"), ("end", "float")])
time_slice = numpy.array([(options.sample_rate, 0, options.template_duration)], dtype = [("rate", "int"),("begin", "float"), ("end", "float")])
# FIXME: psd[options.instrument] is bad
#workspace = cbc_template_fir.templates_workspace(template_table, options.approximant, psd[options.instrument], options.f_low, time_slice, options.f_high)
workspace = cbc_template_fir.templates_workspace(template_table, options.approximant, psd[options.instrument], options.f_low, time_slice, options.f_high)
template, autocorrelation, sigma = workspace.make_whitened_template(row[0])
return options, gw_data_source_info, template, psd
......@@ -284,6 +263,7 @@ if options.mode == 0:
ht_gate_threshold = options.ht_gate_threshold,
veto_segments = options.veto_segments,
track_psd = options.track_psd,
width = options.output_width,
verbose = options.verbose
)
snr_info = lloid_snr(gw_data_source_info, bank, options.instrument, psd = psd)
......@@ -302,11 +282,12 @@ elif options.mode == 1:
average_samples = options.average_samples,
median_samples = options.median_samples,
track_psd = options.track_psd,
rate = options.output_rate,
rate = options.sample_rate,
width = options.output_width,
verbose = options.verbose
)
#FIXME: allow multiple instruments
#FIXME: proper handle for latency
# drop the quadrature phase component
snr_info = fir_snr(gw_data_source_info, template.real, options.instrument, 0, psd = psd)
......
gstlal-inspiral/python/svd_bank_snr.py
+ 17
− 25
View file @ 088e2692
......@@ -4,9 +4,7 @@ Short cutting gstlal inspiral pipeline to produce SNR for template(s)
import sys
import numpy
import struct
import threading
from collections import deque
import gi
gi.require_version('Gst', '1.0')
......@@ -29,8 +27,8 @@ import lal.series
from ligo.lw import array as ligolw_array
from ligo.lw import ligolw
from ligo.lw import utils as ligolw_utils
from ligo.lw import param as ligolw_param
from ligo.lw import utils as ligolw_utils
@ligolw_array.use_in
@ligolw_param.use_in
......@@ -226,29 +224,23 @@ class FIR_SNR(object):
if self.verbose:
sys.stderr.write("Building pipeline to calculate SNR\n")
head, statevector, dqvector = datasource.mkbasicsrc(pipeline, gw_data_source_info, instrument, verbose = self.verbose)
src, statevector, dqvector = datasource.mkbasicsrc(pipeline, gw_data_source_info, instrument, verbose = self.verbose)
if psd is not None:
head = whiten = pipeparts.mkwhiten(pipeline, head, psd_mode = not self.track_psd, fft_length = self.psd_fft_length, zero_pad = self.zero_pad, average_samples = self.average_samples, median_samples = self.median_samples, name = "lal_whiten_%s" % instrument)
def psd_changes(elem, pspec, psd):
n = int(round(elem.get_property("f-nyquist")/elem.get_property("delta-f")) + 1)
psd = reference_psd.interpolate_psd(psd, elem.get_property("delta-f"))
elem.set_property("mean-psd", psd.data.data[:n])
whiten.connect_after("notify::f-nyquist", psd_changes, psd[instrument])
whiten.connect_after("notify::delta-f", psd_changes, psd[instrument])
whiten.connect_after("notify::psd-units", psd_changes, psd[instrument])
else:
head = pipeparts.mkwhiten(pipeline, head, psd_mode = 0, fft_length = self.psd_fft_length, zero_pad = self.zero_pad, average_samples = self.average_samples, median_samples = self.median_samples, name = "lal_whiten_%s" % instrument)
head = pipeparts.mkresample(pipeline, head)
if self.width == 32:
head = pipeparts.mkaudioconvert(pipeline, head, caps_string = "audio/x-raw, rate=%d, format=%s" %(self.rate, GstAudio.AudioFormat.to_string(GstAudio.AudioFormat.F32)))
elif self.width == 64:
head = pipeparts.mkaudioconvert(pipeline, head, caps_string = "audio/x-raw, rate=%d, format=%s" %(self.rate, GstAudio.AudioFormat.to_string(GstAudio.AudioFormat.F64)))
else:
raise ValueError("Invalid width: %d" % self.width)
#FIXME: how to set latency, block_stride, time_domain
head = pipeparts.mkfirbank(pipeline, head, latency = latency, fir_matrix = [template], block_stride = 16 * self.rate, time_domain = False)
hoftdict = multirate_datasource.mkwhitened_multirate_src(
pipeline,
src = src,
rates = [self.rate],
instrument = instrument,
psd = psd[instrument],
psd_fft_length = self.psd_fft_length,
track_psd = self.track_psd,
width = self.width,
statevector = statevector,
dqvector = dqvector
)
#FIXME: how to set latency
head = pipeparts.mkfirbank(pipeline, hoftdict[self.rate], latency = latency, fir_matrix = [template], block_stride = 16 * self.rate, time_domain = False)
appsink = pipeparts.mkappsink(pipeline, head, drop = False)
handler_id = appsink.connect("new-preroll", self.new_preroll_handler)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment