calibration_parts.py 92 KB
Newer Older
Maddie White's avatar
Maddie White committed
1 2
#!/usr/bin/env python
#
Aaron Viets's avatar
Aaron Viets committed
3
# Copyright (C) 2015 Madeline Wade, Aaron Viets
Maddie White's avatar
Maddie White committed
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

19 20
import os

Maddie White's avatar
Maddie White committed
21
from gstlal import pipeparts
Madeline Wade's avatar
Madeline Wade committed
22
import numpy
Maddie White's avatar
Maddie White committed
23

24 25 26 27 28 29 30
import gi
gi.require_version('Gst', '1.0')
from gi.repository import GObject
from gi.repository import Gst
GObject.threads_init()
Gst.init(None)

31 32 33 34
#
# Shortcut functions for common element combos/properties
#

35
def mkqueue(pipeline, head, length = 0, min_length = 0):
36 37 38
	if length < 0:
		return head
	else:
39
		return pipeparts.mkqueue(pipeline, head, max_size_time = int(1000000000 * length), max_size_buffers = 0, max_size_bytes = 0, min_threshold_time = int(1000000000 * min_length))
40

41
def mkcomplexqueue(pipeline, head, length = 0, min_length = 0):
42
	head = pipeparts.mktogglecomplex(pipeline, head)
43
	head = mkqueue(pipeline, head, length = length, min_length = min_length)
44 45
	head = pipeparts.mktogglecomplex(pipeline, head)
	return head
46 47 48 49 50 51 52

def mkcapsfiltersetter(pipeline, head, caps, **properties):
	# Make a capsfilter followed by a capssetter
	head = pipeparts.mkcapsfilter(pipeline, head, caps)
	#head = pipeparts.mkcapssetter(pipeline, head, caps, replace = True, **properties)
	return head

53 54 55 56 57 58 59 60 61 62
def mkinsertgap(pipeline, head, **properties):
	if "bad_data_intervals" in properties:
		# Make sure the array property bad-data-intervals is formatted correctly
		intervals = properties.pop("bad_data_intervals")
		if intervals is not None:
			bad_data_intervals = []
			for i in range(0, len(intervals)):
				bad_data_intervals.append(float(intervals[i]))
			properties["bad_data_intervals"] = bad_data_intervals
	return pipeparts.mkgeneric(pipeline, head, "lal_insertgap", **properties)
63

64 65 66 67 68
#def mkupsample(pipeline, head, new_caps):
#	head = pipeparts.mkgeneric(pipeline, head, "lal_constantupsample")
#	head = pipeparts.mkcapsfilter(pipeline, head, new_caps)
#	return head

69
def mkstockresample(pipeline, head, caps):
70 71
	if type(caps) is int:
		caps = "audio/x-raw,rate=%d,channel-mask=(bitmask)0x0" % caps
72 73 74
	head = pipeparts.mkresample(pipeline, head, quality = 9)
	head = pipeparts.mkcapsfilter(pipeline, head, caps)
	return head
75

76
def mkresample(pipeline, head, quality, zero_latency, caps):
77 78
	if type(caps) is int:
		caps = "audio/x-raw,rate=%d,channel-mask=(bitmask)0x0" % caps
79
	head = pipeparts.mkgeneric(pipeline, head, "lal_resample", quality = quality, zero_latency = zero_latency)
Madeline Wade's avatar
Madeline Wade committed
80
	head = pipeparts.mkcapsfilter(pipeline, head, caps)
Madeline Wade's avatar
Madeline Wade committed
81 82
	return head

83
def mkcomplexfirbank(pipeline, src, latency = None, fir_matrix = None, time_domain = None, block_stride = None):
84 85 86 87 88 89 90 91 92
	if fir_matrix is not None:
		# Make sure the fir matrix is formatted correctly
		matrix = []
		for i in range(0, len(fir_matrix)):
			firfilt = []
			for j in range(0, len(fir_matrix[i])):
				firfilt.append(float(fir_matrix[i][j]))
			matrix.append(firfilt)
		fir_matrix = matrix
93 94 95
	properties = dict((name, value) for name, value in zip(("latency", "fir_matrix", "time_domain", "block_stride"), (latency, fir_matrix, time_domain, block_stride)) if value is not None)
	return pipeparts.mkgeneric(pipeline, src, "lal_complexfirbank", **properties)

96
def mkcomplexfirbank2(pipeline, src, latency = None, fir_matrix = None, time_domain = None, block_stride = None):
97 98 99 100 101 102 103 104 105
	if fir_matrix is not None:
		# Make sure the fir matrix is formatted correctly
		matrix = []
		for i in range(0, len(fir_matrix)):
			firfilt = []
			for j in range(0, len(fir_matrix[i])):
				firfilt.append(float(fir_matrix[i][j]))
			matrix.append(firfilt)
		fir_matrix = matrix
106 107 108
	properties = dict((name, value) for name, value in zip(("latency", "fir_matrix", "time_domain", "block_stride"), (latency, fir_matrix, time_domain, block_stride)) if value is not None)
	return pipeparts.mkgeneric(pipeline, src, "lal_complexfirbank2", **properties)

109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
def mkfccupdate(pipeline, src, **properties):
	if "fir_matrix" in properties:
		# Make sure the fir matrix is formatted correctly
		matrix = properties.pop("fir_matrix")
		if matrix is not None:
			fir_matrix = []
			for i in range(0, len(matrix)):
				firfilt = []
				for j in range(0, len(matrix[i])):
					firfilt.append(float(matrix[i][j]))
				fir_matrix.append(firfilt)
			properties["fir_matrix"] = fir_matrix
	return pipeparts.mkgeneric(pipeline, src, "lal_fcc_update", **properties)

def mktransferfunction(pipeline, src, **properties):
	if "notch_frequencies" in properties:
		# Make sure the array property notch-frequencies is formatted correctly
		freqs = properties.pop("notch_frequencies")
		if freqs is not None:
			notch_frequencies = []
			for i in range(0, len(freqs)):
				notch_frequencies.append(float(freqs[i]))
			properties["notch_frequencies"] = notch_frequencies
	return pipeparts.mkgeneric(pipeline, src, "lal_transferfunction", **properties)

