feature_extraction: change block_stride for whitening, high pass from 1s to...

feature_extraction: change block_stride for whitening, high pass from 1s to 0.25s, allow for a minimum downsampling rate to improve latency at expensive of computational cost

gstlal-burst/bin/gstlal_feature_extractor

@@ -221,6 +221,7 @@ def parse_command_line():
 
 	# check if input sample rate is sensible
 	assert options.sample_rate == 1 or options.sample_rate % 2 == 0
+	assert options.min_downsample_rate % 2 == 0
 
 	# check if persist and save cadence times are sensible
 	assert options.persist_cadence >= options.cadence
@@ -418,24 +419,25 @@ for subset_id, channel_subset in enumerate(data_source_info.channel_subsets, 1):
 			head[channel] = pipeparts.mklatency(pipeline, head[channel], name=utils.latency_name('beforewhitening', 2, channel))
 
 		# whiten auxiliary channel data
-		for rate, thishead in multirate_datasource.mkwhitened_multirate_src(pipeline, head[channel], rates, samp_rate, instrument, channel_name = channel, width=32, nxydump_segment=options.nxydump_segment, psd_fft_length=options.psd_fft_length).items():
+		for rate, thishead in multirate_datasource.mkwhitened_multirate_src(pipeline, head[channel], rates, samp_rate, instrument, channel_name = channel, width=32, nxydump_segment=options.nxydump_segment, psd_fft_length=options.psd_fft_length, min_rate=options.min_downsample_rate).items():
+			thisrate = max(options.min_downsample_rate, rate)
 			if options.latency_output:
 				thishead = pipeparts.mklatency(pipeline, thishead, name=utils.latency_name('afterwhitening', 3, channel, rate))
 
 			# determine whether to do time-domain or frequency-domain convolution
-			time_domain = (waveforms[channel].sample_pts(rate)*rate) < (5*waveforms[channel].sample_pts(rate)*numpy.log2(rate))
+			time_domain = (waveforms[channel].sample_pts(thisrate)*thisrate) < (5*waveforms[channel].sample_pts(thisrate)*numpy.log2(thisrate))
 
 			# create FIR bank of half sine-gaussian templates
-			fir_matrix = numpy.array([waveform for waveform in waveforms[channel].generate_templates(rate)])
+			fir_matrix = numpy.array([waveform for waveform in waveforms[channel].generate_templates(rate, sampling_rate=thisrate)])
 			thishead = pipeparts.mkqueue(pipeline, thishead, max_size_buffers = 0, max_size_bytes = 0, max_size_time = Gst.SECOND * 30)
-			thishead = pipeparts.mkfirbank(pipeline, thishead, fir_matrix = fir_matrix, time_domain = time_domain, block_stride = int(rate), latency = waveforms[channel].latency(rate))
+			thishead = pipeparts.mkfirbank(pipeline, thishead, fir_matrix = fir_matrix, time_domain = time_domain, block_stride = int(thisrate), latency = waveforms[channel].latency(thisrate))
 
 			# add queues, change stream format, add tags
 			if options.latency_output:
 				thishead = pipeparts.mklatency(pipeline, thishead, name=utils.latency_name('afterFIRbank', 4, channel, rate))
 			thishead = pipeparts.mkqueue(pipeline, thishead, max_size_buffers = 1, max_size_bytes = 0, max_size_time = 0)
 			thishead = pipeparts.mktogglecomplex(pipeline, thishead)
-			thishead = pipeparts.mkcapsfilter(pipeline, thishead, caps = "audio/x-raw, format=Z64LE, rate=%i" % rate)
+			thishead = pipeparts.mkcapsfilter(pipeline, thishead, caps = "audio/x-raw, format=Z64LE, rate=%i" % thisrate)
 			thishead = pipeparts.mktaginject(pipeline, thishead, "instrument=%s,channel-name=%s" %( instrument, channel))
 
 			# dump segments to disk if specified
@@ -445,9 +447,9 @@ for subset_id, channel_subset in enumerate(data_source_info.channel_subsets, 1):
 
 			# extract features from time series
 			if options.feature_mode == 'timeseries':
-				thishead = pipeparts.mktrigger(pipeline, tee, int(rate // options.sample_rate), max_snr = True)
+				thishead = pipeparts.mktrigger(pipeline, tee, int(thisrate // options.sample_rate), max_snr = True)
 			elif options.feature_mode == 'etg':
-				thishead = pipeparts.mktrigger(pipeline, tee, rate, snr_thresh = options.snr_threshold)
+				thishead = pipeparts.mktrigger(pipeline, tee, thisrate, snr_thresh = options.snr_threshold)
 
 			if options.latency_output:
 				thishead = pipeparts.mklatency(pipeline, thishead, name=utils.latency_name('aftertrigger', 5, channel, rate))

gstlal-burst/python/fxtools/feature_extractor.py

@@ -534,6 +534,7 @@ def append_options(parser):
 
 	group = optparse.OptionGroup(parser, "Program Behavior")
 	group.add_option("--psd-fft-length", metavar = "seconds", default = 32, type = "int", help = "The length of the FFT used to used to whiten the data (default is 32 s).")
+	group.add_option("--min-downsample-rate", metavar = "Hz", default = 128, type = "int", help = "The minimum sampling rate in which to downsample streams. Default = 128 Hz.")
 	group.add_option("--local-frame-caching", action = "store_true", help = "Pre-reads frame data and stores to local filespace.")
 	group.add_option("--disable-web-service", action = "store_true", help = "If set, disables web service that allows monitoring of PSDS of aux channels.")
 	group.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")

gstlal-burst/python/fxtools/multirate_datasource.py

@@ -57,7 +57,7 @@ NATIVE_RATE_CUTOFF = 128
 # =============================================================================
 #
 
-def mkwhitened_multirate_src(pipeline, src, rates, native_rate, instrument, psd = None, psd_fft_length = PSD_FFT_LENGTH, veto_segments = None, nxydump_segment = None, track_psd = True, block_duration = int(1 * Gst.SECOND), width = 64, channel_name = "hoft"):
+def mkwhitened_multirate_src(pipeline, src, rates, native_rate, instrument, psd = None, psd_fft_length = PSD_FFT_LENGTH, veto_segments = None, nxydump_segment = None, track_psd = True, block_duration = 0.25 * Gst.SECOND, width = 64, channel_name = "hoft", min_rate = 128):
 	"""!
 	Build pipeline stage to whiten and downsample auxiliary channels.
 
@@ -157,7 +157,7 @@ def mkwhitened_multirate_src(pipeline, src, rates, native_rate, instrument, psd
 	# high pass filter
 	#
 
-	block_stride = block_duration * max_rate // Gst.SECOND
+	block_stride = int(block_duration * max_rate // Gst.SECOND)
 	if native_rate >= NATIVE_RATE_CUTOFF:
 		kernel = reference_psd.one_second_highpass_kernel(max_rate, cutoff = 12)
 		assert len(kernel) % 2 == 1, "high-pass filter length is not odd"
@@ -266,7 +266,7 @@ def mkwhitened_multirate_src(pipeline, src, rates, native_rate, instrument, psd
 	for rate in sorted(set(rates))[:-1]:
 		head[rate] = pipeparts.mkqueue(pipeline, tee, max_size_buffers = 0, max_size_bytes = 0, max_size_time = Gst.SECOND * 8)
 		head[rate] = pipeparts.mkaudioamplify(pipeline, head[rate], 1. / math.sqrt(pipeparts.audioresample_variance_gain(quality, max_rate, rate)))
-		head[rate] = pipeparts.mkcapsfilter(pipeline, pipeparts.mkinterpolator(pipeline, head[rate]), caps = "audio/x-raw, rate=%d" % rate)
+		head[rate] = pipeparts.mkcapsfilter(pipeline, pipeparts.mkinterpolator(pipeline, head[rate]), caps = "audio/x-raw, rate=%d" % max(min_rate, rate))
 
 	#
 	# return value is a dictionary of elements indexed by sample rate

gstlal-burst/python/fxtools/waveforms.py

@@ -105,7 +105,7 @@ class TemplateGenerator(object):
 		"""
 		return NotImplementedError
 
-	def generate_templates(self, rate, quadrature = True):
+	def generate_templates(self, rate, quadrature = True, sampling_rate = None):
 		"""
 		Generate all templates corresponding to a parameter range for a given sampling rate.
 		If quadrature is set, yield two templates corresponding to in-phase and quadrature components of the waveform.
@@ -170,17 +170,20 @@ class HalfSineGaussianGenerator(TemplateGenerator):
 		"""
 		return 0.5 * (q/(2.*numpy.pi*f)) * numpy.log(1./self.tolerance)
 
-	def generate_templates(self, rate, quadrature = True):
+	def generate_templates(self, rate, quadrature = True, sampling_rate = None):
 		"""
 		generate all half sine gaussian templates corresponding to a parameter range and template duration
 		for a given sampling rate.
 		"""
+		if not sampling_rate:
+			sampling_rate = rate
+
 		for template in self.parameter_grid[rate]:
 			if quadrature:
 				for phase in self.phases:
-					yield self.waveform(rate, phase, template['frequency'], template['q'])
+					yield self.waveform(sampling_rate, phase, template['frequency'], template['q'])
 			else:
-				yield self.waveform(rate, self.phases[0], template['frequency'], template['q'])
+				yield self.waveform(sampling_rate, self.phases[0], template['frequency'], template['q'])
 
 	def waveform(self, rate, phase, f, q):
 		"""