From 4691beb9d41059bbb12f1d676911f7a845cb90b2 Mon Sep 17 00:00:00 2001 From: Madeline Wade <madeline.wade@ligo.org> Date: Fri, 29 Sep 2017 07:03:49 -0700 Subject: [PATCH] Continued work on gstlal_clean_strain Adding in lal_fcc_update element that computes a new cavity pole filter and sends out message when new filter is computed --- gstlal-calibration/bin/Makefile.am | 3 +- gstlal-calibration/bin/gstlal_clean_strain | 1092 ++--------------- gstlal-calibration/gst/lal/Makefile.am | 3 +- gstlal-calibration/gst/lal/gstlal_fccupdate.c | 680 ++++++++++ gstlal-calibration/gst/lal/gstlal_fccupdate.h | 111 ++ .../gst/lal/gstlalcalibration.c | 2 + 6 files changed, 918 insertions(+), 973 deletions(-) create mode 100644 gstlal-calibration/gst/lal/gstlal_fccupdate.c create mode 100644 gstlal-calibration/gst/lal/gstlal_fccupdate.h diff --git a/gstlal-calibration/bin/Makefile.am b/gstlal-calibration/bin/Makefile.am index 93657f56ca..c5fc06e9c9 100644 --- a/gstlal-calibration/bin/Makefile.am +++ b/gstlal-calibration/bin/Makefile.am @@ -1,2 +1,3 @@ dist_bin_SCRIPTS = \ - gstlal_compute_strain + gstlal_compute_strain \ + gstlal_clean_strain diff --git a/gstlal-calibration/bin/gstlal_clean_strain b/gstlal-calibration/bin/gstlal_clean_strain index 54d6508d95..ccfbff058f 100644 --- a/gstlal-calibration/bin/gstlal_clean_strain +++ b/gstlal-calibration/bin/gstlal_clean_strain @@ -25,7 +25,6 @@ This pipeline is designed to perform various cleaning operations on h(t). Curren """ import sys -import os import numpy import time import resource @@ -92,12 +91,14 @@ parser.add_option("--gps-start-time", metavar = "seconds", help = "Set the start parser.add_option("--gps-end-time", metavar = "seconds", help = "Set the end time of the segment to analyze in GPS seconds. This is required iff --data-source=frames") parser.add_option("--wings", metavar = "seconds", type = "int", help = "Number of seconds to trim off of the beginning and end of the output. Should only be used if --data-source=frames.") parser.add_option("--do-file-checksum", action = "store_true", help = "Set this option to turn on file checksum in the demuxer.") -parser.add_option("--ifo", metavar = "name", help = "Name of the IFO to be calibrated.") +parser.add_option("--ifo", metavar = "name", help = "Name of the IFO strain to be cleaned.") parser.add_option("--raw-shared-memory-partition", metavar = "name", help = "Set the name of the shared memory partition to read from for raw data. This is required iff --data-source=lvshm.") parser.add_option("--hoft-shared-memory-partition", metavar = "name", help = "Set the name of the shared memory partition to read from for h(t) data. This is required iff --data-source=lvshm.") -parser.add_option("--raw-frame-segments-file", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load frame segments for raw data. This is required iff --data-source=frames") -parser.add_option("--hoft-frame-segments-name", metavar = "name", help = "Set the name of the segments to extract from the segment tables for h(t) data. This is required iff --frame-segments-file is given") -parser.add_option("--hoft-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate at which to generate strain data. This should be less than or equal to the sample rate of the error and control signal channels. (Default = 16384 Hz)") +parser.add_option("--raw-frame-segments-file", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load frame segments for raw data. This can only be used if --data-source=frames") +parser.add_option("--hoft-frame-segments-file", metavar = "filename", help = "Set the name of the LIGO light-weight XML file from which to load frame segments for h(t) data. This can only be used if --data-source=frames") +parser.add_option("--hoft-frame-segments-name", metavar = "name", help = "Set the name of the segments to extract from the segment tables for h(t) data. This can only be used if --frame-segments-file is given") +parser.add_option("--raw-frame-segments-name", metavar = "name", help = "Set the name of the segments to extract from the segment tables for raw data. This can only be used if --frame-segments-file is given") +parser.add_option("--hoft-sample-rate", metavar = "Hz", default = 16384, type = "int", help = "Sample rate of the h(t) channel. (Default = 16384 Hz)") parser.add_option("--buffer-length", metavar = "seconds", type = float, default = 1.0, help = "Set the length in seconds of buffers to be used in the pipeline (Default = 1.0)") parser.add_option("--frame-duration", metavar = "seconds", type = "int", default = 4, help = "Set the number of seconds for each frame. (Default = 4)") parser.add_option("--frames-per-file", metavar = "count", type = "int", default = 1, help = "Set the number of frames per frame file. (Default = 1)") @@ -121,10 +122,8 @@ parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose ( # These are options specific to the calibration procedure parser.add_option("--filters-file", metavar="filename", help = "Name of file containing filters (in npz format)") parser.add_option("--factors-from-filters-file", action = "store_true", help = "Compute the calibration factors from reference values contained in the filters file instead of from EPICS channels.") -parser.add_option("--exc-channel-name", metavar = "name", default = "CAL-CS_LINE_SUM_DQ", help = "Set the name of the excitation channel. This is only necessary when the calibration factors computation is turned on, which is the default behavior. (Default = CAL-CS_LINE_SUM_DQ)") -parser.add_option("--tst-exc-channel-name", metavar = "name", default = "SUS-ETMY_L3_CAL_LINE_OUT_DQ", help = "Set the name of the TST excitation channel. This is only necessary when the \kappa_tst factors computation is turned on, which is the default behavior. (Default = SUS-ETMY_L3_CAL_LINE_OUT_DQ)") -parser.add_option("--pcal-channel-name", metavar = "name", default = "CAL-PCALY_RX_PD_OUT_DQ", help = "Set the name of the PCal channel used for calculating the calibration factors. (Default = CAL-PCALY_RX_PD_OUT_DQ)") -parser.add_option("--dewhitening", action = "store_true", help = "Dewhitening should be used on the relevant channels, since the incoming channels are whitened and single precision.") +parser.add_option("--tst-exc-channel-name", metavar = "name", default = "SUS-ETMY_L3_CAL_LINE_OUT_DQ", help = "Set the name of the TST excitation channel. (Default = SUS-ETMY_L3_CAL_LINE_OUT_DQ)") +parser.add_option("--pcal-channel-name", metavar = "name", default = "CAL-PCALY_RX_PD_OUT_DQ", help = "Set the name of the PCal channel used. (Default = CAL-PCALY_RX_PD_OUT_DQ)") parser.add_option("--low-latency", action = "store_true", help = "Run the pipeline in low-latency mode. This uses minimal queueing. Otherwise, maximal queueing is used to prevent the pipeline from locking up.") parser.add_option("--remove-callines", action = "store_true", help = "Remove calibration lines at known freqencies from h(t) using software.") parser.add_option("--remove-powerlines", action = "store_true", help = "Remove 60 Hz spectral lines and some harmonics caused by power lines using witness channel PEM-EY_MAINSMON_EBAY_1_DQ.") @@ -147,37 +146,14 @@ parser.add_option("--remove-length-control", action = "store_true", help = "Remo parser.add_option("--lsc-srcl-channel-name", metavar = "name", default = "LSC-SRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-SRCL_IN1_DQ)") parser.add_option("--lsc-mich-channel-name", metavar = "name", default = "LSC-MICH_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-MICH_IN1_DQ)") parser.add_option("--lsc-prcl-channel-name", metavar = "name", default = "LSC-PRCL_IN1_DQ", help = "Set the name of one of the channels used as input from the LSC to remove length control noise. (Default = LSC-PRCL_IN1_DQ)") +parser.add_option("--kappa-tst-channel-name", metavar = "name", default = "GDS-CALIB_KAPPA_TST_REAL", help = "Set the name of the channel to be used for kappa_tst. (Default = GDS-CALIB_KAPPA_TST_REAL)") +parser.add_option("--hoft-channel-name", metavar = "name", default = "GDS-CALIB_STRAIN", help = "Set the name of the channel to be used for h(t). (Default = GDS-CALIB_STRAIN)") +parser.add_option("--kappa-sample-rate", metavar = "Hz", default = 16, type = "int", help = "Sample rate of the kappa_tst channel (Default = 16 Hz)") # These are all options related to the reference channels used in the calibration factors computation parser.add_option("--ref-channels-sr", metavar = "Hz", default = 16, help = "Set the sample rate for the reference model channels used in the calibration factors calculation. (Default = 16 Hz)") -parser.add_option("--EP4-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the ESD line used for the \kappa_a and \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_TST_REAL)") -parser.add_option("--EP5-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the ESD line used for the \kappa_a calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_REAL)") -parser.add_option("--EP3-real", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL", help = "Set the name of the channel containing the real part of 1/A_pu at the ESD line used for the \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_REAL)") -parser.add_option("--EP4-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the ESD line used for the \kappa_a and \kappa_pu calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_TST_IMAG") -parser.add_option("--EP5-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the ESD line used for the \kappa_A calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_IMAG") -parser.add_option("--EP3-imag", metavar = "name", default = "CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG", help = "Set the name of the channel containing the imaginary part of 1/A_pu at the ESD line used for the \kappa_PU calculation. (Default = CAL-CS_TDEP_DARM_LINE1_REF_A_USUM_INV_IMAG") -parser.add_option("--EP2-real", metavar = "name", default = "CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL", help = "Set the name of the channel containing the real part of the factors used to compute A(f_ctrl). (Default = CAL-CS_TDEP_REF_CLGRATIO_CTRL_REAL)") -parser.add_option("--EP2-imag", metavar = "name", default = "CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG", help = "Set the name of the channel containing the imaginary part of the factors used to compute A(f_ctrl). (Default = CAL-CS_TDEP_REF_CLGRATIO_CTRL_IMAG)") -parser.add_option("--EP6-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL", help = "Set the name of the channel containing the real part of C_res at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_REAL") -parser.add_option("--EP6-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG", help = "Set the name of the channel containing the imaginary part of C_res at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_C_NOCAVPOLE_IMAG") -parser.add_option("--EP7-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL", help = "Set the name of the channel containing the real part of D at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_D_REAL") -parser.add_option("--EP7-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG", help = "Set the name of the channel containing the imaginary part of D at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_D_IMAG") -parser.add_option("--EP8-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_REAL") -parser.add_option("--EP8-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_TST_IMAG") -parser.add_option("--EP9-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_REAL") -parser.add_option("--EP9-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the PCal line used for the \kappa_c and f_cc calculation. (Default = CAL-CS_TDEP_PCALY_LINE2_REF_A_USUM_IMAG") -parser.add_option("--EP1-real", metavar = "name", default = "CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL", help = "Set the name of the channel containing the real part of the \kappa_tst reference factors. (Default = CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_REAL)") -parser.add_option("--EP1-imag", metavar = "name", default = "CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG", help = "Set the name of the channel containing the imaginary part of the \kappa_tst reference factors. (Default = CAL-CS_TDEP_REF_INVA_CLGRATIO_TST_IMAG)") parser.add_option("--EP10-real", metavar = "name", default = "CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_REAL", help = "Set the name of the channel containing the real part of A_tst at the ESD line used for removal of the ESD line. (Default = CAL-CS_TDEP_ESD_LINE1_REF_A_TST_REAL") parser.add_option("--EP10-imag", metavar = "name", default = "CAL-CS_TDEP_ESD_LINE1_REF_A_TST_NOLOCK_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the ESD line used for removal of the ESD line. (Default = CAL-CS_TDEP_ESD_LINE1_REF_A_TST_IMAG") -parser.add_option("--EP11-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL", help = "Set the name of the channel containing the real part of C_res at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_REAL") -parser.add_option("--EP11-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG", help = "Set the name of the channel containing the imaginary part of C_res at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_C_NOCAVPOLE_IMAG") -parser.add_option("--EP12-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL", help = "Set the name of the channel containing the real part of D at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_D_REAL") -parser.add_option("--EP12-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG", help = "Set the name of the channel containing the imaginary part of D at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_D_IMAG") -parser.add_option("--EP13-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL", help = "Set the name of the channel containing the real part of A_tst at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_REAL") -parser.add_option("--EP13-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG", help = "Set the name of the channel containing the imaginary part of A_tst at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_TST_IMAG") -parser.add_option("--EP14-real", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL", help = "Set the name of the channel containing the real part of A_pu at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_REAL") -parser.add_option("--EP14-imag", metavar = "name", default = "CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG", help = "Set the name of the channel containing the imaginary part of A_pu at the PCal line used for the src_Q and f_s calculation. (Default = CAL-CS_TDEP_PCALY_LINE4_REF_A_USUM_IMAG") # Parse options @@ -190,17 +166,20 @@ data_sources = set(("frames", "lvshm")) if options.data_source not in data_sources: raise ValueError("--data-source must be one of %s" % ",".join(data_sources)) -if options.data_source == "frames" and options.frame_cache is None: - raise ValueError("--frame-cache must be specified when using --data-source=frames") +if options.data_source == "frames" and (options.hoft_frame_cache or options.raw_frame_cache) is None: + raise ValueError("--raw-frame-cache and --hoft-frame-cache must be specified when using --data-source=frames") if options.ifo is None: raise ValueError("must specify --ifo") -if options.frame_segments_file is not None and options.data_source != "frames": - raise ValueError("can only give --frame-segments-file if --data-source=frames") +if (options.raw_frame_segments_file or options.hoft_frame_segments_file) is not None and options.data_source != "frames": + raise ValueError("can only give --raw-frame-segments-file or --hoft-frame-segments-file if --data-source=frames") -if options.frame_segments_name is not None and options.frame_segments_file is None: - raise ValueError("can only specify --frame-segments-name if --frame-segments-file is given") +if options.raw_frame_segments_name is not None and options.raw_frame_segments_file is None: + raise ValueError("can only specify --raw-frame-segments-name if --raw-frame-segments-file is given") + +if options.hoft_frame_segments_name is not None and options.hoft_frame_segments_file is None: + raise ValueError("can only specify --hoft-frame-segments-name if --hoft-frame-segments-file is given") if options.data_source == "frames" and (options.gps_start_time is None or options.gps_end_time is None): raise ValueError("must specify --gps-start-time and --gps-end-time when --data-source=frames") @@ -235,26 +214,29 @@ elif options.gps_end_time is not None: instrument = options.ifo # Make segment list if a frame segments file is provided, other set frame_segments to None -if options.frame_segments_file is not None: +if options.raw_frame_segments_file is not None: + # Frame segments from a user defined file + raw_frame_segments = ligolw_segments.segmenttable_get_by_name(utils.load_filename(options.raw_frame_segments_file, contenthandler = datasource.ContentHandler), options.raw_frame_segments_name).coalesce() + if seg is not None: + # clip frame segments to seek segment if it exists (not required, just saves some meory and I/O overhead) + raw_frame_segments = segments.segmentlistdict((instrument, seglist & segments.segmentlist([seg])) for instrument, seglist in raw_frame_segments.items()) +else: + raw_frame_segments = None + +if options.hoft_frame_segments_file is not None: # Frame segments from a user defined file - frame_segments = ligolw_segments.segmenttable_get_by_name(utils.load_filename(options.frame_segments_file, contenthandler = datasource.ContentHandler), options.frame_segments_name).coalesce() + hoft_frame_segments = ligolw_segments.segmenttable_get_by_name(utils.load_filename(options.hoft_frame_segments_file, contenthandler = datasource.ContentHandler), options.hoft_frame_segments_name).coalesce() if seg is not None: # clip frame segments to seek segment if it exists (not required, just saves some meory and I/O overhead) - frame_segments = segments.segmentlistdict((instrument, seglist & segments.segmentlist([seg])) for instrument, seglist in frame_segments.items()) + hoft_frame_segments = segments.segmentlistdict((instrument, seglist & segments.segmentlist([seg])) for instrument, seglist in hoft_frame_segments.items()) else: - frame_segments = None + hoft_frame_segments = None # Set up short-cut names for each of the sample rates used throughout the pipeline and establish caps string shortcuts hoftsr = options.hoft_sample_rate # Sample rate for h(t) -calibstatesr = options.calib_state_sample_rate # Sample rate for the CALIB_STATE_VECTOR -odcsr = options.odc_sample_rate # Sample rate of the ODC channel that is read in -ctrlsr = options.control_sample_rate # Sample rate of the control channel (such as DARM_CTRL or DELTAL_CTRL) -cohsr = options.coh_sample_rate # Sample rate for the coherence uncertainty channels +kappasr = options.kappa_sample_rate # sample rate for kappa_tst hoft_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr -ctrl_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ctrlsr -calibstate_caps = "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % calibstatesr -odc_caps = "audio/x-raw, format=U32LE, rate=%d, channel-mask=(bitmask)0x0" % odcsr -coh_caps = "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % cohsr +kappa_caps = "audio/x-raw, format=F32LE, rate=%d, channel-mask=(bitmask)0x0" % kappasr # Set up string for the channels suffix and prefix as provided by the user if options.chan_suffix is not None: @@ -397,7 +379,10 @@ if options.low_latency: queue_factor = 1 else: queue_factor = 0 -long_queue = queue_factor * max(float(len(reschainfilt) - 1) / hoftsr, float(len(tstfilt) - 1) / tstchainsr, float(len(pumuimfilt) - 1) / pumuimchainsr) + +#FIXME: Replace long queue logic with something real +#long_queue = queue_factor * max(float(len(reschainfilt) - 1) / hoftsr, float(len(tstfilt) - 1) / tstchainsr, float(len(pumuimfilt) - 1) / pumuimchainsr) +long_queue = queue_factor * 16384 short_queue = -1.0 * queue_factor # @@ -406,41 +391,41 @@ short_queue = -1.0 * queue_factor # a dictionary that will be populated with pipeline branch names based on the channel list. # +hoft_head_dict = {} head_dict = {} +hoft_channel_list = [] channel_list = [] +hoft_headkeys = [] headkeys = [] -# If we are computing the CALIB_STATE_VECTOR, we need the ODC state vector +# We will need to read in the h(t) channel +hoft_channel_list.append((instrument, options.hoft_channel_name)) +hoft_headkeys.append("hoft") + +""" # FIXME: Do we need to read in the statevector? Maybe, because maybe the filter settle time is different for cleaned data... But should we really be modifying the CALIB_STATE_VECTOR? probably not... if not options.no_dq_vector: - channel_list.append((instrument, options.dq_channel_name)) - headkeys.append("odcstatevector") - -# Read in the h(t) channel + hoft_channel_list.append((instrument, options.dq_channel_name)) + hoft_headkeys.append("calibstatevector") +""" -# If we are computing the factors in the pipeline, we need the reference model EPICS records -if not options.factors_from_filters_file: - # EP10 is needed to remove the ESD line - if options.remove_callines and remove_esd_act_line and not options.factors_from_filters_file: +if options.remove_callines: + if remove_esd_act_line and not options.factors_from_filters_file: + # EP10 is needed to remove the ESD line channel_list.extend(((instrument, options.EP10_real), (instrument, options.EP10_imag))) headkeys.extend(("EP10_real", "EP10_imag")) -# We need to make sure we have DARM_ERR and the PCAL channel for computing \kappas -if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not options.no_kappapu or not options.no_fs or not options.no_srcQ: + # read in the pcal channel channel_list.append((instrument, options.pcal_channel_name)) headkeys.append("pcal") - if options.partial_calibration: - channel_list.append((instrument, options.darm_err_channel_name)) - headkeys.append("darm_err") - -# For full calibration we need DARM_ERR and DARM_CTRL as our input channels -if options.full_calibration: - channel_list.extend(((instrument, options.darm_err_channel_name), (instrument, options.darm_ctrl_channel_name))) - headkeys.extend(("res", "ctrl")) -# For partial calibration we need DELTAL_TST, DELTAL_PUM, DELTAL_UIM, and DELTAL_RES -elif options.partial_calibration: - channel_list.extend(((instrument, options.deltal_res_channel_name), (instrument, options.deltal_tst_channel_name), (instrument, options.deltal_pum_channel_name), (instrument, options.deltal_uim_channel_name))) - headkeys.extend(("res", "tst", "pum", "uim")) + + # read in the tst excitation channel + channel_list.append((instrument, options.tst_exc_channel_name)) + headkeys.append("tstexc") + + # read in kappa_tst + hoft_channel_list.append((instrument, options.kappa_tst_channel_name)) + hoft_headkeys.append("kappa_tst") # If we are removing additional noise from the spectrum (beam jitter, angular control, 60 Hz lines, etc.), we need more channels if options.remove_powerlines: @@ -464,7 +449,7 @@ if options.remove_length_control: ####################################### Main Pipeline ############################################## #################################################################################################### -pipeline = Gst.Pipeline(name="gstlal_compute_strain") +pipeline = Gst.Pipeline(name="gstlal_clean_strain") mainloop = GObject.MainLoop() handler = simplehandler.Handler(mainloop, pipeline) @@ -481,389 +466,73 @@ if not options.verbose: # if options.data_source == "lvshm": # Data is to be read from shared memory; "low-latency" mode - src = pipeparts.mklvshmsrc(pipeline, shm_name = options.shared_memory_partition, assumed_duration = 1) + hoft_src = pipeparts.mklvshmsrc(pipeline, shm_name = options.hoft_shared_memory_partition, assumed_duration = 4) + raw_src = pipeparts.mklvshmsrc(pipeline, shm_name = options.raw_shared_memory_partition, assumed_duration = 1) elif options.data_source == "frames": # Data is to be read from frame files; "offline" mode - src = pipeparts.mklalcachesrc(pipeline, location = options.frame_cache, cache_dsc_regex = instrument) + hoft_src = pipeparts.mklalcachesrc(pipeline, location = options.hoft_frame_cache, cache_dsc_regex = instrument) + raw_src = pipeparts.mklalcachesrc(pipeline, location = options.raw_frame_cache, cache_dsc_regex = instrument) # # Hook up the relevant channels to the demuxer # if options.data_source == "lvshm": - demux = pipeparts.mkframecppchanneldemux(pipeline, src, do_file_checksum = options.do_file_checksum, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list)) - + hoft_demux = pipeparts.mkframecppchanneldemux(pipeline, hoft_src, do_file_checksum = options.do_file_checksum, skip_bad_files = True, hoft_channel_list = map("%s:%s".__mod__, hoft_channel_list)) + raw_demux = pipeparts.mkframecppchanneldemux(pipeline, raw_src, do_file_checksum = options.do_file_checksum, skip_bad_files = True, channel_list = map("%s:%s".__mod__, channel_list)) elif options.data_source == "frames": - demux = pipeparts.mkframecppchanneldemux(pipeline, src, do_file_checksum = options.do_file_checksum, skip_bad_files = False, channel_list = map("%s:%s".__mod__, channel_list)) + hoft_demux = pipeparts.mkframecppchanneldemux(pipeline, hoft_src, do_file_checksum = options.do_file_checksum, skip_bad_files = False, hoft_channel_list = map("%s:%s".__mod__, hoft_channel_list)) + raw_demux = pipeparts.mkframecppchanneldemux(pipeline, raw_src, do_file_checksum = options.do_file_checksum, skip_bad_files = False, channel_list = map("%s:%s".__mod__, channel_list)) # Write the pipeline graph after pads have been hooked up to the demuxer if options.write_pipeline is not None: - demux.connect("no-more-pads", write_graph) + raw_demux.connect("no-more-pads", write_graph) + hoft_demux.connect("no-more-pads", write_graph) # Get everything hooked up and fill in discontinuities +for key, chan in zip(hoft_headkeys, hoft_channel_list): + hoft_head_dict[key] = calibration_parts.hook_up(pipeline, hoft_demux, chan[1], instrument, options.buffer_length) for key, chan in zip(headkeys, channel_list): - head_dict[key] = calibration_parts.hook_up(pipeline, demux, chan[1], instrument, options.buffer_length) + head_dict[key] = calibration_parts.hook_up(pipeline, raw_demux, chan[1], instrument, options.buffer_length) # When reading from disk, clip the incoming data stream(s) to segment list if one is provided -if options.data_source == "frames" and frame_segments is not None: +if options.data_source == "frames" and raw_frame_segments is not None: for key in headkeys: currenthead = head_dict[key] - head_dict[key] = calibration_parts.mkgate(pipeline, currenthead, pipeparts.mksegmentsrc(pipeline, frame_segments[instrument]), 1, long_queue, long_queue) + head_dict[key] = calibration_parts.mkgate(pipeline, currenthead, pipeparts.mksegmentsrc(pipeline, raw_frame_segments[instrument]), 1, long_queue, long_queue) + +if options.data_source == "frames" and hoft_frame_segments is not None: + for key in hoft_headkeys: + currenthead = hoft_head_dict[key] + hoft_head_dict[key] = calibration_parts.mkgate(pipeline, currenthead, pipeparts.mksegmentsrc(pipeline, hoft_frame_segments[instrument]), 1, long_queue, long_queue) # -# TIME-VARYING FACTORS COMPUTATIONS +# REMOVE CALIBRATION LINES # -for key in headkeys: - if key.startswith("EP"): - head_dict[key] = calibration_parts.caps_and_progress(pipeline, head_dict[key], ref_factors_caps, key) - head_dict[key] = calibration_parts.mkresample(pipeline, head_dict[key], 0, False, compute_calib_factors_caps) - -if not options.no_coherence: - if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: - pcaly_line1_coh = calibration_parts.caps_and_progress(pipeline, head_dict["pcaly_line1_coh"], coh_caps, "pcaly_line1_coh") - pcaly_line1_coh = calibration_parts.mkresample(pipeline, pcaly_line1_coh, 0, False, compute_calib_factors_caps) - sus_coh = calibration_parts.caps_and_progress(pipeline, head_dict["sus_coh"], coh_caps, "sus_coh") - sus_coh = calibration_parts.mkresample(pipeline, sus_coh, 0, False, compute_calib_factors_caps) - darm_coh = calibration_parts.caps_and_progress(pipeline, head_dict["darm_coh"], coh_caps, "darm_coh") - darm_coh = calibration_parts.mkresample(pipeline, darm_coh, 0, False, compute_calib_factors_caps) - if not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ: - pcaly_line2_coh = calibration_parts.caps_and_progress(pipeline, head_dict["pcaly_line2_coh"], coh_caps, "pcaly_line2_coh") - pcaly_line2_coh = calibration_parts.mkresample(pipeline, pcaly_line2_coh, 0, False, compute_calib_factors_caps) - -if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_fs or not options.no_srcQ or (options.remove_callines and remove_esd_act_line): +if options.remove_cal_lines: + EP10_real = calibration_parts.caps_and_progress(pipeline, head_dict["EP10_real"], ref_factors_caps, "EP10_real") + EP10_real = calibration_parts.mkresample(pipeline, EP10_real, 0, False, compute_calib_factors_caps) + EP10_imag = calibration_parts.caps_and_progress(pipeline, head_dict["EP10_imag"], ref_factors_caps, "EP10_imag") + EP10_imag = calibration_parts.mkresample(pipeline, EP10_imag, 0, False, compute_calib_factors_caps) + tstexccaps = "audio/x-raw, format=F64LE, rate=%d" % options.tst_exc_sample_rate tstexc = calibration_parts.caps_and_progress(pipeline, head_dict["tstexc"], tstexccaps, "tstexc") -if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_fs or not options.no_srcQ: - exc = calibration_parts.caps_and_progress(pipeline, head_dict["exc"], hoft_caps, "exc") - -# If we use the coherence channels multiple times in the pipeline, we need to tee the channels -if not options.no_coherence: - lowfreq_coh_use = int(not options.no_kappatst) + int(not options.no_kappapu) + int(not options.no_kappac) + int(not options.no_fcc) + int(not options.no_srcQ) + int(not options.no_fs) + int(not options.no_dq_vector and (not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs)) - highfreq_coh_use = int(not options.no_kappac) + int(not options.no_fcc) + int(not options.no_srcQ) + int(not options.no_fs) + int(not options.no_dq_vector and (not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs)) - if lowfreq_coh_use >= 2: - pcaly_line1_coh = pipeparts.mktee(pipeline, pcaly_line1_coh) - sus_coh = pipeparts.mktee(pipeline, sus_coh) - darm_coh = pipeparts.mktee(pipeline, darm_coh) - if highfreq_coh_use >= 2: - pcaly_line2_coh = pipeparts.mktee(pipeline, pcaly_line2_coh) - -# Set up computations for \kappa_a, \kappa_tst,\kappa_c, \kappa_pu, f_cc, if applicable -if not options.no_kappac or not options.no_fcc or not options.no_kappatst or not options.no_kappapu or not options.no_srcQ or not options.no_fs: - # pcal excitation channel, which will be demodulated pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], hoft_caps, "pcal") pcaltee = pipeparts.mktee(pipeline, pcal) - # DARM_ERR channel, which will have followed different paths if we're doing full vs. partial calibration - if options.full_calibration: - darm_err = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "darm_err") - else: - darm_err = calibration_parts.caps_and_progress(pipeline, head_dict["darm_err"], hoft_caps, "darm_err") - derrtee = pipeparts.mktee(pipeline, darm_err) - # demodulate the PCAL channel and apply the PCAL correction factor at the DARM actuation line frequency pcal_at_darm_act_freq = calibration_parts.demodulate(pipeline, pcaltee, darm_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_darm_act_freq_real, pcal_corr_at_darm_act_freq_imag) - if not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs or options.remove_callines: - pcal_at_darm_act_freq = pipeparts.mktee(pipeline, pcal_at_darm_act_freq) - - # demodulate DARM_ERR at the DARM actuation line frequency - derr_at_darm_act_freq = calibration_parts.demodulate(pipeline, derrtee, darm_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.dewhitening: - # dewhiten DARM_ERR at the DARM actuation line frequency - derr_at_darm_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_darm_act_freq, derr_dewhiten_at_darm_act_freq_real, derr_dewhiten_at_darm_act_freq_imag) - if not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: - derr_at_darm_act_freq = pipeparts.mktee(pipeline, derr_at_darm_act_freq) # demodulate the TST excitation channel at the ESD actuation line frequency tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.remove_callines and remove_esd_act_line: - tstexc_at_esd_act_freq = pipeparts.mktee(pipeline, tstexc_at_esd_act_freq) - - # demodulate DARM_ERR at the ESD actuation line frequency - derr_at_esd_act_freq = calibration_parts.demodulate(pipeline, derrtee, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.dewhitening: - # dewhiten DARM_ERR at the ESD actuation line frequency - derr_at_esd_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_esd_act_freq, derr_dewhiten_at_esd_act_freq_real, derr_dewhiten_at_esd_act_freq_imag) - - # compute kappa_tst, either using reference factors from the filters file or reading them from EPICS channels - if not options.factors_from_filters_file: - EP1 = calibration_parts.merge_into_complex(pipeline, head_dict["EP1_real"], head_dict["EP1_imag"], long_queue, short_queue) - ktst = calibration_parts.compute_kappatst(pipeline, derr_at_esd_act_freq, tstexc_at_esd_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP1, long_queue, short_queue) - elif options.factors_from_filters_file: - ktst = calibration_parts.compute_kappatst_from_filters_file(pipeline, derr_at_esd_act_freq, tstexc_at_esd_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP1_real, EP1_imag, long_queue, short_queue) - - if not options.no_kappatst or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: - ktst = pipeparts.mktee(pipeline, ktst) - if not options.no_kappatst: - smooth_ktst_nogate = pipeparts.mkgeneric(pipeline, ktst, "lal_smoothkappas", default_kappa_re = options.expected_kappatst_real, default_kappa_im = options.expected_kappatst_imag, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - smooth_ktstR_nogate, smooth_ktstI_nogate = calibration_parts.split_into_real(pipeline, smooth_ktst_nogate) - - if not options.no_coherence: - # Gate kappa_tst with the coherence of the PCALY_line1 line - ktst_gated = calibration_parts.mkgate(pipeline, ktst, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - # Gate kappa_tst with the coherence of the suspension line - ktst_gated = calibration_parts.mkgate(pipeline, ktst_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - # Gate kappa_tst with the coherence of the DARM line - ktst_gated = calibration_parts.mkgate(pipeline, ktst_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - ktst_gated = pipeparts.mktee(pipeline, ktst_gated) - smooth_ktstRdq, smooth_ktstIdq = calibration_parts.track_bad_complex_kappas(pipeline, ktst_gated, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - # Smooth kappa_tst - smooth_ktstR, smooth_ktstI = calibration_parts.smooth_complex_kappas(pipeline, ktst_gated, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - else: - if not options.no_dq_vector: - smooth_ktstRdq, smooth_ktstIdq = calibration_parts.track_bad_complex_kappas_no_coherence(pipeline, ktst, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - # Smooth kappa_tst - smooth_ktstR, smooth_ktstI = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, ktst, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, options.expected_kappatst_real, options.expected_kappatst_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - smooth_ktstRtee = pipeparts.mktee(pipeline, smooth_ktstR) - smooth_ktstItee = pipeparts.mktee(pipeline, smooth_ktstI) - -# If we're also computing \kappa_c, f_cc, or \kappa_pu, keep going -if not options.no_kappac or not options.no_fcc or not options.no_kappapu or not options.no_srcQ or not options.no_fs: - # demodulate excitation channel at PU actuation line frequency - exc_at_pu_act_freq = calibration_parts.demodulate(pipeline, exc, pu_act_esd_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - - # demodulate DARM_ERR at PU actuation line frequency - derr_at_pu_act_freq = calibration_parts.demodulate(pipeline, derrtee, pu_act_esd_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.dewhitening: - # dewhiten DARM_ERR at the PU actuation line frequency - derr_at_pu_act_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_pu_act_freq, derr_dewhiten_at_pu_act_freq_real, derr_dewhiten_at_pu_act_freq_imag) - - # compute the factor Afctrl that will be used in the computation of kappa_pu and kappa_a, either using reference factors from the filters file or reading them from EPICS channels - if not options.factors_from_filters_file: - EP2 = calibration_parts.merge_into_complex(pipeline, head_dict["EP2_real"], head_dict["EP2_imag"], long_queue, short_queue) - EP3 = calibration_parts.merge_into_complex(pipeline, head_dict["EP3_real"], head_dict["EP3_imag"], long_queue, short_queue) - EP4 = calibration_parts.merge_into_complex(pipeline, head_dict["EP4_real"], head_dict["EP4_imag"], long_queue, short_queue) - afctrl = calibration_parts.compute_afctrl(pipeline, derr_at_pu_act_freq, exc_at_pu_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP2, long_queue, short_queue) - elif options.factors_from_filters_file: - afctrl = calibration_parts.compute_afctrl_from_filters_file(pipeline, derr_at_pu_act_freq, exc_at_pu_act_freq, pcal_at_darm_act_freq, derr_at_darm_act_freq, EP2_real, EP2_imag, long_queue, short_queue) - - # \kappa_pu calcuation, which needs to happen for any of the other kappas to be computed - if not options.factors_from_filters_file: - kpu = calibration_parts.compute_kappapu(pipeline, EP3, afctrl, ktst, EP4, long_queue, short_queue) - elif options.factors_from_filters_file: - kpu = calibration_parts.compute_kappapu_from_filters_file(pipeline, EP3_real, EP3_imag, afctrl, ktst, EP4_real, EP4_imag, long_queue, short_queue) - - if not options.no_kappapu or not options.no_srcQ or not options.no_fs: - kpu = pipeparts.mktee(pipeline, kpu) - if not options.no_kappapu: - smooth_kpu_nogate = pipeparts.mkgeneric(pipeline, kpu, "lal_smoothkappas", default_kappa_re = options.expected_kappapu_real, default_kappa_im = options.expected_kappapu_imag, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - smooth_kpuR_nogate, smooth_kpuI_nogate = calibration_parts.split_into_real(pipeline, smooth_kpu_nogate) - - if not options.no_coherence: - # Gate kappa_pu with the coherence of the DARM line - kpu_gated = calibration_parts.mkgate(pipeline, kpu, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - # Gate kappa_pu with the coherence of the PCALY_line1 line - kpu_gated = calibration_parts.mkgate(pipeline, kpu_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - # Gate kappa_pu with the coherence of the suspension coherence - kpu_gated = calibration_parts.mkgate(pipeline, kpu_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - kpu_gated = pipeparts.mktee(pipeline, kpu_gated) - smooth_kpuRdq, smooth_kpuIdq = calibration_parts.track_bad_complex_kappas(pipeline, kpu_gated, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - # Smooth kappa_pu - smooth_kpuR, smooth_kpuI = calibration_parts.smooth_complex_kappas(pipeline, kpu_gated, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - else: - if not options.no_dq_vector: - smooth_kpuRdq, smooth_kpuIdq = calibration_parts.track_bad_complex_kappas_no_coherence(pipeline, kpu, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - # Smooth kappa_pu - smooth_kpuR, smooth_kpuI = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, kpu, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, options.expected_kappapu_real, options.expected_kappapu_imag, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_kpuRtee = pipeparts.mktee(pipeline, smooth_kpuR) - smooth_kpuItee = pipeparts.mktee(pipeline, smooth_kpuI) - - # Finally, compute \kappa_c and f_cc - if not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: - # demodulate PCAL channel and apply the PCAL correction factor at optical gain and f_cc line frequency - pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) - if options.remove_callines: - pcal_at_opt_gain_freq = pipeparts.mktee(pipeline, pcal_at_opt_gain_freq) - - # demodulate DARM_ERR at optical gain and f_cc line frequency - derr_at_opt_gain_freq = calibration_parts.demodulate(pipeline, derrtee, opt_gain_fcc_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - if options.dewhitening: - # dewhiten DARM_ERR at optical gain and f_cc line frequency - derr_at_opt_gain_freq = calibration_parts.complex_audioamplify(pipeline, derr_at_opt_gain_freq, derr_dewhiten_at_opt_gain_fcc_freq_real, derr_dewhiten_at_opt_gain_fcc_freq_imag) - - # Compute the factor S which will be used for the kappa_c and f_cc calculations - if not options.factors_from_filters_file: - EP6 = calibration_parts.merge_into_complex(pipeline, head_dict["EP6_real"], head_dict["EP6_imag"], long_queue, short_queue) - EP7 = calibration_parts.merge_into_complex(pipeline, head_dict["EP7_real"], head_dict["EP7_imag"], long_queue, short_queue) - EP8 = calibration_parts.merge_into_complex(pipeline, head_dict["EP8_real"], head_dict["EP8_imag"], long_queue, short_queue) - EP9 = calibration_parts.merge_into_complex(pipeline, head_dict["EP9_real"], head_dict["EP9_imag"], long_queue, short_queue) - S = calibration_parts.compute_S(pipeline, EP6, pcal_at_opt_gain_freq, derr_at_opt_gain_freq, EP7, ktst, EP8, kpu, EP9, long_queue, short_queue) - elif options.factors_from_filters_file: - S = calibration_parts.compute_S_from_filters_file(pipeline, EP6_real, EP6_imag, pcal_at_opt_gain_freq, derr_at_opt_gain_freq, EP7_real, EP7_imag, ktst, EP8_real, EP8_imag, kpu, EP9_real, EP9_imag, long_queue, short_queue) - - S = pipeparts.mktee(pipeline, S) - - SR, SI = calibration_parts.split_into_real(pipeline, S) - - if not options.no_kappac and not options.no_fcc: - SR = pipeparts.mktee(pipeline, SR) - SI = pipeparts.mktee(pipeline, SI) - - # compute kappa_c - if not options.no_kappac or not options.no_srcQ or not options.no_fs: - kc = calibration_parts.compute_kappac(pipeline, SR, SI, long_queue, short_queue) - if not options.no_kappac: - kc = pipeparts.mktee(pipeline, kc) - smooth_kc_nogate = pipeparts.mkgeneric(pipeline, kc, "lal_smoothkappas", default_kappa_re = options.expected_kappac, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - - if not options.no_coherence: - # Gate kappa_c with all four of the calibration lines - kc_gated = calibration_parts.mkgate(pipeline, kc, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - kc_gated = calibration_parts.mkgate(pipeline, kc_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - kc_gated = calibration_parts.mkgate(pipeline, kc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - kc_gated = calibration_parts.mkgate(pipeline, kc_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - kc_gated = pipeparts.mktee(pipeline, kc_gated) - smooth_kcdq = calibration_parts.track_bad_kappas(pipeline, kc_gated, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_kc = calibration_parts.smooth_kappas(pipeline, kc_gated, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - else: - smooth_kc = calibration_parts.smooth_kappas_no_coherence(pipeline, kc, options.kappac_ok_var, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - if not options.no_dq_vector: - smooth_kcdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, kc, options.kappac_ok_var, options.expected_kappac, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_kctee = pipeparts.mktee(pipeline, smooth_kc) - - # compute f_cc - if not options.no_fcc or not options.nosrcQ or not options.no_fs: - fcc = calibration_parts.compute_fcc(pipeline, SR, SI, opt_gain_fcc_line_freq, long_queue, short_queue) - if not options.no_fcc: - fcc = pipeparts.mktee(pipeline, fcc) - smooth_fcc_nogate = pipeparts.mkgeneric(pipeline, fcc, "lal_smoothkappas", default_kappa_re = fcc_default, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - - if not options.no_coherence: - # Gate f_cc with all four of the calibration lines - fcc_gated = calibration_parts.mkgate(pipeline, fcc, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fcc_gated = calibration_parts.mkgate(pipeline, fcc_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - fcc_gated = pipeparts.mktee(pipeline, fcc_gated) - smooth_fccdq = calibration_parts.track_bad_kappas(pipeline, fcc_gated, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_fcc = calibration_parts.smooth_kappas(pipeline, fcc_gated, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - else: - smooth_fcc = calibration_parts.smooth_kappas_no_coherence(pipeline, fcc, options.fcc_ok_var, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - if not options.no_dq_vector: - smooth_fccdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, fcc, options.fcc_ok_var, fcc_default, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_fcctee = pipeparts.mktee(pipeline, smooth_fcc) - - if not options.no_fs or not options.no_srcQ: - # demodulate PCAL channel and apply the PCAL correction factor at SRC detuning line frequency - pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) - if options.remove_callines and remove_src_pcal_line: - pcal_at_src_freq = pipeparts.mktee(pipeline, pcal_at_src_freq) - - # demodulate DARM_ERR at SRC detuning line frequency - derr_at_src_freq = calibration_parts.demodulate(pipeline, derrtee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples) - - # Compute the factor Xi which will be used for the f_s and src_Q calculations - if not options.factors_from_filters_file: - EP11 = calibration_parts.merge_into_complex(pipeline, head_dict["EP11_real"], head_dict["EP11_imag"], long_queue, short_queue) - EP12 = calibration_parts.merge_into_complex(pipeline, head_dict["EP12_real"], head_dict["EP12_imag"], long_queue, short_queue) - EP13 = calibration_parts.merge_into_complex(pipeline, head_dict["EP13_real"], head_dict["EP13_imag"], long_queue, short_queue) - EP14 = calibration_parts.merge_into_complex(pipeline, head_dict["EP14_real"], head_dict["EP14_imag"], long_queue, short_queue) - Xi = calibration_parts.compute_Xi(pipeline, pcal_at_src_freq, derr_at_src_freq, src_pcal_line_freq, EP11, EP12, EP13, EP14, ktst, kpu, kc, fcc, long_queue, short_queue) - elif options.factors_from_filters_file: - Xi = calibration_parts.compute_Xi_from_filters_file(pipeline, pcal_at_src_freq, derr_at_src_freq, src_pcal_line_freq, EP11_real, EP11_imag, EP12_real, EP12_imag, EP13_real, EP13_imag, EP14_real, EP14_imag, ktst, kpu, kc, fcc, long_queue, short_queue) - - XiR, XiI = calibration_parts.split_into_real(pipeline, Xi) - - if options.no_srcQ: - # the imaginary part is only used to compute Q - pipeparts.mkfakesink(pipeline, XiI) - - sqrtXiR = pipeparts.mkpow(pipeline, XiR, exponent = 0.5) - if not options.no_fs and not options.no_srcQ: - sqrtXiR = pipeparts.mktee(pipeline, sqrtXiR) - - # compute f_s - if not options.no_fs: - fs = pipeparts.mkaudioamplify(pipeline, sqrtXiR, src_pcal_line_freq) - fs = pipeparts.mktee(pipeline, fs) - smooth_fs_nogate = pipeparts.mkgeneric(pipeline, fs, "lal_smoothkappas", default_kappa_re = options.expected_fs, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - - if not options.no_coherence: - # Gate f_s with all four of the calibration lines - fs_gated = calibration_parts.mkgate(pipeline, fs, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fs_gated = calibration_parts.mkgate(pipeline, fs_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fs_gated = calibration_parts.mkgate(pipeline, fs_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - fs_gated = calibration_parts.mkgate(pipeline, fs_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - fs_gated = pipeparts.mktee(pipeline, fs_gated) - smooth_fsdq = calibration_parts.track_bad_kappas(pipeline, fs_gated, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_fs = calibration_parts.smooth_kappas(pipeline, fs_gated, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - else: - smooth_fs = calibration_parts.smooth_kappas_no_coherence(pipeline, fs, options.fs_ok_var, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - if not options.no_dq_vector: - smooth_fsdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, fs, options.fs_ok_var, options.expected_fs, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - if not options.no_dq_vector: - smooth_fs = pipeparts.mktee(pipeline, smooth_fs) - - # compute SRC Q_inv - if not options.no_srcQ: - sqrtXiR_inv = pipeparts.mkpow(pipeline, sqrtXiR, exponent = -1.0) - sqrtXiR_inv = pipeparts.mkaudioamplify(pipeline, sqrtXiR_inv, -1.0) - srcQ_inv = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [sqrtXiR_inv, long_queue], [XiI, short_queue])) - srcQ_inv = pipeparts.mktee(pipeline, srcQ_inv) - smooth_srcQ_inv_nogate = pipeparts.mkgeneric(pipeline, srcQ_inv, "lal_smoothkappas", default_kappa_re = 1.0 / options.expected_srcQ, array_size = median_smoothing_samples, avg_array_size = factors_average_samples, default_to_median = options.kappas_default_to_median) - - if not options.no_coherence: - # Gate SRC_Q with all four of the calibration lines - srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv, pcaly_line2_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, darm_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, pcaly_line1_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - srcQ_inv_gated = calibration_parts.mkgate(pipeline, srcQ_inv_gated, sus_coh, options.coherence_uncertainty_threshold, short_queue, long_queue, attack_length = -integration_samples, invert_control = True) - - if not options.no_dq_vector: - srcQ_inv_gated = pipeparts.mktee(pipeline, srcQ_inv_gated) - smooth_srcQdq = calibration_parts.track_bad_kappas(pipeline, srcQ_inv_gated, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - smooth_srcQ_inv = calibration_parts.smooth_kappas(pipeline, srcQ_inv_gated, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - else: - smooth_srcQ_inv = calibration_parts.smooth_kappas_no_coherence(pipeline, srcQ_inv, options.srcQ_ok_var, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - if not options.no_dq_vector: - smooth_srcQdq = calibration_parts.track_bad_kappas_no_coherence(pipeline, srcQ_inv, options.srcQ_ok_var, 1.0 / options.expected_srcQ, median_smoothing_samples, factors_average_samples, options.kappas_default_to_median) - - if not options.no_dq_vector: - smooth_srcQ_inv = pipeparts.mktee(pipeline, smooth_srcQ_inv) - -# -# Calibration Line Removal -# + # demodulate PCAL channel and apply the PCAL correction factor at optical gain and f_cc line frequency + pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) -if options.remove_callines: - # if we didn't compute the kappas, we still need to get the pcal channel - if options.no_kappatst and options.no_kappapu and options.no_kappac and options.no_fcc and options.no_srcQ and options.no_fs: - pcal = calibration_parts.caps_and_progress(pipeline, head_dict["pcal"], hoft_caps, "pcal") - pcaltee = pipeparts.mktee(pipeline, pcal) - # Demodulate pcal at the ~30 Hz pcal line - pcal_at_darm_act_freq = calibration_parts.demodulate(pipeline, pcal_tee, darm_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_darm_act_freq_real, pcal_corr_at_darm_act_freq_imag) + # demodulate PCAL channel and apply the PCAL correction factor at SRC detuning line frequency + pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) # Reconstruct a calibrated (negative) pcal at only the ~30 Hz pcal line pcaly_line1 = calibration_parts.mkresample(pipeline, pcal_at_darm_act_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) @@ -871,9 +540,6 @@ if options.remove_callines: remove_pcaly_line1, trash = calibration_parts.split_into_real(pipeline, pcaly_line1) pipeparts.mkfakesink(pipeline, trash) - # Make sure we have demodulated pcal at the ~300 Hz pcal line - if options.no_kappac and options.no_fcc and options.no_srcQ and options.no_fs: - pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_opt_gain_fcc_freq_real, pcal_corr_at_opt_gain_fcc_freq_imag) # Reconstruct a calibrated (negative) pcal at only the ~300 Hz pcal line pcaly_line2 = calibration_parts.mkresample(pipeline, pcal_at_opt_gain_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) pcaly_line2 = pipeparts.mkgeneric(pipeline, pcaly_line2, "lal_demodulate", line_frequency = -1.0 * opt_gain_fcc_line_freq, prefactor_real = 2.0) @@ -884,19 +550,17 @@ if options.remove_callines: remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_pcaly_line1, long_queue], [remove_pcaly_line2, short_queue])) if remove_esd_act_line: - # Make sure we have demodulated the ESD excitation channel at the ~30 Hz ESD line - if options.no_kappac and options.no_fcc and options.no_kappatst and options.no_kappapu and options.no_srcQ and options.no_fs: - tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_calib_factors_complex_caps, integration_samples) if options.factors_from_filters_file: esd_act_line = calibration_parts.complex_audioamplify(pipeline, tstexc_at_esd_act_freq, EP10_real, EP10_imag) else: # EP10 was read from the frames - EP10 = calibration_parts.merge_into_complex(pipeline, head_dict["EP10_real"], head_dict["EP10_imag"], long_queue, short_queue) + EP10 = calibration_parts.merge_into_complex(pipeline, EP10_real, EP10_imag, long_queue, short_queue) esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [tstexc_at_esd_act_freq, long_queue], [EP10, short_queue])) # Reconstruct a calibrated (negative) ESD injection at the ~30 Hz ESD line - if options.apply_kappatst: + if options.kappatst_applied: # Multiply by the real part of kappa_tst - esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [esd_act_line, long_queue], [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, smooth_ktstRtee, matrix=[[1.0, 0.0]])), short_queue])) + kappa_tst = calibration_parts.caps_and_progress(pipeline, hoft_head_dict["kappa_tst"], kappa_caps, "kappa_tst") + esd_act_line = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [esd_act_line, long_queue], [kappa_tst, short_queue])) esd_act_line = calibration_parts.mkresample(pipeline, esd_act_line, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) esd_act_line_remove = pipeparts.mkgeneric(pipeline, esd_act_line, "lal_demodulate", line_frequency = -1.0 * esd_act_line_freq, prefactor_real = 2.0) esd_act_line_remove, trash = calibration_parts.split_into_real(pipeline, esd_act_line_remove) @@ -927,9 +591,6 @@ if options.remove_callines: remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line4, short_queue])) if remove_src_pcal_line: - # Make sure we have demodulated pcal at the ~8 Hz pcal line - if options.no_fs and options.no_srcQ: - pcal_at_src_freq = calibration_parts.demodulate(pipeline, pcaltee, src_pcal_line_freq, td, compute_calib_factors_complex_caps, integration_samples, pcal_corr_at_src_freq_real, pcal_corr_at_src_freq_imag) # Reconstruct a calibrated (negative) pcal at only the ~3kHz pcal line pcaly_line0 = calibration_parts.mkresample(pipeline, pcal_at_src_freq, 3, False, "audio/x-raw, format=Z128LE, rate=%d, channel-mask=(bitmask)0x0" % hoftsr) pcaly_line0 = pipeparts.mkgeneric(pipeline, pcaly_line0, "lal_demodulate", line_frequency = -1.0 * src_pcal_line_freq, prefactor_real = 2.0) @@ -938,6 +599,11 @@ if options.remove_callines: # Add into the total line removal stream remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [remove_from_strain, long_queue], [remove_pcaly_line0, short_queue])) + +# +# REMOVE POWER LINES +# + if options.remove_powerlines: powerlines = calibration_parts.caps_and_progress(pipeline, head_dict["powerlines"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % powerlinessr, "powerlines") powerlines = pipeparts.mkfirbank(pipeline, powerlines, latency = int(powerlinesdelay), fir_matrix = [powerlinesfilt[::-1]], time_domain = td) @@ -947,6 +613,10 @@ if options.remove_powerlines: else: remove_from_strain = powerlines +# +# REMOVE IMC JITTER +# + if options.remove_jitter_imc: imc_a_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcapitsr, "imc_a_pitch") imc_a_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["imc_a_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % imcayawsr, "imc_a_yaw") @@ -986,6 +656,10 @@ if options.remove_jitter_psl: else: remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [bullseye_width, long_queue], [bullseye_pitch, long_queue], [bullseye_yaw, long_queue])) +# +# REMOVE ANGULAR CONTROL +# + if options.remove_angular_control: asc_dhard_pitch = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_pitch"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdpitsr, "asc_dhard_pitch") asc_dhard_yaw = calibration_parts.caps_and_progress(pipeline, head_dict["asc_dhard_yaw"], "audio/x-raw, format=F64LE, rate=%d, channel-mask=(bitmask)0x0" % ascdyawsr, "asc_dhard_yaw") @@ -1025,504 +699,20 @@ if options.remove_length_control: else: remove_from_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [lsc_srcl, long_queue], [lsc_mich, long_queue], [lsc_prcl, long_queue])) -# -# CONTROL BRANCH -# - -# zero out filter settling samples -tst_filter_settle_time = 0.0 -tst_filter_latency = 0.0 -pumuim_filter_settle_time = 0.0 -pumuim_filter_latency = 0.0 - -# The reverse of the filters will be used in all filtering below due to the definition of the filtering procedure employed by lal_firbank -if options.partial_calibration: - # enforce caps on actuation channels and set up progress report if verbose is on - tst = calibration_parts.caps_and_progress(pipeline, head_dict["tst"], ctrl_caps, "tst") - tsttee = pipeparts.mktee(pipeline, tst) - pum = calibration_parts.caps_and_progress(pipeline, head_dict["pum"], ctrl_caps, "pum") - pumtee = pipeparts.mktee(pipeline, pum) - uim = calibration_parts.caps_and_progress(pipeline, head_dict["uim"], ctrl_caps, "uim") - uimtee = pipeparts.mktee(pipeline, uim) - - # add together the PUM and UIM actuation channels; this may change in the future... - pumuim = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [pumtee, long_queue], [uimtee, short_queue])) - - # if you need to, dewhiten the TST and PUM/UIM chains - if options.dewhitening: - pumuim = calibration_parts.mkstockresample(pipeline, pumuim, "audio/x-raw, format=F64LE, rate=%d" % pumuimdewhitensr) - pumuim = pipeparts.mkfirbank(pipeline, pumuim, latency = int(pumuimdewhitendelay), fir_matrix = [pumuimdewhiten[::-1]], time_domain = td) - pumuim_filter_settle_time += float(len(pumuimdewhiten)-pumuimdewhitendelay)/pumuimdewhitensr - pumuim_filter_latency += float(pumuimdewhitendelay)/pumuimdewhitensr - tst = calibration_parts.mkstockresample(pipeline, tsttee, "audio/x-raw, format=F64LE, rate=%d" % tstdewhitensr) - tst = pipeparts.mkfirbank(pipeline, tst, latency = int(tstdewhitendelay), fir_matrix = [tstdewhiten[::-1]], time_domain = td) - tst_filter_settle_time += float(len(tstdewhiten)-tstdewhitendelay)/tstdewhitensr - tst_filter_latency += float(tstdewhitendelay)/tstdewhitensr - else: - tst = tsttee - -if options.full_calibration: - # enforce caps on actuation channels and set up progress report, if verbose is on - ctrl = calibration_parts.caps_and_progress(pipeline, head_dict["ctrl"], hoft_caps, "ctrl") - darmctrltee = pipeparts.mktee(pipeline, ctrl) - - if options.dewhitening: - # dewhiten the DARM_CTRL channel - ctrl = calibration_parts.mkstockresample(pipeline, darmctrltee, "audio/x-raw, format=F64LE, rate=%d" % ctrldewhitensr) - ctrl = pipeparts.mkfirbank(pipeline, ctrl, latency = int(ctrldewhitendelay), fir_matrix = [ctrldewhiten[::-1]], time_domain = td) - tst_filter_settle_time += float(len(ctrldewhiten)-ctrldewhitendelay)/ctrldewhitensr - tst_filter_latency += float(ctrldewhitendelay)/ctrldewhitensr - pumuim_filter_settle_time += float(len(ctrldewhiten)-ctrldewhitendelay)/ctrldewhitensr - pumuim_filter_latency += float(ctrldewhitendelay)/ctrldewhitensr - # tee off DARM_CTRL to be filtered with PUM/UIM and TST filters separately - ctrltee = pipeparts.mktee(pipeline, ctrl) - else: - ctrltee = pipeparts.mktee(pipeline, darmctrltee) - tst = ctrltee - pumuim = ctrltee - -# resample what will become the TST actuation chain to the TST FIR filter sample rate -tst = calibration_parts.mkstockresample(pipeline, tst, "audio/x-raw, format=F64LE, rate=%d" % tstchainsr) -# filter TST chain with the TST acutaiton filter -tst = pipeparts.mkfirbank(pipeline, tst, latency = int(tstdelay), fir_matrix = [tstfilt[::-1]], time_domain = td) -if options.low_latency: - tst = calibration_parts.mkinsertgap(pipeline, tst, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) -tst_filter_settle_time += float(len(tstfilt)-tstdelay)/tstchainsr -tst_filter_latency += float(tstdelay)/tstchainsr -# resample the TST actuation chain to the full sample rate -if tstchainsr != pumuimchainsr or options.apply_kappatst or options.apply_kappapu: - tst = calibration_parts.mkstockresample(pipeline, tst, hoft_caps) - -# resample what will become the PUM/UIM actuation chain to the PUM/UIM FIR filter sample rate -pumuim = calibration_parts.mkstockresample(pipeline, pumuim, "audio/x-raw, format=F64LE, rate=%d" % pumuimchainsr) -# filter the PUM/UIM chain with the PUM/UIM actuation filter -pumuim = pipeparts.mkfirbank(pipeline, pumuim, latency = int(pumuimdelay), fir_matrix = [pumuimfilt[::-1]], time_domain = td) -if options.low_latency: - pumuim = calibration_parts.mkinsertgap(pipeline, pumuim, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) -pumuim_filter_settle_time += float(len(pumuimfilt)-pumuimdelay)/pumuimchainsr -pumuim_filter_latency += float(pumuimdelay)/pumuimchainsr -# resample the PUM/UIM actuation chain to the full sample rate -if tstchainsr != pumuimchainsr or options.apply_kappapu or options.apply_kappatst: - pumuim = calibration_parts.mkstockresample(pipeline, pumuim, hoft_caps) - -# apply kappa_tst -if options.apply_kappatst: - # Only apply the real part of \kappa_tst as a correction to A_tst - ktst_for_tst = calibration_parts.mkresample(pipeline, smooth_ktstRtee, 3, False, hoft_caps) - tst = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [ktst_for_tst, ktst_queue_length], [tst, tst_queue_length])) -# apply kappa_pu -if options.apply_kappapu: - # Only apply the real part of \kappa_pu as a correction to A_pu - kpu_for_pu = calibration_parts.mkresample(pipeline, smooth_kpuRtee, 3, False, hoft_caps) - pumuim = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [kpu_for_pu, kpu_queue_length], [pumuim, pumuim_queue_length])) - -# Add the TST and PUM/UIM chains together to form the full actuation chain -ctrl = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [tst, long_queue], [pumuim, short_queue])) -if tstchainsr != hoftsr or not options.apply_kappatst or not options.apply_kappapu: - ctrl = calibration_parts.mkstockresample(pipeline, ctrl, hoft_caps) - -# -# RESIDUAL BRANCH -# - -# zero out res filter settle time -res_filter_settle_time = 0.0 -res_filter_latency = 0.0 - -# The reverse of the filters will be used in all filtering below due to the definition of the filtering procedure employed by lal_firbank - -# enforce caps on the residual branch and hook up progress report if verbose is on -if options.full_calibration: - if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: - restee = derrtee - else: - res = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "res") - restee = pipeparts.mktee(pipeline, res) -if options.partial_calibration: - res = calibration_parts.caps_and_progress(pipeline, head_dict["res"], hoft_caps, "res") - restee = pipeparts.mktee(pipeline, res) - -# apply the residual chain filter -res = pipeparts.mkfirbank(pipeline, restee, latency = int(reschaindelay), fir_matrix = [reschainfilt[::-1]], time_domain = td) -if options.low_latency: - res = calibration_parts.mkinsertgap(pipeline, res, bad_data_intervals = [-1e35, 1e35], block_duration = 1000000000 * options.buffer_length, remove_gap = False) -res_filter_settle_time += float(len(reschainfilt)-reschaindelay)/hoftsr -res_filter_latency += float(reschaindelay)/hoftsr -if options.dewhitening: - res = pipeparts.mkfirbank(pipeline, res, latency = int(resdewhitendelay), fir_matrix = [resdewhiten[::-1]], time_domain = td) - res_filter_settle_time += float(len(resdewhiten)-resdewhitendelay)/hoftsr - res_filter_latency += float(resdewhitendelay)/hoftsr - -# Apply factors to actuation and sensing chains, if applicable -if options.apply_kappac: - kc_modify_res = calibration_parts.mkresample(pipeline, smooth_kctee, 3, False, hoft_caps) - res = calibration_parts.mkmultiplier(pipeline, calibration_parts.list_srcs(pipeline, [res, reskc_queue_length], [pipeparts.mkpow(pipeline, kc_modify_res, exponent = -1.0), kc_queue_length])) - -filter_settle_time = max(res_filter_settle_time, tst_filter_settle_time, pumuim_filter_settle_time) -filter_latency = max(res_filter_latency, tst_filter_latency, pumuim_filter_latency) # -# CONTROL + RESIDUAL = H(T) +# REMOVE EVERYTHING FROM THE STRAIN # -# Add control and residual chains and divide by L to make h(t) -strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [res, res_queue_length], [ctrl, ctrl_queue_length])) -# Remove the calibration lines, if we want -if options.remove_callines or options.remove_powerlines or options.remove_jitter_imc or options.remove_jitter_psl or options.remove_angular_control or options.remove_length_control: - remove_from_strain = pipeparts.mkaudioamplify(pipeline, remove_from_strain, -1.0) - strain = pipeparts.mktee(pipeline, strain) - cleaned_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [strain, short_queue], [remove_from_strain, long_queue])) -# Divide by L in a way that is compatitble with old and new filters files, since old filter files don't recored "arm length" -try: - strain = pipeparts.mkaudioamplify(pipeline, strain, 1.0/float(filters["arm_length"])) -except KeyError: - strain = pipeparts.mkaudioamplify(pipeline, strain, 1.0/3994.5) - -strain = pipeparts.mkprogressreport(pipeline, strain, "progress_hoft_%s" % instrument) +strain = calibration_parts.caps_and_progress(pipeline, hoft_head_dict["hoft"], hoft_caps, "hoft") +remove_from_strain = pipeparts.mkaudioamplify(pipeline, remove_from_strain, -1.0) +cleaned_strain = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [strain, short_queue], [remove_from_strain, long_queue])) -if options.remove_callines: - try: - cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/float(filters["arm_length"])) - except KeyError: - cleaned_strain = pipeparts.mkaudioamplify(pipeline, cleaned_strain, 1.0/3994.5) +cleaned_strain = pipeparts.mkprogressreport(pipeline, cleaned_strain, "progress_hoft_cleand_%s" % instrument) - cleaned_strain = pipeparts.mkprogressreport(pipeline, cleaned_strain, "progress_hoft_cleaned_%s" % instrument) # Put the units back to strain before writing to frames -straintagstr = "units=strain,channel-name=%sCALIB_STRAIN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) cleaned_straintagstr = "units=strain,channel-name=%sCALIB_STRAIN_CLEAN%s,instrument=%s" % (chan_prefix, chan_suffix, instrument) -if not options.no_dq_vector: - straintee = pipeparts.mktee(pipeline, strain) - strain = pipeparts.mktaginject(pipeline, straintee, straintagstr) -else: - strain = pipeparts.mktaginject(pipeline, strain, straintagstr) -if options.remove_callines: - cleaned_strain = pipeparts.mktaginject(pipeline, cleaned_strain, cleaned_straintagstr) - -# -# CALIB_STATE_VECTOR BRANCH -# - -#FIXME: Add more comments! - -if not options.no_dq_vector: - # FIXME: When the ODC is written as unsigned ints, this piece can be removed - odcstatevector = calibration_parts.caps_and_progress(pipeline, head_dict["odcstatevector"], odc_caps, "odc_%s" % instrument) - odctagstr = "channel-name=%s:%s, instrument=%s" % (instrument, options.dq_channel_name, instrument) - odcstatevector = pipeparts.mktaginject(pipeline, odcstatevector, odctagstr) - odcstatevectortee = pipeparts.mktee(pipeline, odcstatevector) - - # - # GAP BIT BRANCH - # - - nogap = pipeparts.mkbitvectorgen(pipeline, odcstatevectortee, threshold=1, bit_vector = 1) - nogap = pipeparts.mkcapsfilter(pipeline, nogap, odc_caps) - nogap = pipeparts.mkgeneric(pipeline, nogap, "lal_logicalundersample", required_on = 1, status_out = 512) - nogap = pipeparts.mkcapsfilter(pipeline, nogap, calibstate_caps) - - # - # OBSERVATION-INTENT BIT BRANCH - # - - obsintent = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = options.obs_intent_bitmask, status_out = 2) - obsintent = pipeparts.mkcapsfilter(pipeline, obsintent, calibstate_caps) - obsintenttee = pipeparts.mktee(pipeline, obsintent) - - # - # OBSERVATION-READY BIT BRANCH - # - - obsready = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = options.obs_ready_bitmask, status_out = 4) - obsready = pipeparts.mkcapsfilter(pipeline, obsready, calibstate_caps) - obsreadytee = pipeparts.mktee(pipeline, obsready) - - # - # H(t)-PRODUCED BIT BRANCH - # - - htproduced = pipeparts.mkbitvectorgen(pipeline, straintee, bit_vector = 8, threshold = 0) - htproduced = pipeparts.mkcapsfilter(pipeline, htproduced, "audio/x-raw, format=U32LE, rate=%d" % hoftsr) - htproduced = pipeparts.mkgeneric(pipeline, htproduced, "lal_logicalundersample", required_on = 8, status_out = 8) - htproduced = pipeparts.mkcapsfilter(pipeline, htproduced, calibstate_caps) - - # - # FILTERS-OK BIT BRANCH - # - - # Set the FILTERS-OK bit based on observation-ready transitions - filtersok = pipeparts.mkbitvectorgen(pipeline, obsintenttee, bit_vector=16, threshold=2) - filtersok = pipeparts.mkcapsfilter(pipeline, filtersok, calibstate_caps) - filtersok = calibration_parts.mkgate(pipeline, filtersok, obsreadytee, 4, long_queue, short_queue, attack_length = -int(filter_settle_time * calibstatesr), hold_length = -int(filter_latency * calibstatesr)) - filtersok = pipeparts.mkbitvectorgen(pipeline, filtersok, bit_vector = 16, nongap_is_control = True) - filtersok = pipeparts.mkcapsfilter(pipeline, filtersok, calibstate_caps) - - # - # NO-INVALID-INPUT BRANCH - # - - # Check if any of the input data channels had to be replaced by zeroes because they were < 1e-35 - resok = pipeparts.mkcapsfilter(pipeline, restee, hoft_caps) - resok = pipeparts.mkbitvectorgen(pipeline, resok, threshold=1e-35, bit_vector=1) - resok = pipeparts.mkcapsfilter(pipeline, resok, "audio/x-raw, format=U32LE, rate=%d" % hoftsr) - resok = pipeparts.mkgeneric(pipeline, resok, "lal_logicalundersample", required_on = 1, status_out = 1) - resok = pipeparts.mkcapsfilter(pipeline, resok, calibstate_caps) - if options.partial_calibration: - tstok = pipeparts.mkcapsfilter(pipeline, tsttee, ctrl_caps) - tstok = pipeparts.mkbitvectorgen(pipeline, tstok, threshold=1e-35, bit_vector=1) - tstok = pipeparts.mkcapsfilter(pipeline, tstok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) - tstok = pipeparts.mkgeneric(pipeline, tstok, "lal_logicalundersample", required_on = 1, status_out = 1) - tstok = pipeparts.mkcapsfilter(pipeline, tstok, calibstate_caps) - pumok = pipeparts.mkcapsfilter(pipeline, pumtee, ctrl_caps) - pumok = pipeparts.mkbitvectorgen(pipeline, pumok, threshold=1e-35, bit_vector=1) - pumok = pipeparts.mkcapsfilter(pipeline, pumok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) - pumok = pipeparts.mkgeneric(pipeline, pumok, "lal_logicalundersample", required_on = 1, status_out = 1) - pumok = pipeparts.mkcapsfilter(pipeline, pumok, calibstate_caps) - uimok = pipeparts.mkcapsfilter(pipeline, uimtee, ctrl_caps) - uimok = pipeparts.mkbitvectorgen(pipeline, uimok, threshold=1e-35, bit_vector=1) - uimok = pipeparts.mkcapsfilter(pipeline, uimok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) - uimok = pipeparts.mkgeneric(pipeline, uimok, "lal_logicalundersample", required_on = 1, status_out = 1) - uimok = pipeparts.mkcapsfilter(pipeline, uimok, calibstate_caps) - noinvalidinput = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [resok, short_queue], [tstok, long_queue], [pumok, long_queue], [uimok, long_queue])) - noinvalidinput = pipeparts.mkbitvectorgen(pipeline, noinvalidinput, threshold=4, bit_vector=33554432) - if options.full_calibration: - ctrlok = pipeparts.mkbitvectorgen(pipeline, darmctrltee, threshold=1e-35, bit_vector=1) - ctrlok = pipeparts.mkcapsfilter(pipeline, ctrlok, "audio/x-raw, format=U32LE, rate=%d" % ctrlsr) - ctrlok = pipeparts.mkgeneric(pipeline, ctrlok, "lal_logicalundersample", required_on = 1, status_out = 1) - ctrlok = pipeparts.mkcapsfilter(pipeline, ctrlok, calibstate_caps) - noinvalidinput = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [resok, long_queue], [ctrlok, long_queue])) - noinvalidinput = pipeparts.mkbitvectorgen(pipeline, noinvalidinput, threshold=2, bit_vector=33554432) - noinvalidinput = pipeparts.mkcapsfilter(pipeline, noinvalidinput, calibstate_caps) - noinvalidinput = pipeparts.mktee(pipeline, noinvalidinput) - # inputs that are replaced with zeros affect h(t) for a short time before and after the zeros, so we also must account for this corrupted time. - noinvalidinput = calibration_parts.mkgate(pipeline, noinvalidinput, noinvalidinput, 33554432, long_queue, short_queue, attack_length = -int(filter_settle_time * calibstatesr), hold_length = -int(filter_latency * calibstatesr)) - - # - # KAPPA-SMOOTHING-SETTLED BIT BRANCH - # - if not options.no_kappac or not options.no_kappatst or not options.no_kappapu or not options.no_fcc: - smoothingok = pipeparts.mkbitvectorgen(pipeline, obsreadytee, bit_vector=1024, threshold=4) - smoothingok = pipeparts.mkcapsfilter(pipeline, smoothingok, calibstate_caps) - smoothingok = calibration_parts.mkgate(pipeline, smoothingok, obsreadytee, 4, long_queue, short_queue, attack_length =-(median_smoothing_samples+factors_average_samples + integration_samples)) - smoothingok = pipeparts.mkbitvectorgen(pipeline, smoothingok, bit_vector = 1024, nongap_is_control = True) - smoothingok = pipeparts.mkcapsfilter(pipeline, smoothingok, calibstate_caps) - - # - # KAPPATST BITS BRANCH - # - if not options.no_kappatst: - ktstSmoothInRange, ktstMedianUncorrupt = calibration_parts.compute_kappa_bits(pipeline, smooth_ktstRtee, smooth_ktstItee, smooth_ktstRdq, smooth_ktstIdq, options.expected_kappatst_real, options.expected_kappatst_imag, options.kappatst_real_ok_var, options.kappatst_imag_ok_var, long_queue, short_queue, status_out_smooth = 2048, status_out_median = 4096, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # KAPPAPU BITS BRANCH - # - if not options.no_kappapu: - kpuSmoothInRange, kpuMedianUncorrupt = calibration_parts.compute_kappa_bits(pipeline, smooth_kpuRtee, smooth_kpuItee, smooth_kpuRdq, smooth_kpuIdq, options.expected_kappapu_real, options.expected_kappapu_imag, options.kappapu_real_ok_var, options.kappapu_imag_ok_var, long_queue, short_queue, status_out_smooth = 8192, status_out_median = 16384, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # KAPPAC BITS BRANCH - # - if not options.no_kappac: - kcSmoothInRange, kcMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_kctee, smooth_kcdq, options.expected_kappac, options.kappac_ok_var, status_out_smooth = 131072, status_out_median = 262144, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # FCC BITS BRANCH - # - if not options.no_fcc: - fccSmoothInRange, fccMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fcctee, smooth_fccdq, fcc_default, options.fcc_ok_var, status_out_smooth = 524288, status_out_median = 1048576, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # FS BITS BRANCH - # - if not options.no_fs: - fsSmoothInRange, fsMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_fs, smooth_fsdq, options.expected_fs, options.fs_ok_var, status_out_smooth = 67108864, status_out_median = 134217728, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # SRCQ BITS BRANCH - # - if not options.no_srcQ: - srcQSmoothInRange, srcQMedianUncorrupt = calibration_parts.compute_kappa_bits_only_real(pipeline, smooth_srcQ_inv, smooth_srcQdq, options.expected_srcQ, options.srcQ_ok_var, status_out_smooth = 268435456, status_out_median = 536870912, starting_rate = options.compute_factors_sr, ending_rate = calibstatesr) - - # - # COHERENCE BITS BRANCH - # - if not options.no_coherence: - if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: - pcaly_line1_coh_ok = pipeparts.mkbitvectorgen(pipeline, pcaly_line1_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = 8388608, invert_control = True) - pcaly_line1_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line1_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) - pcaly_line1_coh_ok = pipeparts.mkgeneric(pipeline, pcaly_line1_coh_ok, "lal_logicalundersample", required_on = 8388608, status_out = 8388608) - pcaly_line1_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line1_coh_ok, calibstate_caps) - - sus_coh_ok = pipeparts.mkbitvectorgen(pipeline, sus_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = 2097152, invert_control = True) - sus_coh_ok = pipeparts.mkcapsfilter(pipeline, sus_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) - sus_coh_ok = pipeparts.mkgeneric(pipeline, sus_coh_ok, "lal_logicalundersample", required_on = 2097152, status_out = 2097152) - sus_coh_ok = pipeparts.mkcapsfilter(pipeline, sus_coh_ok, calibstate_caps) - - darm_coh_ok = pipeparts.mkbitvectorgen(pipeline, darm_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = 4194304, invert_control = True) - darm_coh_ok = pipeparts.mkcapsfilter(pipeline, darm_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) - darm_coh_ok = pipeparts.mkgeneric(pipeline, darm_coh_ok, "lal_logicalundersample", required_on = 4194304, status_out = 4194304) - darm_coh_ok = pipeparts.mkcapsfilter(pipeline, darm_coh_ok, calibstate_caps) - coherence_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [pcaly_line1_coh_ok, short_queue], [sus_coh_ok, long_queue], [darm_coh_ok, long_queue])) - if not options.no_kappac or not options.no_fcc or not options.no_srcQ or not options.no_fs: - pcaly_line2_coh_ok = pipeparts.mkbitvectorgen(pipeline, pcaly_line2_coh, threshold = options.coherence_uncertainty_threshold, bit_vector = 16777216, invert_control = True) - pcaly_line2_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line2_coh_ok, "audio/x-raw, format=U32LE, rate=%d" % cohsr) - pcaly_line2_coh_ok = pipeparts.mkgeneric(pipeline, pcaly_line2_coh_ok, "lal_logicalundersample", required_on = 16777216, status_out = 16777216) - pcaly_line2_coh_ok = pipeparts.mkcapsfilter(pipeline, pcaly_line2_coh_ok, calibstate_caps) - coherence_bits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [coherence_bits, short_queue], [pcaly_line2_coh_ok, long_queue])) - - # - # H(T)-OK BIT BRANCH - # - - # First combine higher order bits to determine h(t)-OK - higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [filtersok, long_queue], [htproduced, short_queue], [obsreadytee, long_queue], [noinvalidinput, long_queue])) - htok_threshold = 28+33554432 - if options.apply_kappatst: - higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [ktstSmoothInRange, long_queue])) - htok_threshold += 2048 - if options.apply_kappapu: - higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [kpuSmoothInRange, long_queue])) - htok_threshold += 8192 - if options.apply_kappac: - higherbits = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [higherbits, short_queue], [kcSmoothInRange, long_queue])) - htok_threshold += 131072 - higherbitstee = pipeparts.mktee(pipeline, higherbits) - - # Now calculate h(t)-OK bit - htok = pipeparts.mkbitvectorgen(pipeline, higherbitstee, bit_vector = 1, threshold = htok_threshold) - htok = pipeparts.mkcapsfilter(pipeline, htok, calibstate_caps) - - # - # HW INJECTION BITS - # - - hwinjcbc = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_cbc_bitmask), status_out = 64) - hwinjcbc = pipeparts.mkcapsfilter(pipeline, hwinjcbc, calibstate_caps) - - hwinjburst = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_burst_bitmask), status_out = 128) - hwinjburst = pipeparts.mkcapsfilter(pipeline, hwinjburst, calibstate_caps) - - hwinjdetchar = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_detchar_bitmask), status_out = 256) - hwinjdetchar = pipeparts.mkcapsfilter(pipeline, hwinjdetchar, calibstate_caps) - - hwinjstoch = pipeparts.mkgeneric(pipeline, odcstatevectortee, "lal_logicalundersample", required_on = int(options.hw_inj_stoch_bitmask), status_out = 32) - hwinjstoch = pipeparts.mkcapsfilter(pipeline, hwinjstoch, calibstate_caps) - - - # - # COMBINE ALL BITS TO MAKE GDS-CALIB_STATE_VECTOR - # - - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [nogap, long_queue], [higherbitstee, short_queue], [obsintenttee, long_queue], [htok, long_queue], [hwinjcbc, long_queue], [hwinjburst, long_queue], [hwinjdetchar, long_queue], [hwinjstoch, long_queue])) - if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or not options.no_fcc: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [smoothingok, long_queue])) - if not options.no_coherence: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [coherence_bits, long_queue])) - if not options.no_kappatst: - if not options.apply_kappatst: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [ktstSmoothInRange, long_queue], [ktstMedianUncorrupt, long_queue])) - else: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [ktstMedianUncorrupt, long_queue])) - if not options.no_kappapu: - if not options.apply_kappapu: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kpuSmoothInRange, long_queue], [kpuMedianUncorrupt, long_queue])) - else: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kpuMedianUncorrupt, long_queue])) - if not options.no_kappac: - if not options.apply_kappac: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kcSmoothInRange, long_queue], [kcMedianUncorrupt, long_queue])) - else: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [kcMedianUncorrupt, long_queue])) - if not options.no_fcc: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [fccSmoothInRange, long_queue], [fccMedianUncorrupt, long_queue])) - if not options.no_fs: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [fsSmoothInRange, long_queue], [fsMedianUncorrupt, long_queue])) - if not options.no_srcQ: - calibstatevector = calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, [calibstatevector, short_queue], [srcQSmoothInRange, long_queue], [srcQMedianUncorrupt, long_queue])) - - calibstatevector = pipeparts.mkprogressreport(pipeline, calibstatevector, "progress_calibstatevec_%s" % instrument) - dqtagstr = "channel-name=%s:GDS-CALIB_STATE_VECTOR, instrument=%s" % (instrument, instrument) - calibstatevector = pipeparts.mktaginject(pipeline, calibstatevector, dqtagstr) - -# Resample the \kappa_a channels at the specified recording sample rate and change them to single precision channels -record_kappa_caps = "audio/x-raw, format=F32LE, rate=%d" % options.record_factors_sr - -# Resample the \kappa_pu channels at the specified recording sample rate and change them to single precision channels -if not options.no_kappapu: - - kpuRout = pipeparts.mkaudioconvert(pipeline, smooth_kpuRtee) - kpuRout = calibration_parts.mkresample(pipeline, kpuRout, 1, False, record_kappa_caps) - kpuRout = pipeparts.mkprogressreport(pipeline, kpuRout, "progress_kappa_pu_real_%s" % instrument) - - kpuIout = pipeparts.mkaudioconvert(pipeline, smooth_kpuItee) - kpuIout = calibration_parts.mkresample(pipeline, kpuIout, 1, False, record_kappa_caps) - kpuIout = pipeparts.mkprogressreport(pipeline, kpuIout, "progress_kappa_pu_imag_%s" % instrument) - - smooth_kpuR_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kpuR_nogate) - smooth_kpuR_nogate = calibration_parts.mkresample(pipeline, smooth_kpuR_nogate, 1, False, record_kappa_caps) - smooth_kpuR_nogate = pipeparts.mkprogressreport(pipeline, smooth_kpuR_nogate, "progress_kappa_pu_real_nogate_%s" % instrument) - - smooth_kpuI_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kpuI_nogate) - smooth_kpuI_nogate = calibration_parts.mkresample(pipeline, smooth_kpuI_nogate, 1, False, record_kappa_caps) - smooth_kpuI_nogate = pipeparts.mkprogressreport(pipeline, smooth_kpuI_nogate, "progress_kappa_pu_imag_nogate_%s" % instrument) - -# Resample the \kappa_tst channels at the specified recording sample rate and change them to single precision channels -if not options.no_kappatst: - - ktstRout = pipeparts.mkaudioconvert(pipeline, smooth_ktstRtee) - ktstRout = calibration_parts.mkresample(pipeline, ktstRout, 1, False, record_kappa_caps) - ktstRout = pipeparts.mkprogressreport(pipeline, ktstRout, "progress_kappa_tst_real_%s" % instrument) - - ktstIout = pipeparts.mkaudioconvert(pipeline, smooth_ktstItee) - ktstIout = calibration_parts.mkresample(pipeline, ktstIout, 1, False, record_kappa_caps) - ktstIout = pipeparts.mkprogressreport(pipeline, ktstIout, "progress_kappa_tst_imag_%s" % instrument) - - smooth_ktstR_nogate = pipeparts.mkaudioconvert(pipeline, smooth_ktstR_nogate) - smooth_ktstR_nogate = calibration_parts.mkresample(pipeline, smooth_ktstR_nogate, 1, False, record_kappa_caps) - smooth_ktstR_nogate = pipeparts.mkprogressreport(pipeline, smooth_ktstR_nogate, "progress_kappa_tst_real_nogate_%s" % instrument) - - smooth_ktstI_nogate = pipeparts.mkaudioconvert(pipeline, smooth_ktstI_nogate) - smooth_ktstI_nogate = calibration_parts.mkresample(pipeline, smooth_ktstI_nogate, 1, False, record_kappa_caps) - smooth_ktstI_nogate = pipeparts.mkprogressreport(pipeline, smooth_ktstI_nogate, "progress_kappa_tst_imag_nogate_%s" % instrument) - -# Resample the \kappa_c channels at the specified recording sample rate and change it to a single precision channel -if not options.no_kappac: - kcout = pipeparts.mkaudioconvert(pipeline, smooth_kctee) - kcout = calibration_parts.mkresample(pipeline, kcout, 1, False, record_kappa_caps) - kcout = pipeparts.mkprogressreport(pipeline, kcout, "progress_kappa_c_%s" % instrument) - - smooth_kc_nogate = pipeparts.mkaudioconvert(pipeline, smooth_kc_nogate) - smooth_kc_nogate = calibration_parts.mkresample(pipeline, smooth_kc_nogate, 1, False, record_kappa_caps) - smooth_kc_nogate = pipeparts.mkprogressreport(pipeline, smooth_kc_nogate, "progress_kappa_c_nogate_%s" % instrument) - -# Resample the f_cc channels at the specified recording sample rate and change it to a single precision channel -if not options.no_fcc: - fccout = pipeparts.mkaudioconvert(pipeline, smooth_fcctee) - fccout = calibration_parts.mkresample(pipeline, fccout, 1, False, record_kappa_caps) - fccout = pipeparts.mkprogressreport(pipeline, fccout, "progress_f_cc_%s" % instrument) - - smooth_fcc_nogate = pipeparts.mkaudioconvert(pipeline, smooth_fcc_nogate) - smooth_fcc_nogate = calibration_parts.mkresample(pipeline, smooth_fcc_nogate, 1, False, record_kappa_caps) - smooth_fcc_nogate = pipeparts.mkprogressreport(pipeline, smooth_fcc_nogate, "progress_f_cc_nogate_%s" % instrument) - -# Resample the f_s channels at the specified recording sample rate and change it to a single precision channel -if not options.no_fs: - fsout = pipeparts.mkaudioconvert(pipeline, smooth_fs) - fsout = calibration_parts.mkresample(pipeline, fsout, 1, False, record_kappa_caps) - fsout = pipeparts.mkprogressreport(pipeline, fsout, "progress_f_s_%s" % instrument) - - smooth_fs_nogate = pipeparts.mkaudioconvert(pipeline, smooth_fs_nogate) - smooth_fs_nogate = calibration_parts.mkresample(pipeline, smooth_fs_nogate, 1, False, record_kappa_caps) - smooth_fs_nogate = pipeparts.mkprogressreport(pipeline, smooth_fs_nogate, "progress_f_s_nogate_%s" % instrument) - -# Resample the f_s channels at the specified recording sample rate and change it to a single precision channel -if not options.no_srcQ: - srcQ_inv_out = pipeparts.mkaudioconvert(pipeline, smooth_srcQ_inv) - srcQ_inv_out = calibration_parts.mkresample(pipeline, srcQ_inv_out, 1, False, record_kappa_caps) - srcQ_inv_out = pipeparts.mkprogressreport(pipeline, srcQ_inv_out, "progress_SRC_Q_%s" % instrument) - - smooth_srcQ_inv_nogate = pipeparts.mkaudioconvert(pipeline, smooth_srcQ_inv_nogate) - smooth_srcQ_inv_nogate = calibration_parts.mkresample(pipeline, smooth_srcQ_inv_nogate, 1, False, record_kappa_caps) - smooth_srcQ_inv_nogate = pipeparts.mkprogressreport(pipeline, smooth_srcQ_inv_nogate, "progress_SRC_Q_nogate_%s" % instrument) +cleaned_strain = pipeparts.mktaginject(pipeline, cleaned_strain, cleaned_straintagstr) # # CREATE MUXER AND HOOK EVERYTHING UP TO IT @@ -1537,53 +727,12 @@ if options.frames_per_file is not None: mux.set_property("compression-scheme", options.compression_scheme) mux.set_property("compression-level", options.compression_level) -# Link the output DQ vectors up to the muxer, if applicable -if not options.no_dq_vector: - calibration_parts.mkqueue(pipeline, calibstatevector, short_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STATE_VECTOR%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, odcstatevectortee, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%s" % (instrument, options.dq_channel_name))) - strain_queue_length = long_queue -else: - strain_queue_length = short_queue # Link the strain branch to the muxer -calibration_parts.mkqueue(pipeline, strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN%s" % (instrument, chan_prefix, chan_suffix))) -if options.remove_callines: - calibration_parts.mkqueue(pipeline, cleaned_strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN_CLEAN%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the real and imaginary parts of \kappa_tst to the muxer -if not options.no_kappatst: - calibration_parts.mkqueue(pipeline, ktstRout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_REAL%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, ktstIout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_IMAGINARY%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_ktstR_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_REAL_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_ktstI_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_TST_IMAGINARY_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the real and imaginary parts of \kappa_pu to the muxer -if not options.no_kappapu: - calibration_parts.mkqueue(pipeline, kpuRout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_REAL%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, kpuIout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_IMAGINARY%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_kpuR_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_REAL_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_kpuI_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_PU_IMAGINARY_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the \kappa_c to the muxer -if not options.no_kappac: - calibration_parts.mkqueue(pipeline, kcout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_C%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_kc_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_KAPPA_C_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the f_cc to the muxer -if not options.no_fcc: - calibration_parts.mkqueue(pipeline, fccout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_CC%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_fcc_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_CC_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the f_s to the muxer -if not options.no_fs: - calibration_parts.mkqueue(pipeline, fsout, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_S%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_fs_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_F_S_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) - -# Link the src_Q to the muxer -if not options.no_srcQ: - calibration_parts.mkqueue(pipeline, srcQ_inv_out, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_SRC_Q_INVERSE%s" % (instrument, chan_prefix, chan_suffix))) - calibration_parts.mkqueue(pipeline, smooth_srcQ_inv_nogate, long_queue).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_SRC_Q_INVERSE_NOGATE%s" % (instrument, chan_prefix, chan_suffix))) +calibration_parts.mkqueue(pipeline, cleaned_strain, strain_queue_length).get_static_pad("src").link(mux.get_request_pad("%s:%sCALIB_STRAIN_CLEAN%s" % (instrument, chan_prefix, chan_suffix))) +# FIXME: Need to addd this back in when the filter latency is sorted out from the cleaning +""" # Check that all frames are long enough, that they have all of the channels by requring a certain amount of time from start-up, and that frames aren't written for times requested by the wings option def check_complete_frames(pad, info, (output_start, frame_duration, wings_start, wings_end)): buf = info.get_buffer() @@ -1609,6 +758,7 @@ if not options.no_kappatst or not options.no_kappapu or not options.no_kappac or output_start = start + max(int(filter_settle_time), options.demodulation_filter_time + options.median_smoothing_time + options.factors_averaging_time) else: output_start = start + int(filter_settle_time) +""" # If the wings option is set, need to also check that frames aren't written during the requested wing time if options.wings is not None: diff --git a/gstlal-calibration/gst/lal/Makefile.am b/gstlal-calibration/gst/lal/Makefile.am index fb2cbb39aa..c54707c73f 100644 --- a/gstlal-calibration/gst/lal/Makefile.am +++ b/gstlal-calibration/gst/lal/Makefile.am @@ -12,7 +12,8 @@ libgstlalcalibration_la_SOURCES = \ gstlal_resample.c gstlal_resample.h \ gstlal_logicalundersample.c gstlal_logicalundersample.h \ gstlal_demodulate.c gstlal_demodulate.h \ - gstlal_insertgap.c gstlal_insertgap.h + gstlal_insertgap.c gstlal_insertgap.h \ + gstlal_fccupdate.c gstlal_fccupdate.h libgstlalcalibration_la_CPPFLAGS = $(AM_CPPFLAGS) $(PYTHON_CPPFLAGS) libgstlalcalibration_la_CFLAGS = $(AM_CFLAGS) $(LAL_CFLAGS) $(GSTLAL_CFLAGS) $(gstreamer_CFLAGS) $(gstreamer_audio_CFLAGS) libgstlalcalibration_la_LDFLAGS = $(AM_LDFLAGS) $(LAL_LIBS) $(GSTLAL_LIBS) $(PYTHON_LIBS) $(gstreamer_LIBS) $(gstreamer_audio_LIBS) $(GSTLAL_PLUGIN_LDFLAGS) diff --git a/gstlal-calibration/gst/lal/gstlal_fccupdate.c b/gstlal-calibration/gst/lal/gstlal_fccupdate.c new file mode 100644 index 0000000000..17a81d513c --- /dev/null +++ b/gstlal-calibration/gst/lal/gstlal_fccupdate.c @@ -0,0 +1,680 @@ +/*Copyright (C) 2016 Kipp Cannon <kipp.cannon@ligo.org>, Theresa Chmiel, + * Madeline Wade + * + *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. + */ + +/* + *============================================================================ + * + * Preamble + * + *============================================================================ + */ + +/* + *stuff from gobject/gstreamer + */ +//#include <gstlal_fccupdate.h> +#include <stdlib.h> +#include <glib.h> +#include <gst/gst.h> +#include <gst/audio/audio.h> +#include <gst/base/gstbasetransform.h> + +//extra stuff I added + + +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_vector.h> +#include <gsl/gsl_blas.h> + +#include <gstlal/gstaudioadapter.h> +#include <gstlal/gstlal_debug.h> +#include <gstlal/gstlal.h> + +#include <complex.h> +#include <string.h> +#include <fftw3.h> +#include <math.h> + +#include <lal/LALAtomicDatatypes.h> +#include <stddef.h> + +/* + *our own stuff + */ + + +#include <gstlal_fccupdate.h> + + +/* + *============================================================================ + * + * GStreamer Boiler Plate + * + *============================================================================ + */ + +#define GST_CAT_DEFAULT gstlal_fcc_update_debug +GST_DEBUG_CATEGORY_STATIC(GST_CAT_DEFAULT); + + +G_DEFINE_TYPE_WITH_CODE( + GSTLALFccUpdate, + gstlal_fcc_update, + GST_TYPE_BASE_TRANSFORM, + GST_DEBUG_CATEGORY_INIT(GST_CAT_DEFAULT, "lal_fcc_update", 0, "lal_fcc_update element") +); + +/* + *============================================================================ + * + * Parameters + * + *============================================================================ + */ + +#define DEFAULT_FIR_MATRIX 100.0 +#define Pi 3.14159265358979323846 + +/* + *============================================================================ + * + * Utilities + * + *============================================================================ + */ +//finds the length of a 1D gsl_matrix +int num_rows(gsl_matrix * A){ + int r = A->size2; + return r; +} + +//Computes the hanning window +float *hanning(int N, short itype) +{ + int half, i, idx, n; + float *w; + + w = (float*) calloc(N, sizeof(float)); + memset(w, 0, N*sizeof(float)); + + if(itype==1) //periodic function + n = N-1; + else + n = N; + + if(n%2==0) + { + half = n/2; + for(i=0; i<half; i++) //CALC_HANNING Calculates Hanning window samples. + w[i] = 0.5 * (1 - cos(2*Pi*(i+1) / (n+1))); + + idx = half-1; + for(i=half; i<n; i++) { + w[i] = w[idx]; + idx--; + } + } + else + { + half = (n+1)/2; + for(i=0; i<half; i++) //CALC_HANNING Calculates Hanning window samples. + w[i] = 0.5 * (1 - cos(2*Pi*(i+1) / (n+1))); + + idx = half-2; + for(i=half; i<n; i++) { + w[i] = w[idx]; + idx--; + } + } + + if(itype==1) //periodic function + { + for(i=N-1; i>=1; i--) + w[i] = w[i-1]; + w[0] = 0.0; + } + return(w); +} + +//Computes a tukey window +float *tukey(int N) +{ + float * HannWindow=malloc(sizeof(float)*N/2); + HannWindow=hanning(N/2,0); + float * TukeyWindow=malloc(sizeof(float)*N); + int i=0; + int j=N/4; + for(i=0;i<N/4;i++) + TukeyWindow[i]=HannWindow[i]; + for(i=N/4;i<3*N/4;i++) + TukeyWindow[i]=1.0; + for(i=3*N/4;i<N;i++) { + TukeyWindow[i]=HannWindow[j]; + j++; + } + + return(TukeyWindow); + +} + + + +//constructs the new filter with the computed average +long double* MakeFilter(GSTLALFccUpdate *element) { + long double FcAverage=element->currentaverage; + + printf("FcAverage=%Lf\n,",FcAverage); + + long double FiltDur=(long double) element->filterduration; + long double Filtdf=1.0L/(long double)FiltDur; + long double Filtdt=1.0L/(long double) element->datarate; + int filtlength= (int) (element->datarate/(2.0*Filtdf))+1; + + int i=0; + long double Filtf[filtlength]; + while(i<filtlength) { + Filtf[i]=(i*Filtdf); + i++; + } + i=0; + long double complex CavPoleFiltForTD[filtlength]; + double instrument_cavity_pole_frequency = element->fcmodel; + + while(i<filtlength) { + CavPoleFiltForTD[i]=(1.0L+(I*(Filtf[i]/FcAverage)))/(1.0L+(I*(Filtf[i]/instrument_cavity_pole_frequency))); + i++; + } + + //adds delay to the filter + i=0; + double delaysamples=filtlength; + long double complex Delay[filtlength]; + while(i<filtlength) { + Delay[i]=cexpl(-2.0L*Pi*I*Filtf[i]*delaysamples*Filtdt); + i++; + } + i=0; + long double complex DelayedFilt[filtlength]; + while(i<filtlength) { + DelayedFilt[i]=CavPoleFiltForTD[i]*Delay[i]; + i++; + } + + //adds the negative frequencies to the filter to prepare for ifft + i=0; + long double complex NegFrequencies[filtlength]; + while(i<filtlength) { + NegFrequencies[i]=conjl(DelayedFilt[(filtlength-1-i)]); + i++; + } + + long double complex CavPoleFiltTotal[filtlength*2-2]; + i=0; + while(i<filtlength) { + CavPoleFiltTotal[i]=DelayedFilt[i]; + i++; + } + while(i<(filtlength*2)-2) { + CavPoleFiltTotal[i]=NegFrequencies[i-filtlength+1]; + i++; + } + + + // ifft (to make faster, don't create plan each time) + int N=(filtlength*2)-2; + fftw_complex *in, *out; + fftw_plan p; + in=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N); + out=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N); + p=fftw_plan_dft_1d(N,in,out,FFTW_BACKWARD,FFTW_ESTIMATE); + + i=0; + while(i<N) { + in[i]=CavPoleFiltTotal[i]; + i++; + } + + fftw_execute(p); + + i=0; + long double * CavPoleFiltTD=malloc(sizeof(long double)*N); + + while(i<N) { + CavPoleFiltTD[i]=creall(out[i]/N); + i++; + } + + + //adds Tukey window + float * TukeyWindow=malloc(sizeof(float)*N); + TukeyWindow=tukey(N); + i=0; + while(i<N) { + CavPoleFiltTD[i]=CavPoleFiltTD[i]*TukeyWindow[i]; + i++; + } + + element->fir_matrix=gsl_matrix_alloc(1,filtlength); + for(i=0;i<filtlength;i++) { + gsl_matrix_set(element->fir_matrix,0,i,CavPoleFiltTD[i]); + } + + + fftw_destroy_plan(p); + fftw_free(in); + fftw_free(out); + + +/* + //prints out CavPoleFiltTD + for(i=0;i<N;i++) { + printf("%Lf,",CavPoleFiltTD[i]); + } +*/ + + return CavPoleFiltTD; + +} + +//computes the running average +void FindAverage(GSTLALFccUpdate *element, double newpoint, int i) { + + double updatedaverage; + updatedaverage=(((i)/(i+1.0))*(element->currentaverage))+((1.0/(i+1.0))*newpoint); + element->currentaverage=updatedaverage; +} + + +//sets up signaling + +GstMessage *gstlal_fcc_update_message_fir_new(GSTLALFccUpdate *element,int filtlength) +{ + GArray *va=g_array_sized_new(FALSE,TRUE,sizeof(double),filtlength); + + int i=0; + double data; + for(i=0;i<filtlength;i++) { + data=gsl_matrix_get(element->fir_matrix,0,i); + g_array_append_val(va,data); + } + GstStructure *s = gst_structure_new( + "new_fir_matrix", + "magnitude", G_TYPE_ARRAY, va, + NULL + ); + GstMessage *m = gst_message_new_element(GST_OBJECT(element), s); + g_array_free(va,TRUE); + + printf("signal sent\n"); + + //GST_MESSAGE_TIMESTAMP(m) = XLALGPSToINT8NS(&psd->epoch); + + return m; +} + + +static void rebuild_workspace_and_reset(GObject *object) +{ + return; +} +/* +*============================================================================ +* +* Signals +* +*============================================================================ +*/ + + + +/* + *============================================================================ + * + * GstBaseTransform Method Overrides + * + *============================================================================ + */ + + +static GstFlowReturn transform_ip(GstBaseTransform *trans, GstBuffer *buf) { + + GSTLALFccUpdate *element = GSTLAL_FCC_UPDATE(trans); + GstMapInfo mapinfo; + GstFlowReturn result = GST_FLOW_OK; + + + GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP); + gst_buffer_map(buf, &mapinfo, GST_MAP_READ); + g_assert(mapinfo.size % sizeof(gdouble) == 0); + + gdouble *data, *end; + data = (gdouble*) mapinfo.data; + end = (gdouble*) (mapinfo.data+mapinfo.size); + int averaging_length=floor(element->averaging_time*element->fccrate); + int i=element->index; + + while(data<end) { + + if(i<averaging_length) { + //continuously updated the average fcc_filter value + FindAverage(element,*data,i); + i++; + *data++; + + } + else { + printf("reached update length\n"); + //makes the new fir_matrix using the current average fcc_filter value + int filtlength=((int)(element->datarate*element->filterduration/2.0))*2; + long double *FccFilter=MakeFilter(element); + + //updated the values of the fir_matrix + element->fir_matrix=gsl_matrix_alloc(1,filtlength); + for(i=0;i<filtlength;i++) { + gsl_matrix_set(element->fir_matrix,0,i,FccFilter[i]); + } + + //sends a signal that the fir-matrix has been updated + int messagesent; + g_object_notify(G_OBJECT(element), "fir-matrix"); + messagesent=gst_element_post_message( + GST_ELEMENT(element), + gstlal_fcc_update_message_fir_new(element,filtlength) + ); + + i=0; + } + + } + element->index=i; + + gst_buffer_unmap(buf,&mapinfo); + return result; +} + + +/* + *============================================================================ + * _ + * GObject Method Overrides + * + *============================================================================ + */ + + +enum property { + FIR_MATRIX = 1, + DATA_RATE, + FCC_RATE, + FILTER_DURATION, + FC_MODEL, + AVERAGING_TIME +}; + + +//Set FIR_MATRIX property + +static void set_property(GObject *object, enum property prop_id, const GValue *value, GParamSpec *pspec) +{ + GSTLALFccUpdate *element = GSTLAL_FCC_UPDATE(object); + + GST_OBJECT_LOCK(element); + + switch (prop_id) { + case FIR_MATRIX: + if(element->fir_matrix) + gsl_matrix_free(element->fir_matrix); + element->fir_matrix = gstlal_gsl_matrix_from_g_value_array(g_value_get_boxed(value)); + break; + case DATA_RATE: + element->datarate=g_value_get_int(value); + break; + case FCC_RATE: + element->fccrate=g_value_get_int(value); + break; + case FILTER_DURATION: + element->filterduration=g_value_get_double(value); + break; + case FC_MODEL: + element->fcmodel = g_value_get_double(value); + break; + case AVERAGING_TIME: + element->averaging_time = g_value_get_double(value); + break; + default: + G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +//get FIR_MATRIX property + +static void get_property(GObject *object, enum property prop_id, GValue *value, GParamSpec *pspec) +{ + GSTLALFccUpdate *element = GSTLAL_FCC_UPDATE(object); + + GST_OBJECT_LOCK(element); + + switch (prop_id) { + case FIR_MATRIX: + if(element->fir_matrix) + g_value_take_boxed(value, gstlal_g_value_array_from_gsl_matrix(element->fir_matrix)); + break; + case DATA_RATE: + g_value_set_int(value,element->datarate); + break; + case FCC_RATE: + g_value_set_int(value,element->fccrate); + break; + case FILTER_DURATION: + g_value_set_double(value,element->filterduration); + break; + case FC_MODEL: + g_value_set_double(value, element->fcmodel); + break; + case AVERAGING_TIME: + g_value_set_double(value, element->averaging_time); + break; + default: + G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec); + break; + } + + GST_OBJECT_UNLOCK(element); +} + +//finalize + +static void finalize(GObject *object) +{ + GSTLALFccUpdate *element = GSTLAL_FCC_UPDATE(object); + if(element->fir_matrix) { + gsl_matrix_free(element->fir_matrix); + element->fir_matrix = NULL; + } + G_OBJECT_CLASS(gstlal_fcc_update_parent_class)->finalize(object); +} + + + + +//Class_init() +static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE( + GST_BASE_TRANSFORM_SINK_NAME, + GST_PAD_SINK, + GST_PAD_ALWAYS, + GST_STATIC_CAPS( + "audio/x-raw, " \ + "rate = " GST_AUDIO_RATE_RANGE ", " \ + "channels = (int) 1, " \ + "format = (string) {" GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) "}, " \ + "layout = (string) interleaved, " \ + "channel-mask = (bitmask) 0" + ) +); + + +static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE( + GST_BASE_TRANSFORM_SRC_NAME, + GST_PAD_SRC, + GST_PAD_ALWAYS, + GST_STATIC_CAPS( + GST_AUDIO_CAPS_MAKE("{" GST_AUDIO_NE(F32) ", " GST_AUDIO_NE(F64) "}") ", " \ + "layout = (string) interleaved, " \ + "channel-mask = (bitmask) 0" + ) +); + + + +static void gstlal_fcc_update_class_init(GSTLALFccUpdateClass *klass) +{ + + GstBaseTransformClass *transform_class = GST_BASE_TRANSFORM_CLASS(klass); + GstElementClass *element_class = GST_ELEMENT_CLASS(klass); + GObjectClass *gobject_class = G_OBJECT_CLASS(klass); + + gst_element_class_set_details_simple( + element_class, + "Update the Fcc Filter", + "Filter/Audio", + "Makes a new fcc filter based on a new cavity pole frequency value", + "Theresa Chmiel <theresa.chmiel@ligo.org>" + ); + + gobject_class->set_property = GST_DEBUG_FUNCPTR(set_property); + gobject_class->get_property = GST_DEBUG_FUNCPTR(get_property); + gobject_class->finalize = GST_DEBUG_FUNCPTR(finalize); + + + g_object_class_install_property( + gobject_class, + FIR_MATRIX, + g_param_spec_value_array( + "fir-matrix", + "FIR Matrix", + "Array of the cavity pole filter information in the time domain", + g_param_spec_value_array( + "response", + "Impulse Response", + "Array of amplitudes.", + g_param_spec_double( + "amplitude", + "Amplitude", + "Impulse response sample", + -G_MAXDOUBLE, G_MAXDOUBLE, DEFAULT_FIR_MATRIX, + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS + ), + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS + ), + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | GST_PARAM_CONTROLLABLE + ) + ); + + + g_object_class_install_property( + gobject_class, + DATA_RATE, + g_param_spec_int( + "data-rate", + "Data rate", + "The rate of the incoming data to be filtered (not the incoming fcc data).", + 0, + G_MAXINT, + 16384, //Default rate + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); + + g_object_class_install_property( + gobject_class, + FCC_RATE, + g_param_spec_int( + "fcc-rate", + "Fcc sample rate", + "The rate of the incoming fcc data (not the incoming data to be filtered).", + 0, + G_MAXINT, + 16, //Default rate + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); + + g_object_class_install_property( + gobject_class, + FILTER_DURATION, + g_param_spec_double( + "filter-duration", + "Filter duration", + "The the length of the desired filter to be generated in seconds.", + 0.0, + G_MAXDOUBLE, + 0.01, //Default filter duration + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); + + g_object_class_install_property( + gobject_class, + FC_MODEL, + g_param_spec_double( + "fcc-model", + "F_cc model value", + "The cavity pole frequency value from the static calibration model.", + 0.0, + G_MAXDOUBLE, + 360.0, + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); + + g_object_class_install_property( + gobject_class, + AVERAGING_TIME, + g_param_spec_double( + "averaging-time", + "Averaging time", + "The amount of time to averaging computed f_c values before constructing a new FIR filter.", + 1.0, + G_MAXDOUBLE, + 1024.0, + G_PARAM_READWRITE | G_PARAM_STATIC_STRINGS | G_PARAM_CONSTRUCT + ) + ); + + gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&src_factory)); + gst_element_class_add_pad_template(element_class, gst_static_pad_template_get(&sink_factory)); + + transform_class->transform_ip = GST_DEBUG_FUNCPTR(transform_ip); +} + + +//init + +static void gstlal_fcc_update_init(GSTLALFccUpdate *element) +{ + g_signal_connect(G_OBJECT(element), "notify::fir-matrix",G_CALLBACK(rebuild_workspace_and_reset), NULL); + gst_base_transform_set_gap_aware(GST_BASE_TRANSFORM(element), TRUE); + gst_base_transform_set_prefer_passthrough (GST_BASE_TRANSFORM(element), TRUE); +} + + + diff --git a/gstlal-calibration/gst/lal/gstlal_fccupdate.h b/gstlal-calibration/gst/lal/gstlal_fccupdate.h new file mode 100644 index 0000000000..676705c17b --- /dev/null +++ b/gstlal-calibration/gst/lal/gstlal_fccupdate.h @@ -0,0 +1,111 @@ +/* + *Copyright (C) 2016 Kipp Cannon <kipp.cannon@ligo.org>, Theresa Chmiel, Madeline Wade + * + *Permission is hereby granted, free of charge, to any person obtaining a + *copy of this software and associated documentation files (the + *"Software"), to deal in the Software without restriction, including + *without limitation the rights to use, copy, modify, merge, publish, + *distribute, sublicense, and/or sell copies of the Software, and to + *permit persons to whom the Software is furnished to do so, subject to + *the following conditions: + * + *The above copyright notice and this permission notice shall be included + *in all copies or substantial portions of the Software. + * + *THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + *OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + *MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + *IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + *CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + *TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + *SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * + *Alternatively, the contents of this file may be used under the GNU + *Lesser General Public License Version 2.1 (the "LGPL"), in which case + *the following provisions apply instead of the ones mentioned above: + * + *This library is free software; you can redistribute it and/or modify it + *under the terms of the GNU Library General Public License as published + *by the Free Software Foundation; either version 2 of the License, or (at + *your option) any later version. + * + *This library 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 + *Library General Public License for more details. + * + *You should have received a copy of the GNU Library General Public + *License along with this library; if not, write to the Free Software + *Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, + *USA. + */ + +#ifndef __GST_LAL_FCC_UPDATE_H__ +#define __GST_LAL_FCC_UPDATE_H__ + + +#include <glib.h> +#include <gst/gst.h> +#include <gst/base/gstbasetransform.h> +#include <gsl/gsl_matrix.h> + +#include <gst/audio/audio.h> +#include <gstlal/gstaudioadapter.h> +#include <complex.h> +#include <fftw3.h> + +G_BEGIN_DECLS + +#define GSTLAL_FCC_UPDATE_TYPE \ + (gstlal_fcc_update_get_type()) +#define GSTLAL_FCC_UPDATE(obj) \ + (G_TYPE_CHECK_INSTANCE_CAST((obj), GSTLAL_FCC_UPDATE_TYPE, GSTLALFccUpdate)) +#define GSTLAL_FCC_UPDATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_CAST((klass), GSTLAL_FCC_UPDATE_TYPE, GSTLALFccUpdateClass)) +#define GST_IS_GSTLAL_FCC_UPDATE(obj) \ + (G_TYPE_CHECK_INSTANCE_TYPE((obj), GSTLAL_FCC_UPDATE_TYPE)) +#define GST_IS_GSTLAL_FCC_UPDATE_CLASS(klass) \ + (G_TYPE_CHECK_CLASS_TYPE((klass), GSTLAL_FCC_UPDATE_TYPE)) + + + +typedef struct _GSTLALFccUpdate GSTLALFccUpdate; +typedef struct _GSTLALFccUpdateClass GSTLALFccUpdateClass; + + + +struct _GSTLALFccUpdate { + GstBaseTransform element; + gsl_matrix *fir_matrix; + int index; + double currentaverage; + gint datarate; + gint fccrate; + gdouble filterduration; + gdouble fcmodel; + gdouble averaging_time; +}; + + + + + +struct _GSTLALFccUpdateClass { + GstBaseTransformClass parent_class; +}; + + +GType gstlal_fcc_update_get_type(void); + +G_END_DECLS + + +#endif /* __GST_LAL_FCC_UPDATE_H__ */ + + + + + + + + diff --git a/gstlal-calibration/gst/lal/gstlalcalibration.c b/gstlal-calibration/gst/lal/gstlalcalibration.c index 12a2969b1e..ef299ab4ff 100644 --- a/gstlal-calibration/gst/lal/gstlalcalibration.c +++ b/gstlal-calibration/gst/lal/gstlalcalibration.c @@ -61,6 +61,7 @@ #include <gstlal_logicalundersample.h> #include <gstlal_demodulate.h> #include <gstlal_insertgap.h> +#include <gstlal_fccupdate.h> /* @@ -89,6 +90,7 @@ static gboolean plugin_init(GstPlugin *plugin) {"lal_logicalundersample", GSTLAL_LOGICALUNDERSAMPLE_TYPE}, {"lal_demodulate", GSTLAL_DEMODULATE_TYPE}, {"lal_insertgap", GSTLAL_INSERTGAP_TYPE}, + {"lal_fcc_update", GSTLAL_FCC_UPDATE_TYPE}, {NULL, 0}, }; -- GitLab