def mkadaptivefirfilt(pipeline, src, **properties):
	# Make sure each array property is formatted correctly
	if "static_filter" in properties:
		staticfilt = properties.pop("static_filter")
		if staticfilt is not None:
			static_filter = []
			for i in range(0, len(staticfilt)):
				static_filter.append(float(staticfilt[i]))
			properties["static_filter"] = static_filter
	if "static_zeros" in properties:
		staticz = properties.pop("static_zeros")
		if staticz is not None:
			static_zeros = []
			for i in range(0, len(staticz)):
				static_zeros.append(float(staticz[i]))
			properties["static_zeros"] = static_zeros
	if "static_poles" in properties:
		staticp = properties.pop("static_poles")
		if staticp is not None:
			static_poles = []
			for i in range(0, len(staticp)):
				static_poles.append(float(staticp[i]))
			properties["static_poles"] = static_poles
	return pipeparts.mkgeneric(pipeline, src, "lal_adaptivefirfilt", **properties)

159 160 161
def mkpow(pipeline, src, **properties):
	return pipeparts.mkgeneric(pipeline, src, "cpow", **properties)

162
def mkmultiplier(pipeline, srcs, sync = True, queue_length = [0]):
163
	elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync=sync, mix_mode="product")
Maddie White's avatar
Maddie White committed
164
	if srcs is not None:
165 166
		for i in range(0, len(srcs)):
			mkqueue(pipeline, srcs[i], length = queue_length[min(i, len(queue_length) - 1)]).link(elem)
167 168
	return elem

169
def mkadder(pipeline, srcs, sync = True, queue_length = [0]):
170 171
	elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync=sync)
	if srcs is not None:
172 173
		for i in range(0, len(srcs)):
			mkqueue(pipeline, srcs[i], length = queue_length[min(i, len(queue_length) - 1)]).link(elem)
174 175
	return elem

176 177
def mkgate(pipeline, src, control, threshold, queue_length = 0, **properties):
	elem = pipeparts.mkgate(pipeline, mkqueue(pipeline, src, length = queue_length), control = mkqueue(pipeline, control, length = queue_length), threshold = threshold, **properties)
178 179
	return elem

180
def mkinterleave(pipeline, srcs, complex_data = False, queue_length = [0]):
181 182
	complex_factor = 1 + int(complex_data)
	num_srcs = complex_factor * len(srcs)
183 184 185 186 187 188
	i = 0
	mixed_srcs = []
	for src in srcs:
		matrix = [numpy.zeros(num_srcs)]
		matrix[0][i] = 1
		mixed_srcs.append(pipeparts.mkmatrixmixer(pipeline, src, matrix=matrix))
189
		i += complex_factor
190
	elem = mkadder(pipeline, tuple(mixed_srcs), queue_length = queue_length)
191 192 193

	#chan1 = pipeparts.mkmatrixmixer(pipeline, src1, matrix=[[1,0]])
	#chan2 = pipeparts.mkmatrixmixer(pipeline, src2, matrix=[[0,1]])
194
	#elem = mkadder(pipeline, list_srcs(pipeline, chan1, chan2)) 
195

196 197 198 199
	#elem = pipeparts.mkgeneric(pipeline, None, "interleave")
	#if srcs is not None:
	#	for src in srcs:
	#		pipeparts.mkqueue(pipeline, src).link(elem)
200 201
	return elem

202 203
def mkdeinterleave(pipeline, src, num_channels, complex_data = False):
	complex_factor = 1 + int(complex_data)
204 205 206
	head = pipeparts.mktee(pipeline, src)
	streams = []
	for i in range(0, num_channels):
207
		matrix = numpy.zeros((num_channels, complex_factor))
208 209 210 211 212
		matrix[i][0] = 1.0
		streams.append(pipeparts.mkmatrixmixer(pipeline, head, matrix = matrix))

	return tuple(streams)

Maddie White's avatar
Maddie White committed
213

214 215 216 217 218 219 220 221 222 223 224
#
# Write a pipeline graph function
#

def write_graph(demux, pipeline, name):
	pipeparts.write_dump_dot(pipeline, "%s.%s" % (name, "PLAYING"), verbose = True)

#
# Common element combo functions
#

225
def hook_up(pipeline, demux, channel_name, instrument, buffer_length, element_name_suffix = ""):
226
	if channel_name.endswith("UNCERTAINTY"):
227
		head = mkinsertgap(pipeline, None, bad_data_intervals = [-1e35, -1e-35, 1e-35, 1e35], insert_gap = False, remove_gap = True, fill_discont = True, block_duration = int(1000000000 * buffer_length), replace_value = 1, name = "insertgap_%s%s" % (channel_name, element_name_suffix))
228
	else:
229
		head = mkinsertgap(pipeline, None, bad_data_intervals = [-1e35, -1e-35, 1e-35, 1e35], insert_gap = False, remove_gap = True, fill_discont = True, block_duration = int(1000000000 * buffer_length), replace_value = 0, name = "insertgap_%s%s" % (channel_name, element_name_suffix))
230
	pipeparts.src_deferred_link(demux, "%s:%s" % (instrument, channel_name), head.get_static_pad("sink"))
231

232 233
	return head

234
def caps_and_progress(pipeline, head, caps, progress_name):
235 236
	head = pipeparts.mkgeneric(pipeline, head, "lal_typecast")
	head = pipeparts.mkcapsfilter(pipeline, head, caps)
237
	head = pipeparts.mkprogressreport(pipeline, head, "progress_src_%s" % progress_name)
238

239 240 241 242 243 244 245
	return head


#
# Function to make a list of heads to pass to, i.e. the multiplier or adder
#

Maddie White's avatar
Maddie White committed
246 247 248
def list_srcs(pipeline, *args):
	out = []
	for src in args:
249
		out.append(src)
Maddie White's avatar
Maddie White committed
250 251
	return tuple(out)

252 253 254 255
#
# Various filtering functions
#

