Commit 3a1b20c7 authored by Aaron Viets's avatar Aaron Viets
Browse files

gstlal-calibration: ASD plotting tools

parent 1261a767
Pipeline #213137 passed with stages
in 37 minutes and 47 seconds
......@@ -1062,6 +1062,33 @@ def DolphChebyshev(N, alpha, return_double = False):
return win / win[N // 2]
def Blackman(N, return_double = False):
win = np.zeros(N, dtype = np.float128)
for i in range(N):
win[i] = 0.42 - 0.5 * np.cos((two_pi * i) / (N - 1)) + 0.08 * np.cos((2 * two_pi * i) / (N - 1))
if return_double:
return np.float64(win)
else:
return win
def Hann(N, return_double = False):
win = np.zeros(N, dtype = np.float128)
for i in range(N):
a = np.sin((pi * i) / (N - 1))
win[i] = a * a
if return_double:
return np.float64(win)
else:
return win
#
# A resampler
#
......@@ -1215,3 +1242,42 @@ def freqresp(filt, delay_samples = 0, samples_per_lobe = 8, return_double = Fals
return rfft(filt_prime, return_double = return_double)
#
# Take an ASD of a time series
#
def asd(data, sr, fft_samples, fft_spacing, window = 'blackman', freq_res = 1.0):
# How many FFTs will we take?
if len(data) < fft_samples:
print("WARNING: Not enough data to take an ASD of requested length")
fft_samples = len(data)
num_ffts = (len(data) - fft_samples) // fft_spacing + 1
# Window function
if window == 'dpss' or window == 'DPSS' or window == 'slepian' or window == 'Slepian':
alpha = max(freq_res / sr, 1.0 / fft_samples) * fft_samples
win = DPSS(fft_samples, alpha)
elif window == 'kaiser' or window == 'Kaiser':
alpha = max(freq_res / sr, 1.0 / fft_samples) * fft_samples
win = kaiser(fft_samples, pi * alpha)
elif window == 'DC' or window == 'DolphChebyshev' or window == 'dolph_chebyshev' or window == 'dolphchebyshev' or window == 'Dolph_Chebyshev':
alpha = max(freq_res / sr, 1.0 / fft_samples) * fft_samples
win = DolphChebyshev(fft_samples, pi * alpha)
elif window == 'Hann' or window == 'hann' or window == 'Hanning' or window == 'hanning':
win = Hann(fft_samples)
else:
if window != 'Blackman' and window != 'blackman':
print("WARNING: Unrecognized window type '%s'. Using Blackman window." % window)
win = Blackman(fft_samples)
# Compute the ASD
asd = abs(rfft(win * data[:fft_samples]))
for i in range(1, num_ffts):
asd += abs(rfft(win * data[i * fft_spacing : i * fft_spacing + fft_samples]))
asd /= num_ffts * sr * np.sqrt((float(fft_samples) / sr))
return np.float64(asd)
......@@ -8,11 +8,11 @@ IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O3
START = $(shell echo 1256655618 | bc)
START = $(shell echo 1263106818 | bc)
# 1237831461
# 1238288418
# 1269090018
END = $(shell echo 1269367218 | bc)
END = $(shell echo 1263116818 | bc)
# 1269367218
# 1238338818
# 1269133218
......@@ -35,6 +35,8 @@ DCSEXACTKAPPASCONFIGS = Filters/O3/GDSFilters/H1DCS_test_1237831461_exactKappas.
DCSAPPROXKAPPASCONFIGS = Filters/O3/GDSFilters/H1DCS_test_1237831461_approxKappas.ini
#Filters/O3/GDSFilters/H1DCS_test_1237831461_approxKappas.ini
#Filters/O3/GDSFilters/H1DCS_test_1256655618_v2_approxKappas.ini
DCSSTATSUBCONFIGS = Filters/O3/GDSFilters/H1DCS_test_1256655618_v2_stationarySubtraction.ini
DCSNONSTATSUBCONFIGS = Filters/O3/GDSFilters/H1DCS_test_1256655618_v2_nonstationarySubtraction.ini
DCSO31CONFIGS = Filters/O3/GDSFilters/H1DCS_KappaErrorTest_1237831461.ini
DCSO32CONFIGS = Filters/O3/GDSFilters/H1DCS_KappaErrorTest_1244307456.ini
DCSO33CONFIGS = Filters/O3/GDSFilters/H1DCS_KappaErrorTest_1251057623.ini
......@@ -44,7 +46,8 @@ DCSTESTCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1252083618_AllCorrectionsTest.i
GDSSHMCONFIGS = Filters/O3/GDSFilters/H1GDS_1258216456_testing.ini
GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini
all: approxKappasErrorPlot
all: noise_subtraction_ASD_DCS
#$(IFO)1_hoft_DCS_STATSUB_frames.cache $(IFO)1_hoft_DCS_NONSTATSUB_frames.cache
###############################################
### These commands should change less often ###
......@@ -119,6 +122,14 @@ $(IFO)1_hoft_DCS_APPROXKAPPAS_frames.cache: $(IFO)1_easy_raw_frames.cache filter
#123828
#126912
$(IFO)1_hoft_DCS_STATSUB_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(DCSSTATSUBCONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS_STATSUB-*.gwf | lalapps_path2cache > $@
$(IFO)1_hoft_DCS_NONSTATSUB_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(DCSNONSTATSUBCONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS_NONSTATSUB-*.gwf | lalapps_path2cache > $@
$(IFO)1_hoft_DCS_TDCFUNC_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(DCS_TDCF_UNCERTAINTY_CONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS_TDCFUNC-*.gwf | lalapps_path2cache > $@
......@@ -283,8 +294,11 @@ kappastimeseries: $(IFO)1_hoft_DCS_frames.cache $(IFO)1_C01_frames.cache
SRC_detuning_C00: $(IFO)1_C00_frames.cache
python3 timeserieskappas.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_C00_frames.cache --channel-list 'GDS-CALIB_F_S_SQUARED,GDS-CALIB_SRC_Q_INVERSE'
noise_subtraction_ASD_DCS: $(IFO)1_ALL_CORR_frames.cache
./ASD_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --frame-cache-list '$(IFO)1_ALL_CORR_frames.cache,$(IFO)1_ALL_CORR_frames.cache' --channel-list 'DCS-CALIB_STRAIN,DCS-CALIB_STRAIN_CLEAN' --ASD-fmin 1 --ASD-fmax 8192 --ratio-fmin 10 --ratio-fmax 1000 --ASD-min 1e-24 --ASD-max 1e-18 --ratio-min 0
noise_subtraction_ASD_DCS: $(IFO)1_hoft_DCS_STATSUB_frames.cache
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_hoft_DCS_STATSUB_frames.cache,$(IFO)1_hoft_DCS_STATSUB_frames.cache --channel-list DCS-CALIB_STRAINSTATSUB,DCS-CALIB_STRAIN_CLEANSTATSUB --sample-rate 16384 --fft-time 32 --fft-spacing 16 --freq-res 0.2 --frequency-max=90 --labels 'No Cleaning,Stationary Sub'
noise_channels_ASD_DCS:
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_easy_raw_frames.cache,$(IFO)1_easy_raw_frames.cache,$(IFO)1_easy_raw_frames.cache,$(IFO)1_easy_raw_frames.cache,$(IFO)1_easy_raw_frames.cache,$(IFO)1_easy_raw_frames.cache --channel-list ASC-DHARD_P_OUT_DQ,ASC-DHARD_Y_OUT_DQ,ASC-CHARD_P_OUT_DQ,ASC-CHARD_Y_OUT_DQ,ASC-CSOFT_P_OUT_DQ,ASC-DSOFT_P_OUT_DQ --sample-rate 512 --fft-time 8 --fft-spacing 4 --freq-res 1 --frequency-max=256 --labels 'ASC-DHARD_P,ASC-DHARD_Y,ASC-CHARD_P,ASC-CHARD_Y,ASC-CSOFT_P,ASC-DSOFT_P' --ASD-min 0.0001 --ASD-max 10
noise_subtraction_ASD_DCS_LINES: $(IFO)1_hoft_DCS_frames.cache
./ASD_plots --ifo $(IFO)1 --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --frame-cache-list '$(IFO)1_hoft_DCS_frames.cache,$(IFO)1_hoft_DCS_frames.cache' --channel-list 'DCS-CALIB_STRAIN,DCS-CALIB_STRAIN_CLEAN' --ASD-fmin 0.3
......
......@@ -240,6 +240,28 @@ if not any(channel_list):
channel_list.append("ASC-DHARD_Y_OUT_DQ")
channel_list.append("ASC-CHARD_P_OUT_DQ")
channel_list.append("ASC-CHARD_Y_OUT_DQ")
channel_list.append("ASC-CSOFT_P_OUT_DQ")
channel_list.append("ASC-DSOFT_P_OUT_DQ")
channel_list.append("ASC-INP1_P_INMON")
channel_list.append("ASC-INP1_Y_INMON")
channel_list.append("ASC-MICH_P_INMON")
channel_list.append("ASC-MICH_Y_INMON")
channel_list.append("ASC-PRC1_P_INMON")
channel_list.append("ASC-PRC1_Y_INMON")
channel_list.append("ASC-PRC2_P_INMON")
channel_list.append("ASC-PRC2_Y_INMON")
channel_list.append("ASC-SRC1_P_INMON")
channel_list.append("ASC-SRC1_Y_INMON")
channel_list.append("ASC-SRC2_P_INMON")
channel_list.append("ASC-SRC2_Y_INMON")
channel_list.append("ASC-DHARD_P_INMON")
channel_list.append("ASC-DHARD_Y_INMON")
channel_list.append("ASC-CHARD_P_INMON")
channel_list.append("ASC-CHARD_Y_INMON")
channel_list.append("ASC-DSOFT_P_INMON")
channel_list.append("ASC-DSOFT_Y_INMON")
channel_list.append("ASC-CSOFT_P_INMON")
channel_list.append("ASC-CSOFT_Y_INMON")
channel_list.append("CAL-CS_TDEP_PCAL_LINE1_UNCERTAINTY")
channel_list.append("CAL-CS_TDEP_PCAL_LINE2_REF_A_PUM_IMAG")
channel_list.append("CAL-CS_TDEP_PCAL_LINE2_REF_A_PUM_REAL")
......
#!/usr/bin/env python3
#
# Copyright (C) 2021 Aaron Viets
#
# 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
#
# =============================================================================
#
import matplotlib; matplotlib.use('Agg')
import numpy as np
import matplotlib.patches as mpatches
from matplotlib import rc
rc('text', usetex = True)
matplotlib.rcParams['font.family'] = 'Times New Roman'
matplotlib.rcParams['font.size'] = 22
matplotlib.rcParams['legend.fontsize'] = 18
matplotlib.rcParams['mathtext.default'] = 'regular'
import matplotlib.pyplot as plt
from optparse import OptionParser, Option
import gi
gi.require_version('Gst', '1.0')
from gi.repository import GObject, Gst
GObject.threads_init()
Gst.init(None)
import lal
from lal import LIGOTimeGPS
from gstlal import pipeparts
from gstlal import calibration_parts
from gstlal import simplehandler
from gstlal import datasource
from ligo import segments
from gstlal import test_common
from gstlal import FIRtools as fir
from ticks_and_grid import ticks_and_grid
parser = OptionParser()
parser.add_option("--gps-start-time", metavar = "seconds", type = int, help = "GPS time at which to start processing data")
parser.add_option("--gps-end-time", metavar = "seconds", type = int, help = "GPS time at which to stop processing data")
parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the interferometer (IFO), e.g., H1, L1")
parser.add_option("--frame-cache-list", metavar = "name", type = str, help = "Comma-separated list of frame cache files that contain numerators of transfer functions")
parser.add_option("--channel-list", metavar = "list", type = str, default = None, help = "Comma-separated list of channel-names of numerators")
parser.add_option("--zeros", metavar = "list", type = str, default = None, help = "Semicolon-separated list of comma-separated lists of real and imaginary parts of zeros with which to filter the ASD.")
parser.add_option("--poles", metavar = "list", type = str, default = None, help = "Semicolon-separated list of comma-separated list of real and imaginary parts of poles with which which to filter the ASD.")
parser.add_option("--gain", type = str, default = None, help = "Gain factor(s) to apply to the ASD")
parser.add_option("--sample-rate", metavar = "Hz", type = str, default = '16384', help = "Sample rate(s) of input data.")
parser.add_option("--fft-time", metavar = "seconds", type = float, default = 4, help = "Duration in seconds of FFTs used to compute ASD")
parser.add_option("--fft-spacing", metavar = "seconds", type = float, default = 2, help = "Spacing in seconds of FFTs used to compute ASD")
parser.add_option("--window", type = str, default = 'blackman', help = "Which window function to use to window the data when taking FFTs. Options are 'blackman', (default), 'kaiser', 'dpss', 'dolph_chebyshev', and 'hann'.")
parser.add_option("--freq-res", metavar = "Hz", type = float, default = 1.0, help = "Frequency resolution of ASD. Only applies when using Kaiser, DPSS, and Dolph Chebyshev windows.")
parser.add_option("--frequency-min", type = float, default = 10, help = "Minimum frequency for plot")
parser.add_option("--frequency-max", type = float, default = 8192, help = "Maximum frequency for plot")
parser.add_option("--ASD-min", type = float, default = 1e-24, help = "Minimum for ASD plot")
parser.add_option("--ASD-max", type = float, default = 1e-18, help = "Maximum for ASD plot")
parser.add_option("--labels", metavar = "list", type = str, default = None, help = "Comma-separated list of labels corresponding to each ASD, to be added to plot legend.")
options, filenames = parser.parse_args()
# Get frame cache list and channel list
ifo = options.ifo
frame_cache_list = options.frame_cache_list.split(',')
channel_list = options.channel_list.split(',')
chan_list = []
for i in range(len(channel_list)):
chan_list.append((ifo, channel_list[i]))
labels = options.labels.split(',')
if len(frame_cache_list) != len(channel_list):
raise ValueError('Number of frame cache files must equal number of channels: %d != %d' % (len(frame_cache_list), len(channel_list)))
if len(labels) != len(channel_list):
raise ValueError('Number of labels must equal number of channels: %d != %d' % (len(labels), len(channel_list)))
sample_rate = options.sample_rate.split(',')
sr = []
for i in range(len(labels)):
if i < len(sample_rate):
sr.append(int(sample_rate[i]))
else:
sr.append(sr[-1])
max_sr = 1
while max_sr < options.frequency_max * 2.1:
max_sr *= 2
for i in range(len(sr)):
sr[i] = sr[i] if sr[i] < max_sr else max_sr
#
# =============================================================================
#
# Pipelines
#
# =============================================================================
#
#
# This pipeline reads in time series data and writes it to file.
#
def read_data(pipeline, name):
for i in range(len(labels)):
data = pipeparts.mklalcachesrc(pipeline, location = frame_cache_list[i], cache_dsc_regex = ifo)
data = pipeparts.mkframecppchanneldemux(pipeline, data, do_file_checksum = False, skip_bad_files = True, channel_list = list(map("%s:%s".__mod__, chan_list)))
data = calibration_parts.hook_up(pipeline, data, channel_list[i], ifo, 1.0, element_name_suffix = "%d" % i)
data = calibration_parts.caps_and_progress(pipeline, data, "audio/x-raw,format=F64LE", labels[i])
data = calibration_parts.mkresample(pipeline, data, 5, False, sr[i])
pipeparts.mknxydumpsink(pipeline, data, '%s_%s_ASD_data_%d-%d.txt' % (ifo, labels[i].replace(' ', '_'), options.gps_start_time, options.gps_end_time - options.gps_start_time))
#
# done
#
return pipeline
#
# =============================================================================
#
# Main
#
# =============================================================================
#
# Run pipeline
test_common.build_and_run(read_data, "read_data", segment = segments.segment((LIGOTimeGPS(0, 1000000000 * options.gps_start_time), LIGOTimeGPS(0, 1000000000 * options.gps_end_time))))
# Frequency vector
fft_samples = int(options.fft_time * max(sr))
fvec = np.arange(0, max(sr) / 2.0 + max(sr) / 4.0 / (fft_samples // 2), max(sr) / 2.0 / (fft_samples // 2))
# Compute ASDs
ASDs = []
for i in range(len(labels)):
data = np.loadtxt('%s_%s_ASD_data_%d-%d.txt' % (ifo, labels[i].replace(' ', '_'), options.gps_start_time, options.gps_end_time - options.gps_start_time))
data = np.transpose(data)[1]
ASDs.append(fir.asd(data, sr[i], int(options.fft_time * sr[i]), int(options.fft_spacing * sr[i]), window = options.window, freq_res = options.freq_res))
# Get any zeros and poles that we want to use to filter the ASD
if options.zeros is not None:
zeros = []
strzeros = options.zeros.split(';')
while(len(strzeros) < len(labels)):
strzeros.append(strzeros[-1])
for i in range(len(strzeros)):
strzeros[i] = strzeros[i].split(',')
zeros.append([])
for j in range(len(strzeros[i]) // 2):
z = float(strzeros[i][2 * j]) + 1j * float(strzeros[i][2 * j + 1])
for k in range(len(ASDs[i])):
ASDs[i][k] = ASDs[i][k] * (1.0 + 1j * fvec[k] / z)
if options.poles is not None:
poles = []
strpoles = options.poles.split(';')
while(len(strpoles) < len(labels)):
strpoles.append(strpoles[-1])
for i in range(len(strpoles)):
strpoles[i] = strpoles[i].split(',')
poles.append([])
for j in range(len(strpoles[i]) // 2):
p = float(strpoles[i][2 * j]) + 1j * float(strpoles[i][2 * j + 1])
for k in range(len(ASDs[i])):
ASDs[i][k] = ASDs[i][k] / (1.0 + 1j * fvec[k] / p)
if options.gain is not None:
gain = options.gain.split(',')
for i in range(len(gain)):
gain[i] = float(gain[i])
while len(gain) < len(labels):
gain.append(gain[-1])
for i in range(len(labels)):
ASDs[i] *= gain[i]
# Decide a color scheme.
colors = []
# Start with defaults to use if we don't recognize what we are plotting.
default_colors = ['red', 'limegreen', 'gold', 'orange', 'aqua', 'magenta', 'blue']
# Use specific colors for known version of calibration
C02_labels = ['C02', 'c02']
C01_labels = ['C01', 'c01', 'DCS', 'dcs', 'high-latency', 'High-Latency', 'High-latency', 'HIGH-LATENCY', 'high_latency', 'High_Latency', 'High_latency', 'HIGH_LATENCY', 'high latency', 'High Latency', 'High latency', 'HIGH LATENCY', 'offline', 'Offline', 'OFFLINE']
C00_labels = ['C00', 'c00', 'GDS', 'gds', 'low-latency', 'Low-Latency', 'Low-latency', 'LOW-LATENCY', 'low_latency', 'Low_Latency', 'Low_latency', 'LOW_LATENCY', 'low latency', 'Low Latency', 'Low latency', 'LOW LATENCY', 'online', 'Online', 'ONLINE']
CALCS_labels = ['CALCS', 'calcs', 'Calcs', 'CalCS', 'CAL-CS', 'Cal-CS', 'cal-cs', 'CAL_CS', 'Cal_CS', 'cal_cs', 'front', 'FRONT', 'Front']
for i in range(len(labels)):
if any ([x for x in C02_labels if x in labels[i]]):
colors.append('springgreen')
elif any ([x for x in C01_labels if x in labels[i]]):
colors.append('maroon')
elif any ([x for x in C00_labels if x in labels[i]]):
colors.append('royalblue')
elif any([x for x in CALCS_labels if x in labels[i]]):
colors.append('silver')
else:
colors.append(default_colors[i % 7])
# Make plots
freq_scale = 'log' if options.frequency_min > 0.0 and options.frequency_max / options.frequency_min > 10 else 'linear'
ASD_scale = 'log' if options.ASD_min > 0.0 and options.ASD_max / options.ASD_min > 10 else 'linear'
plt.figure(figsize = (10, 6))
for i in range(0, len(labels)):
plt.plot(fvec[:len(ASDs[i])], ASDs[i], colors[i % 6], linewidth = 0.75) #, label = r'$%s$' % labels[i])
if i == 0:
patches = [mpatches.Patch(color = colors[j], label = r'$%s$' % labels[j]) for j in range(len(labels))]
plt.legend(handles = patches, loc = 'upper right', ncol = 1)
#leg = plt.legend(fancybox = True, loc = "lower right")
#leg.get_frame().set_alpha(0.8)
if i == 0:
#plt.title(r'${\rm %s} \ %s \ / \ %s$' % (ifo, options.name, options.denominator_name))
plt.ylabel(r'${\rm ASD}\ \left[{\rm strain / }\sqrt{\rm Hz}\right]$')
plt.xlabel(r'${\rm Frequency \ [Hz]}$')
ticks_and_grid(plt.gca(), xmin = options.frequency_min, xmax = options.frequency_max, ymin = options.ASD_min, ymax = options.ASD_max, xscale = freq_scale, yscale = ASD_scale)
plt.savefig('%s_%s_ASD_%d-%d.png' % (ifo, labels[i].replace(' ', '_'), options.gps_start_time, options.gps_end_time - options.gps_start_time))
plt.savefig('%s_%s_ASD_%d-%d.pdf' % (ifo, labels[i].replace(' ', '_'), options.gps_start_time, options.gps_end_time - options.gps_start_time))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment