From e73569d6d037a73d8783db280462a41d5856f068 Mon Sep 17 00:00:00 2001
From: Aaron Viets <aaron.viets@ligo.org>
Date: Thu, 26 Sep 2024 18:34:17 -0700
Subject: [PATCH] plot_spectrogram.py: better default behavior for plot
parameters
---
tests/check_calibration/Makefile | 6 ++--
tests/check_calibration/plot_spectrogram.py | 39 ++++++++++++++-------
2 files changed, 30 insertions(+), 15 deletions(-)
diff --git a/tests/check_calibration/Makefile b/tests/check_calibration/Makefile
index 8b6146445..1a58247c1 100644
--- a/tests/check_calibration/Makefile
+++ b/tests/check_calibration/Makefile
@@ -10,7 +10,7 @@ OBSRUN = O3
# Determines whether to find calibration filters in the old SVN or the git repo. Starting with PreER15 and O4, we switched to the git repo. Set to 'svn' or 'git'.
FILTERSRC = svn
-START = $(shell echo 1410782418 | bc)
+START = $(shell echo 1409576018 | bc)
# 1409576018 or 1410782418 gstlal-calibration-1.5.8 linesub issue
# 1399057926 line sub not working
# 1371862818 oversubtraction study
@@ -21,7 +21,7 @@ START = $(shell echo 1410782418 | bc)
# 1238288418
# 1269090018
# 1269114645-21449
-END = $(shell echo 1410839118 | bc)
+END = $(shell echo 1409604100 | bc)
# 1409604100 or 1410839118 gstlal-calibration-1.5.8 linesub issue
# 1399086062 line sub not working
# 1371911734 oversubtraction study
@@ -423,7 +423,7 @@ noise_subtraction_ASD_GDS: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_hoft_GDS_frames
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_hoft_GDS_frames.cache,$(IFO)1_hoft_GDS_frames_shortTF.cache,$(IFO)1_hoft_GDS_frames_longTF.cache,$(IFO)1_hoft_GDS_frames.cache --channel-list GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 512 --fft-spacing 256 --window=0 --freq-res 0.05 --frequency-min 1 --frequency-max=8192 --ASD-max=1e-17 --labels 'strain,128s\ estimate,4096s\ estimate,switch'
noise_subtraction_spectrogram_GDS: $(IFO)1_hoft_GDS_frames.cache
- python3 plot_spectrogram.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_GDS_frames.cache --channel-name GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 1 --fft-spacing 4 --window=3 --frequency-min 10 --frequency-max=1000 --update-time 16 --amplitude-min 0.01 --amplitude-max 100
+ python3 plot_spectrogram.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache $(IFO)1_hoft_GDS_frames.cache --channel-name GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --window=3 --frequency-min 10 --frequency-max=1000 --amplitude-min 0.01 --amplitude-max 100
noise_subtraction_ASD_GDS1: $(IFO)1_hoft_GDS_frames.cache $(IFO)1_hoft_GDS_frames_shortTF.cache $(IFO)1_hoft_GDS_frames_longTF.cache
python3 plot_ASD.py --gps-start-time $(PLOT_START) --gps-end-time $(PLOT_END) --ifo $(IFO)1 --frame-cache-list $(IFO)1_hoft_GDS_frames.cache,$(IFO)1_hoft_GDS_frames_shortTF.cache,$(IFO)1_hoft_GDS_frames_longTF.cache,$(IFO)1_hoft_GDS_frames.cache --channel-list GDS-CALIB_STRAIN,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES,GDS-CALIB_STRAIN_NOLINES --sample-rate 16384 --fft-time 4096 --fft-spacing 2048 --freq-res 0.001 --frequency-min 410.27 --frequency-max=410.33 --ASD-max=1e-19 --labels 'strain,128s\ estimate,4096s\ estimate,switch'
diff --git a/tests/check_calibration/plot_spectrogram.py b/tests/check_calibration/plot_spectrogram.py
index 13ff97c7e..3f084b9b7 100644
--- a/tests/check_calibration/plot_spectrogram.py
+++ b/tests/check_calibration/plot_spectrogram.py
@@ -67,11 +67,11 @@ parser.add_option("--ifo", metavar = "name", type = str, help = "Name of the int
parser.add_option("--frame-cache", metavar = "name", type = str, help = "Frame cache file that contains data from which to compute the spectrogram")
parser.add_option("--channel-name", metavar = "name", type = str, default = None, help = "Names of channel that contains data from which to compute the spectrogram")
parser.add_option("--sample-rate", metavar = "Hz", type = str, default = '16384', help = "Sample rate of input data.")
-parser.add_option("--fft-time", metavar = "seconds", type = float, default = 1, help = "Duration in seconds of FFTs used to compute spectrogram")
-parser.add_option("--fft-spacing", metavar = "seconds", type = float, default = 0.25, help = "Spacing in seconds of FFTs used to compute spectrogram")
+parser.add_option("--fft-time", metavar = "seconds", type = float, default = 4, help = "Duration in seconds of FFTs used to compute spectrogram")
+parser.add_option("--fft-spacing", metavar = "seconds", type = float, default = None, help = "Spacing in seconds of FFTs used to compute spectrogram. Default is auto.")
parser.add_option("--window", type = int, default = 3, help = "Which window function to use to window the data when taking FFTs. Options are 0 (DPSS), 1 (kaiser), 2 (Dolph-Chebyshev), 3 (Blackman, default), 4 (Hann), 5 (None)")
parser.add_option("--freq-res", metavar = "Hz", type = float, default = 1.0, help = "Frequency resolution of spectrogram. Only applies when using Kaiser, DPSS, and Dolph Chebyshev windows.")
-parser.add_option("--update-time", type = float, default = 4, help = "Amount of time between consecutive ASDs on the spectrogram")
+parser.add_option("--update-time", type = float, default = None, help = "Amount of time between consecutive ASDs on the spectrogram. Default is auto.")
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")
@@ -88,6 +88,21 @@ frame_cache = options.frame_cache
channel = options.channel_name
chan_list = [(ifo, channel)]
sr = sample_rate = int(options.sample_rate)
+update_time = options.update_time
+dur = options.gps_end_time - options.gps_start_time
+if update_time is None:
+ # ~1000 samples is about right
+ update_time = dur / 1000.0
+ # Make it a power of 2
+ update_time = 2**int(round(np.log2(update_time)))
+ # It should be more than the fft_time
+ while update_time < options.fft_time:
+ update_time *= 2
+fft_spacing = options.fft_spacing
+if fft_spacing is None:
+ fft_spacing = update_time / 8
+ while fft_spacing < options.fft_time / 4 and fft_spacing < update_time:
+ fft_spacing *= 2
#
# =============================================================================
@@ -112,9 +127,9 @@ def compute_spectrogram(pipeline, name):
data = calibration_parts.mkresample(pipeline, data, 5, False, sr)
data = pipeparts.mktee(pipeline, data)
# Compute ASDs for the spectrogram
- pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - options.fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = int(options.update_time), use_td_median = True, filename = '%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+ pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = update_time, use_td_median = True, filename = '%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# Compute one median ASD to normalize the above ASDs
- pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - options.fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = 2 * int(options.gps_end_time - options.gps_start_time), use_td_median = True, filename = '%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+ pipeparts.mkgeneric(pipeline, data, "lal_asd", fft_samples = int(options.fft_time * sr), overlap_samples = int((options.fft_time - fft_spacing) * sr), window_type = int(options.window), frequency_resolution = float(options.freq_res), update_time = 2 * int(dur), use_td_median = True, filename = '%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
#
# done
@@ -135,20 +150,20 @@ def compute_spectrogram(pipeline, name):
test_common.build_and_run(compute_spectrogram, "compute_spectrogram", segment = segments.segment((LIGOTimeGPS(0, 1000000000 * options.gps_start_time), LIGOTimeGPS(0, 1000000000 * options.gps_end_time))))
# Unnormalized spectrogram data
-spec = np.loadtxt('%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+spec = np.loadtxt('%s_%s_spectrogram_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# ASD to normalize it
-asd = np.loadtxt('%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, options.gps_end_time - options.gps_start_time))
+asd = np.loadtxt('%s_%s_ASD_%d-%d.txt' % (ifo, channel, options.gps_start_time, dur))
# Frequency vector
fvec = np.arange(0, sr / 2.0 + sr / 4.0 / len(asd), sr / 2.0 / (len(asd) - 1))
# Time vector
num_t = len(spec) // len(asd)
-max_t = options.update_time * (num_t - 1)
-tvec = np.arange(0, max_t + options.update_time / 2, options.update_time)
+max_t = update_time * (num_t - 1)
+tvec = np.arange(0, max_t + update_time / 2, update_time)
-t0 = options.gps_start_time + options.update_time / 2
+t0 = options.gps_start_time + update_time / 2
# Normalize spectrogram
asd = np.transpose(asd)[1]
@@ -189,7 +204,7 @@ plt.title(r'${\rm %s:%s \ Spectrogram}$' % (ifo, channel.replace('_', '\_')))
ticks_and_grid(ax, ymin = options.frequency_min, ymax = options.frequency_max, yscale = freq_scale, xscale = 'linear')
-plt.savefig('%s_%s_spectrogram_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, options.gps_end_time - options.gps_start_time))
+plt.savefig('%s_%s_spectrogram_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, dur))
plt.close()
# Median ASD for reference
@@ -201,6 +216,6 @@ 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.tight_layout()
-plt.savefig('%s_%s_ASD_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, options.gps_end_time - options.gps_start_time))
+plt.savefig('%s_%s_ASD_%d-%dHz_%d-%d.png' % (ifo, channel, int(options.frequency_min), int(options.frequency_max), options.gps_start_time, dur))
--
GitLab