256
def demodulate(pipeline, head, freq, td, rate, filter_time, filter_latency, prefactor_real = 1.0, prefactor_imag = 0.0, freq_update = None):
257 258 259
	# demodulate input at a given frequency freq

	head = pipeparts.mkgeneric(pipeline, head, "lal_demodulate", line_frequency = freq, prefactor_real = prefactor_real, prefactor_imag = prefactor_imag)
260
	if type(freq_update) is list:
261 262 263
		freq_update[0].connect("notify::current-average", update_property_simple, head, "current_average", "line_frequency", 1)
		freq_update[1].connect("notify::current-average", update_property_simple, head, "current_average", "prefactor_real", 1)
		freq_update[2].connect("notify::current-average", update_property_simple, head, "current_average", "prefactor_imag", 1)
264
	elif freq_update is not None:
265
		freq_update.connect("notify::current-average", update_property_simple, head, "current_average", "line_frequency", 1)
266
	head = mkresample(pipeline, head, 4, filter_latency == 0.0, rate)
267 268 269
	if filter_latency != 0:
		# Remove the first several seconds of output, which depend on start time
		head = pipeparts.mkgeneric(pipeline, head, "lal_insertgap", chop_length = 7000000000)
270
	head = lowpass(pipeline, head, rate, length = filter_time, fcut = 0, filter_latency = filter_latency, td = td)
271 272 273

	return head

274
def remove_harmonics(pipeline, signal, f0, num_harmonics, f0_var, filter_latency, compute_rate = 16, rate_out = 16384):
275 276 277
	# remove any line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable
	# function argument caps must be complex caps

Aaron Viets's avatar
Aaron Viets committed
278
	filter_param = 0.0625
279 280 281
	head = pipeparts.mktee(pipeline, head)
	elem = pipeparts.mkgeneric(pipeline, None, "lal_adder", sync = True)
	mkqueue(pipeline, head).link(elem)
282 283
	for i in range(1, num_harmonics + 1):
		line = pipeparts.mkgeneric(pipeline, head, "lal_demodulate", line_frequency = i * f0)
284
		line = mkresample(pipeline, line, 4, filter_latency == 0, compute_rate)
Aaron Viets's avatar
Aaron Viets committed
285
		line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = filter_param / (f0_var * i), fcut = 0, filter_latency = filter_latency)
286
		line = mkresample(pipeline, line, 3, filter_latency == 0.0, rate_out)
287
		line = pipeparts.mkgeneric(pipeline, line, "lal_demodulate", line_frequency = -1.0 * i * f0, prefactor_real = -2.0)
288 289
		line = pipeparts.mkgeneric(pipeline, line, "creal")
		mkqueue(pipeline, line).link(elem)
290 291 292

	return elem

293
def remove_lines_with_witnesses(pipeline, signal, witnesses, freqs, freq_vars, freq_channels, filter_latency = 0, compute_rate = 16, rate_out = 16384, num_median = 2048, num_avg = 160, noisesub_gate_bit = None):
294 295 296
	# remove line(s) from a spectrum. filter length for demodulation (given in seconds) is adjustable
	# function argument caps must be complex caps

297 298
	# Re-format inputs if necessary
	if type(witnesses) is not list and type(witnesses) is not tuple and type(witnesses) is not numpy.ndarray:
299
		print("remove_lines_with_witnesses(): argument 3 should be type list.  Converting %s to list" % type(witnesses))
300 301
		witnesses = [[witnesses]]
	if type(freqs) is not list and type(freqs) is not tuple and type(freqs) is not numpy.ndarray:
302
		print("remove_lines_with_witnesses(): argument 4 should be type list.  Converting %s to list" % type(freqs))
303 304
		freqs = [[freqs]]
	if type(freq_vars) is not list and type(freq_vars) is not tuple and type(freq_vars) is not numpy.ndarray:
305
		print("remove_lines_with_witnesses(): argument 5 should be type list.  Converting %s to list" % type(freq_vars))
306 307
		freq_vars = [freq_vars]
	for i in range(0, len(witnesses) - len(freqs)):
308
		print("remove_lines_with_witnesses(): Warning: not enough elements in argument 4")
309 310
		freqs.append(freqs[-1])
	for i in range(0, len(witnesses) - len(freq_vars)):
311
		print("remove_lines_with_witnesses(): Warning: not enough elements in argument 5")
312 313
		freq_vars.append(freq_vars[-1])
	if len(freqs) > len(witnesses):
314
		print("remove_lines_with_witnesses(): Warning: too many elements in argument 4")
315 316
		freqs = freqs[:len(witnesses)]
	if len(freq_vars) > len(witnesses):
317
		print("remove_lines_with_witnesses(): Warning: too many elements in argument 5")
318 319 320
		freq_vars = freq_vars[:len(witnesses)]
	for i in range(0, len(witnesses)):
		if type(witnesses[i]) is not list and type(witnesses[i]) is not tuple and type(witnesses[i]) is not numpy.ndarray:
321
			print("remove_lines_with_witnesses(): argument 3 should be list of lists.  Converting %s to list" % type(witnesses[i]))
322 323
			witnesses[i] = [witnesses[i]]
		if type(freqs[i]) is not list and type(freqs[i]) is not tuple and type(freqs[i]) is not numpy.ndarray:
324
			print("remove_lines_with_witnesses(): argument 4 should be list of lists.  Converting %s to list" % type(freqs[i]))
325 326
			freqs[i] = [freqs[i]]

Aaron Viets's avatar
Aaron Viets committed
327
	filter_param = 0.0625
328
	downsample_quality = 4
Aaron Viets's avatar
Aaron Viets committed
329
	upsample_quality = 4
330
	resample_shift = 16.0 + 16.5
Aaron Viets's avatar
Aaron Viets committed
331
	zero_latency = filter_latency == 0.0
332

333
	for i in range(0, len(witnesses)):
334 335
		for j in range(0, len(witnesses[i])):
			witnesses[i][j] = pipeparts.mktee(pipeline, witnesses[i][j])
336
	signal = pipeparts.mktee(pipeline, signal)
