plot_spectrogram.py: better default behavior for plot parameters

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'

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))