337 338
	signal_minus_lines = [signal]

339 340 341 342 343 344 345 346 347 348 349 350 351
	for m in range(0, len(witnesses)):
		# If freqs[m][0] strays from its nominal value and there is a timestamp shift in the signal
		# (e.g., to achieve zero latency), we need to correct the phase in the reconstructed
		# signal. To do this, we measure the frequency in the witness and find the beat
		# frequency between that and the nominal frequency freqs[m][0].
		if filter_latency != 0.5 and freq_vars[m]:
			# The low-pass and resampling filters are not centered in time
			f0_measured = pipeparts.mkgeneric(pipeline, witnesses[m][0], "lal_trackfrequency", num_halfcycles = int(round((filter_param / freq_vars[m] / 2 + resample_shift / compute_rate) * freqs[m][0])))
			f0_measured = mkresample(pipeline, f0_measured, 3, zero_latency, compute_rate)
			f0_measured = pipeparts.mkgeneric(pipeline, f0_measured, "lal_smoothkappas", array_size = 1, avg_array_size = int(round((filter_param / freq_vars[m] / 2 * compute_rate + resample_shift) / 2)), default_kappa_re = 0, default_to_median = True, filter_latency = filter_latency)
			f0_beat_frequency = pipeparts.mkgeneric(pipeline, f0_measured, "lal_add_constant", value = -freqs[m][0])
			f0_beat_frequency = pipeparts.mktee(pipeline, f0_beat_frequency)

352
		for n in range(len(freqs[m])):
353
			# Length of low-pass filter
354
			filter_length = filter_param / (max(freq_vars[m], 0.003) * freqs[m][n] / freqs[m][0])
355 356 357 358
			filter_samples = int(filter_length * compute_rate) + (1 - int(filter_length * compute_rate) % 2)
			sample_shift = filter_samples / 2 - int((filter_samples - 1) * filter_latency + 0.5)
			# shift of timestamp relative to data
			time_shift = float(sample_shift) / compute_rate + zero_latency * resample_shift / compute_rate
359
			two_n_pi_delta_t = 2 * freqs[m][n] / freqs[m][0] * numpy.pi * time_shift
360 361 362 363 364 365 366 367 368 369

			# Only do this if we have to
			if filter_latency != 0.5 and freq_vars[m]:
				# Find phase shift due to timestamp shift for each harmonic
				phase_shift = pipeparts.mkmatrixmixer(pipeline, f0_beat_frequency, matrix=[[0, two_n_pi_delta_t]])
				phase_shift = pipeparts.mktogglecomplex(pipeline, phase_shift)
				phase_factor = pipeparts.mkgeneric(pipeline, phase_shift, "cexp")
				phase_factor = pipeparts.mktee(pipeline, phase_factor)

			# Find amplitude and phase of line in signal
370 371 372 373 374 375 376 377 378
			line_in_signal = pipeparts.mkgeneric(pipeline, signal, "lal_demodulate", line_frequency = freqs[m][n])
			# Connect to line frequency updater if given
			if freq_channels[m][n] is not None:
				if type(freq_channels[m][n]) is float:
					# It's a harmonic of the frequency in freq_channels[m][0]
					freq_channels[m][0].connect("notify::current-average", update_property_simple, line_in_signal, "current_average", "line_frequency", freq_channels[m][n])
				else:
					# The channel carries the correct frequency
					freq_channels[m][n].connect("notify::current-average", update_property_simple, line_in_signal, "current_average", "line_frequency", 1)
379 380 381 382 383
			line_in_signal = mkresample(pipeline, line_in_signal, downsample_quality, zero_latency, compute_rate)
			line_in_signal = lowpass(pipeline, line_in_signal, compute_rate, length = filter_length, fcut = 0, filter_latency = filter_latency)
			line_in_signal = pipeparts.mktee(pipeline, line_in_signal)

			# Make ones for use in matrix equation
384
			if m == 0 and n == 0:
385 386 387 388 389 390
				ones = pipeparts.mktee(pipeline, mkpow(pipeline, line_in_signal, exponent = 0.0))

			line_in_witnesses = []
			tfs_at_f = [None] * len(witnesses[m]) * (len(witnesses[m]) + 1)
			for i in range(0, len(witnesses[m])):
				# Find amplitude and phase of each harmonic in each witness channel
391 392 393 394 395 396 397 398 399
				line_in_witness = pipeparts.mkgeneric(pipeline, witnesses[m][i], "lal_demodulate", line_frequency = freqs[m][n])
				# Connect to line frequency updater if given
				if freq_channels[m][n] is not None:
					if type(freq_channels[m][n]) is float:
						# It's a harmonic of the frequency in freq_channels[m][0]
						freq_channels[m][0].connect("notify::current-average", update_property_simple, line_in_witness, "current_average", "line_frequency", freq_channels[m][n])
					else:
						# The channel carries the correct frequency
						freq_channels[m][n].connect("notify::current-average", update_property_simple, line_in_witness, "current_average", "line_frequency", 1)
400 401 402 403 404 405 406 407 408 409
				line_in_witness = mkresample(pipeline, line_in_witness, downsample_quality, zero_latency, compute_rate)
				line_in_witness = lowpass(pipeline, line_in_witness, compute_rate, length = filter_length, fcut = 0, filter_latency = filter_latency)
				line_in_witness = pipeparts.mktee(pipeline, line_in_witness)
				line_in_witnesses.append(line_in_witness)

				# Find transfer function between witness channel and signal at this frequency
				tf_at_f = complex_division(pipeline, line_in_signal, line_in_witness)

				# Remove worthless data from computation of transfer function if we can
				if noisesub_gate_bit is not None:
410
					tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples), name = "powerlines_gate_%d_%d_%d" % (m, n, i))
411
				tfs_at_f[i] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0.0, default_kappa_im = 0.0, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
412 413 414 415 416 417 418 419 420 421
				tfs_at_f[(i + 1) * len(witnesses[m]) + i] = ones

			for i in range(0, len(witnesses[m])):
				for j in range(0, len(witnesses[m])):
					if(i != j):
						# Find transfer function between 2 witness channels at this frequency
						tf_at_f = complex_division(pipeline, line_in_witnesses[j], line_in_witnesses[i])

						# Remove worthless data from computation of transfer function if we can
						if noisesub_gate_bit is not None:
422
							tf_at_f = mkgate(pipeline, tf_at_f, noisesub_gate_bit, 1, attack_length = -((1.0 - filter_latency) * filter_samples), name = "powerlines_gate_%d_%d_%d_%d" % (m, n, i, j))
423
						tfs_at_f[(i + 1) * len(witnesses[m]) + j] = pipeparts.mkgeneric(pipeline, tf_at_f, "lal_smoothkappas", default_kappa_re = 0.0, default_kappa_im = 0.0, array_size = num_median, avg_array_size = num_avg, default_to_median = True, filter_latency = filter_latency)
424 425 426 427 428 429 430 431 432 433 434 435

			tfs_at_f = mkinterleave(pipeline, tfs_at_f, complex_data = True)
			tfs_at_f = pipeparts.mkgeneric(pipeline, tfs_at_f, "lal_matrixsolver")
			tfs_at_f = mkdeinterleave(pipeline, tfs_at_f, len(witnesses[m]), complex_data = True)

			for i in range(0, len(witnesses[m])):
				# Use gated, averaged transfer function to reconstruct the sinusoid as it appears in the signal from the witness channel
				if filter_latency == 0.5 or not freq_vars[m]:
					reconstructed_line_in_signal = mkmultiplier(pipeline, list_srcs(pipeline, tfs_at_f[i], line_in_witnesses[i]))
				else:
					reconstructed_line_in_signal = mkmultiplier(pipeline, list_srcs(pipeline, tfs_at_f[i], line_in_witnesses[i], phase_factor))
				reconstructed_line_in_signal = mkresample(pipeline, reconstructed_line_in_signal, upsample_quality, zero_latency, rate_out)
436 437 438 439 440 441 442 443 444
				reconstructed_line_in_signal = pipeparts.mkgeneric(pipeline, reconstructed_line_in_signal, "lal_demodulate", line_frequency = -1.0 * freqs[m][n], prefactor_real = -2.0)
				# Connect to line frequency updater if given
				if freq_channels[m][n] is not None:
					if type(freq_channels[m][n]) is float:
						# It's a harmonic of the frequency in freq_channels[m][0]
						freq_channels[m][0].connect("notify::current-average", update_property_simple, reconstructed_line_in_signal, "current_average", "line_frequency", -1.0 * freq_channels[m][n])
					else:
						# The channel carries the correct frequency
						freq_channels[m][n].connect("notify::current-average", update_property_simple, reconstructed_line_in_signal, "current_average", "line_frequency", -1.0)
445 446 447
				reconstructed_line_in_signal = pipeparts.mkgeneric(pipeline, reconstructed_line_in_signal, "creal")

				signal_minus_lines.append(reconstructed_line_in_signal)
448 449

	clean_signal = mkadder(pipeline, tuple(signal_minus_lines))
450

451
	return clean_signal
452

453
def removeDC(pipeline, head, rate):
454
	head = pipeparts.mktee(pipeline, head)
455
	DC = mkresample(pipeline, head, 4, True, 16)
456
	#DC = pipeparts.mkgeneric(pipeline, DC, "lal_smoothkappas", default_kappa_re = 0, array_size = 1, avg_array_size = 64)
457
	DC = mkresample(pipeline, DC, 4, True, rate)
458 459
	DC = pipeparts.mkaudioamplify(pipeline, DC, -1)

460
	return mkadder(pipeline, list_srcs(pipeline, head, DC))
461

462
def lowpass(pipeline, head, rate, length = 1.0, fcut = 500, filter_latency = 0.5, td = True):
463
	length = int(length * rate)
464 465 466 467 468 469 470 471 472
	if not length % 2:
		length += 1 # Make sure the filter length is odd

	# Compute a low-pass filter.
	lowpass = numpy.sinc(2 * float(fcut) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Now apply the filter
Aaron Viets's avatar
Aaron Viets committed
473
	return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency + 0.5), fir_matrix = [lowpass], time_domain = td)
474

475
def highpass(pipeline, head, rate, length = 1.0, fcut = 9.0, filter_latency = 0.5, td = True):
476
	length = int(length * rate)
477 478 479 480 481 482 483 484 485 486
	if not length % 2:
		length += 1 # Make sure the filter length is odd

	# Compute a low-pass filter.
	lowpass = numpy.sinc(2 * float(fcut) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Create a high-pass filter from the low-pass filter through spectral inversion.
	highpass = -lowpass
487
	highpass[int((length - 1) / 2)] += 1
488 489

	# Now apply the filter
Aaron Viets's avatar
Aaron Viets committed
490
	return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * filter_latency + 0.5), fir_matrix = [highpass], time_domain = td)
491

492
def bandpass(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filter_latency = 0.5, td = True):
493
	length = int(length * rate / 2)
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
	if not length % 2:
		length += 1 # Make sure the filter length is odd

	# Compute a temporary low-pass filter.
	lowpass = numpy.sinc(2 * float(f_low) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Create the high-pass filter from the low-pass filter through spectral inversion.
	highpass = -lowpass
	highpass[(length - 1) / 2] += 1

	# Compute the low-pass filter.
	lowpass = numpy.sinc(2 * float(f_high) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Convolve the high-pass and low-pass filters to make a band-pass filter
	bandpass = numpy.convolve(highpass, lowpass)

	# Now apply the filter
Aaron Viets's avatar
Aaron Viets committed
515
	return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency + 0.5), fir_matrix = [bandpass], time_domain = td)
516

517
def bandstop(pipeline, head, rate, length = 1.0, f_low = 100, f_high = 400, filter_latency = 0.5, td = True):
518
	length = int(length * rate / 2)
519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543
	if not length % 2:
		length += 1 # Make sure the filter length is odd

	# Compute a temporary low-pass filter.
	lowpass = numpy.sinc(2 * float(f_low) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Create a high-pass filter from the low-pass filter through spectral inversion.
	highpass = -lowpass
	highpass[(length - 1) / 2] += 1

	# Compute a low-pass filter.
	lowpass = numpy.sinc(2 * float(f_high) / rate * (numpy.arange(length) - (length - 1) / 2))
	lowpass *= numpy.blackman(length)
	lowpass /= numpy.sum(lowpass)

	# Convolve the high-pass and low-pass filters to make a temporary band-pass filter
	bandpass = numpy.convolve(highpass, lowpass)

	# Create a band-stop filter from the band-pass filter through spectral inversion.
	bandstop = -bandpass
	bandstop[length - 1] += 1

	# Now apply the filter
Aaron Viets's avatar
Aaron Viets committed
544
	return mkcomplexfirbank(pipeline, head, latency = int((length - 1) * 2 * filter_latency + 0.5), fir_matrix = [bandstop], time_domain = td)
545

546
def linear_phase_filter(pipeline, head, shift_samples, num_samples = 256, gain = 1.0, filter_update = None, sample_rate = 2048, update_samples = 320, average_samples = 1, phase_measurement_frequency = 100, taper_length = 320, kernel_endtime = None, filter_timeshift = 0):
547

548 549 550 551 552
	# Apply a linear-phase filter to shift timestamps.  shift_samples is the number
	# of samples of timestamp shift.  It need not be an integer.  A positive value
	# advances the output data relative to the timestamps, and a negative value
	# delays the output.

553 554
	# Compute filter using odd filter length
	odd_num_samples = int(num_samples) - (1 - int(num_samples) % 2)
555

556
	filter_latency_samples = int(num_samples / 2) + int(numpy.floor(shift_samples))
557 558 559
	fractional_shift_samples = shift_samples % 1

	# Make a filter using a sinc table, slightly shifted relative to the samples
560
	sinc_arg = numpy.arange(-int(odd_num_samples / 2), 1 + int(odd_num_samples / 2)) + fractional_shift_samples
561 562
	sinc_filter = numpy.sinc(sinc_arg)
	# Apply a Blackman window
563
	sinc_filter *= numpy.blackman(odd_num_samples)
564
	# Normalize the filter
565
	sinc_filter *= gain / numpy.sum(sinc_filter)
566 567 568
	# In case filter length is actually even
	if not int(num_samples) % 2:
		sinc_filter = numpy.insert(sinc_filter, 0, 0.0)
569 570

	# Filter the data
571 572 573 574 575 576 577
	if filter_update is None:
		# Static filter
		head =  mkcomplexfirbank(pipeline, head, latency = filter_latency_samples, fir_matrix = [sinc_filter[::-1]], time_domain = True)
	else:
		# Filter gets updated with variable time delay and gain
		if kernel_endtime is None:
			# Update filter as soon as new filter is available, and do it with minimal latency
578 579
			head = pipeparts.mkgeneric(pipeline, head, "lal_tdwhiten", kernel = sinc_filter[::-1], latency = filter_latency_samples, taper_length = taper_length)
			filter_update = mkadaptivefirfilt(pipeline, filter_update, variable_filter_length = num_samples, adaptive_filter_length = num_samples, update_samples = update_samples, average_samples = average_samples, filter_sample_rate = sample_rate, phase_measurement_frequency = phase_measurement_frequency, tukey_param = 0.5)
580 581 582
			filter_update.connect("notify::adaptive-filter", update_filter, head, "adaptive_filter", "kernel")
		else:
			# Update filters at specified timestamps to ensure reproducibility
583 584
			head = pipeparts.mkgeneric(pipeline, mkqueue(pipeline, head), "lal_tdwhiten", kernel = sinc_filter[::-1], latency = filter_latency_samples, taper_length = taper_length, kernel_endtime = kernel_endtime)
			filter_update = mkadaptivefirfilt(pipeline, filter_update, variable_filter_length = num_samples, adaptive_filter_length = num_samples, update_samples = update_samples, average_samples = average_samples, filter_sample_rate = sample_rate, phase_measurement_frequency = phase_measurement_frequency, filter_timeshift = filter_timeshift, tukey_param = 0.5)
585
			filter_update.connect("notify::adaptive-filter", update_filter, head, "adaptive_filter", "kernel")
586
			filter_update.connect("notify::filter-endtime", update_property_simple, head, "filter_endtime", "kernel_endtime", 1)
587
	return head
588

589
def compute_rms(pipeline, head, rate, average_time, f_min = None, f_max = None, filter_latency = 0.5, rate_out = 16, td = True):
590 591
	# Find the root mean square amplitude of a signal between two frequencies
	# Downsample to save computational cost
592
	head = mkresample(pipeline, head, 4, filter_latency == 0.0, rate)
593 594 595

	# Remove any frequency content we don't care about
	if (f_min is not None) and (f_max is not None):
596
		head = bandpass(pipeline, head, rate, f_low = f_min, f_high = f_max, filter_latency = filter_latency, td = td)
597
	elif f_min is not None:
598
		head = highpass(pipeline, head, rate, fcut = f_min, filter_latency = filter_latency, td = td)
599
	elif f_max is not None:
600
		head = lowpass(pipeline, head, rate, fcut = f_max, filter_latency = filter_latency, td = td)
601 602

	# Square it
603
	head = mkpow(pipeline, head, exponent = 2.0)
604 605

	# Downsample again to save computational cost
606
	head = mkresample(pipeline, head, 4, filter_latency == 0.0, rate_out)
607 608

	# Compute running average
609
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = 0.0, array_size = 1, avg_array_size = average_time * rate_out, filter_latency = filter_latency)
610

611 612 613
	# Take the square root
	head = mkpow(pipeline, head, exponent = 0.5)

614 615
	return head

616 617 618 619
#
# Calibration factor related functions
#

620
def smooth_kappas_no_coherence(pipeline, head, var, expected, N, Nav, default_to_median, filter_latency):
621
	# Find median of calibration factors array with size N and smooth out medians with an average over Nav samples
622
	# Use the maximum_offset_re property to determine whether input kappas are good or not
623
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", maximum_offset_re = var, default_kappa_re = expected, array_size = N, avg_array_size = Nav, default_to_median = default_to_median, filter_latency = filter_latency)
624
	return head
625

626
def smooth_complex_kappas_no_coherence(pipeline, head, real_var, imag_var, real_expected, imag_expected, N, Nav, default_to_median, filter_latency):
627 628
	# Find median of complex calibration factors array with size N, split into real and imaginary parts, and smooth out medians with an average over Nav samples
	# Use the maximum_offset_re and maximum_offset_im properties to determine whether input kappas are good or not
629
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", maximum_offset_re = real_var, maximum_offset_im = imag_var, default_kappa_re = real_expected, default_kappa_im = imag_expected, array_size = N, avg_array_size = Nav, default_to_median = default_to_median, filter_latency = filter_latency)
630
	return head
631

632
def smooth_kappas(pipeline, head, expected, N, Nav, default_to_median, filter_latency):
633 634
	# Find median of calibration factors array with size N and smooth out medians with an average over Nav samples
	# Assume input was previously gated with coherence uncertainty to determine if input kappas are good or not
635
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = expected, array_size = N, avg_array_size = Nav, default_to_median = default_to_median, filter_latency = filter_latency)
636 637
	return head

638
def smooth_complex_kappas(pipeline, head, real_expected, imag_expected, N, Nav, default_to_median, filter_latency):
639 640 641
	# Find median of complex calibration factors array with size N and smooth out medians with an average over Nav samples
	# Assume input was previously gated with coherence uncertainty to determine if input kappas are good or not

642
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = real_expected, default_kappa_im = imag_expected, array_size = N, avg_array_size = Nav, default_to_median = default_to_median, filter_latency = filter_latency)
643
	return head
644

645
def track_bad_kappas_no_coherence(pipeline, head, var, expected, N, Nav, default_to_median, filter_latency):
646
	# Produce output of 1's or 0's that correspond to median not corrupted (1) or corrupted (0) based on whether median of input array is defualt value.
647
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", maximum_offset_re = var, default_kappa_re = expected, array_size = N, avg_array_size = Nav if default_to_median else 1, track_bad_kappa = True, default_to_median = default_to_median, filter_latency = filter_latency)
648 649
	return head

650
def track_bad_complex_kappas_no_coherence(pipeline, head, real_var, imag_var, real_expected, imag_expected, N, Nav, default_to_median, filter_latency):
651
	# Produce output of 1's or 0's that correspond to median not corrupted (1) or corrupted (0) based on whether median of input array is defualt value.
652
	# Real and imaginary parts are done separately (outputs of lal_smoothkappas can be 1+i, 1, i, or 0)
653
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", maximum_offset_re = real_var, maximum_offset_im = imag_var, default_kappa_re = real_expected, default_kappa_im = imag_expected, array_size = N, avg_array_size = Nav if default_to_median else 1, track_bad_kappa = True, default_to_median = default_to_median, filter_latency = filter_latency)
654 655 656
	re, im = split_into_real(pipeline, head)
	return re, im

657
def track_bad_kappas(pipeline, head, expected, N, Nav, default_to_median, filter_latency):
658
	# Produce output of 1's or 0's that correspond to median not corrupted (1) or corrupted (0) based on whether median of input array is defualt value.
659
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = expected, array_size = N, avg_array_size = Nav if default_to_median else 1, track_bad_kappa = True, default_to_median = default_to_median, filter_latency = filter_latency)
660 661
	return head

662
def track_bad_complex_kappas(pipeline, head, real_expected, imag_expected, N, Nav, default_to_median, filter_latency):
663
	# Produce output of 1's or 0's that correspond to median not corrupted (1) or corrupted (0) based on whether median of input array is defualt value.
664
	# Real and imaginary parts are done separately (outputs of lal_smoothkappas can be 1+i, 1, i, or 0)
665

666
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", default_kappa_re = real_expected, default_kappa_im = imag_expected, array_size = N, avg_array_size = Nav if default_to_median else 1, track_bad_kappa = True, default_to_median = default_to_median, filter_latency = filter_latency)
667
	re, im = split_into_real(pipeline, head)
668
	return re, im
669

670
def smooth_kappas_no_coherence_test(pipeline, head, var, expected, N, Nav, default_to_median, filter_latency):
671 672 673
	# Find median of calibration factors array with size N and smooth out medians with an average over Nav samples
	head = pipeparts.mktee(pipeline, head)
	pipeparts.mknxydumpsink(pipeline, head, "raw_kappatst.txt")
674
	head = pipeparts.mkgeneric(pipeline, head, "lal_smoothkappas", maximum_offset_re = var, default_kappa_re = expected, array_size = N, avg_array_size = Nav, default_to_median = default_to_median, filter_latency = filter_latency)
675 676 677
	head = pipeparts.mktee(pipeline, head)
	pipeparts.mknxydumpsink(pipeline, head, "smooth_kappatst.txt")
	return head
678

679 680 681 682 683 684 685 686 687 688
def compute_kappa_bits(pipeline, smoothR, smoothI, expected_real, expected_imag, real_ok_var, imag_ok_var, median_samples, avg_samples, status_out_smooth = 1, starting_rate=16, ending_rate=16):

	# Compensate for digital error in the running average
	expected_real_sum = 0.0
	expected_imag_sum = 0.0
	for i in range(0, avg_samples):
		expected_real_sum = expected_real_sum + expected_real
		expected_imag_sum = expected_imag_sum + expected_imag
	expected_real = expected_real_sum / avg_samples
	expected_imag = expected_imag_sum / avg_samples
689

690
	if type(real_ok_var) is list:
691
		smoothRInRange = mkinsertgap(pipeline, smoothR, bad_data_intervals = [real_ok_var[0], expected_real, expected_real, real_ok_var[1]], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
692
	else:
693
		smoothRInRange = mkinsertgap(pipeline, smoothR, bad_data_intervals = [expected_real - real_ok_var, expected_real, expected_real, expected_real + real_ok_var], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
694
	smoothRInRange = pipeparts.mkbitvectorgen(pipeline, smoothRInRange, nongap_is_control = True, bit_vector = status_out_smooth)
695 696 697 698
	smoothRInRange = pipeparts.mkcapsfilter(pipeline, smoothRInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % starting_rate)
	if starting_rate != ending_rate:
		smoothRInRange = pipeparts.mkgeneric(pipeline, smoothRInRange, "lal_logicalundersample", required_on = status_out_smooth, status_out = status_out_smooth)
		smoothRInRange = pipeparts.mkcapsfilter(pipeline, smoothRInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % ending_rate)
699 700
	smoothRInRangetee = pipeparts.mktee(pipeline, smoothRInRange)

701
	if type(imag_ok_var) is list:
702
		smoothIInRange = mkinsertgap(pipeline, smoothI, bad_data_intervals = [imag_ok_var[0], expected_imag, expected_imag, imag_ok_var[1]], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
703
	else:
704
		smoothIInRange = mkinsertgap(pipeline, smoothI, bad_data_intervals = [expected_imag - imag_ok_var, expected_imag, expected_imag, expected_imag + imag_ok_var], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
705
	smoothIInRange = pipeparts.mkbitvectorgen(pipeline, smoothIInRange, nongap_is_control = True, bit_vector = status_out_smooth)
706 707 708 709
	smoothIInRange = pipeparts.mkcapsfilter(pipeline, smoothIInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % starting_rate)
	if starting_rate != ending_rate:
		smoothIInRange = pipeparts.mkgeneric(pipeline, smoothIInRange, "lal_logicalundersample", required_on = status_out_smooth, status_out = status_out_smooth)
		smoothIInRange = pipeparts.mkcapsfilter(pipeline, smoothIInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % ending_rate)
710

711
	min_samples = int(median_samples / 2) + avg_samples
712
	smoothInRange = mkgate(pipeline, smoothRInRangetee, mkadder(pipeline, list_srcs(pipeline, smoothIInRange, smoothRInRangetee)), status_out_smooth * 2, attack_length = -min_samples)
713
	smoothInRange = pipeparts.mkbitvectorgen(pipeline, smoothInRange, nongap_is_control = True, bit_vector = status_out_smooth)
714
	smoothInRange = pipeparts.mkcapsfilter(pipeline, smoothInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % ending_rate)
715

716
	return smoothInRange
717

718 719 720 721 722 723 724
def compute_kappa_bits_only_real(pipeline, smooth, expected, ok_var, median_samples, avg_samples, status_out_smooth = 1, starting_rate=16, ending_rate=16):

	# Compensate for digital error in the running average
	expected_sum = 0.0
	for i in range(0, avg_samples):
		expected_sum = expected_sum + expected
	expected = expected_sum / avg_samples
725

726
	if type(ok_var) is list:
727
		smoothInRange = mkinsertgap(pipeline, smooth, bad_data_intervals = [ok_var[0], expected, expected, ok_var[1]], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
728
	else:
729
		smoothInRange = mkinsertgap(pipeline, smooth, bad_data_intervals = [expected - ok_var, expected, expected, expected + ok_var], insert_gap = True, remove_gap = False, replace_value = 0, fill_discont = True, block_duration = Gst.SECOND)
730
	smoothInRange = pipeparts.mkbitvectorgen(pipeline, smoothInRange, nongap_is_control = True, bit_vector = status_out_smooth)
731 732 733 734
	smoothInRange = pipeparts.mkcapsfilter(pipeline, smoothInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % starting_rate)
	if starting_rate != ending_rate:
		smoothInRange = pipeparts.mkgeneric(pipeline, smoothInRange, "lal_logicalundersample", required_on = status_out_smooth, status_out = status_out_smooth)
		smoothInRange = pipeparts.mkcapsfilter(pipeline, smoothInRange, "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % ending_rate)
735
	smoothInRangetee = pipeparts.mktee(pipeline, smoothInRange)
736
	min_samples = int(median_samples / 2) + avg_samples
737
	smoothInRange = mkgate(pipeline, smoothInRangetee, smoothInRangetee, status_out_smooth, attack_length = -min_samples)
738
	smoothInRange = pipeparts.mkbitvectorgen(pipeline, smoothInRange, nongap_is_control = True, bit_vector = status_out_smooth)
739

740
	return smoothInRange
741

742
def merge_into_complex(pipeline, real, imag):
743
	# Merge real and imag into one complex channel with complex caps
744
	head = mkinterleave(pipeline, list_srcs(pipeline, real, imag))
745
	head = pipeparts.mktogglecomplex(pipeline, head)
746 747
	return head

748
def split_into_real(pipeline, complex_chan):
749
	# split complex channel with complex caps into two channels (real and imag) with real caps
750

751 752 753
	complex_chan = pipeparts.mktee(pipeline, complex_chan)
	real = pipeparts.mkgeneric(pipeline, complex_chan, "creal")
	imag = pipeparts.mkgeneric(pipeline, complex_chan, "cimag")
754 755 756 757 758 759 760

#	elem = pipeparts.mkgeneric(pipeline, elem, "deinterleave", keep_positions=True)
#	real = pipeparts.mkgeneric(pipeline, None, "identity")
#	pipeparts.src_deferred_link(elem, "src_0", real.get_static_pad("sink"))
#	imag = pipeparts.mkgeneric(pipeline, None, "identity")
#	pipeparts.src_deferred_link(elem, "src_1", imag.get_static_pad("sink"))

761 762
	return real, imag

763
def complex_audioamplify(pipeline, chan, WR, WI):
764
	# Multiply a complex channel chan by a complex number WR+I WI
765 766
	# Re[out] = -chanI*WI + chanR*WR
	# Im[out] = chanR*WI + chanI*WR
Maddie White's avatar
Maddie White committed
767

768 769 770
	head = pipeparts.mktogglecomplex(pipeline, chan)
	head = pipeparts.mkmatrixmixer(pipeline, head, matrix=[[WR, WI],[-WI, WR]])
	head = pipeparts.mktogglecomplex(pipeline, head)