diff --git a/tupak/utils.py b/tupak/utils.py index 6b6f575ae53843ff279d8d6a076d224b1c4b467b..381265541233cd01724481dc4d8b2070ee97f70b 100644 --- a/tupak/utils.py +++ b/tupak/utils.py @@ -4,6 +4,8 @@ import os import numpy as np from math import fmod from gwpy.timeseries import TimeSeries +from gwpy.signal import filter_design +from scipy import signal import argparse # Constants @@ -514,6 +516,110 @@ def get_open_strain_data( return strain +def read_frame_file(file_name, t1, t2, channel=None, **kwargs): + """ A function which accesses the open strain data + + This uses `gwpy` to download the open data and then saves a cached copy for + later use + + Parameters + ---------- + file_name: str + The name of the frame to read + t1, t2: float + The GPS time of the start and end of the data + channel: str + The name of the channel being searched for, some standard channel names are attempted + if channel is not specified or if specified channel is not found. + **kwargs: + Passed to `gwpy.timeseries.TimeSeries.fetch_open_data` + + Returns + ----------- + strain: gwpy.timeseries.TimeSeries + + """ + loaded = False + if channel is not None: + try: + strain = TimeSeries.read(source=file_name, channel=channel, start=t1, end=t2, **kwargs) + loaded = True + logging.info('Successfully loaded {}.'.format(channel)) + except RuntimeError: + logging.warning('Channel {} not found. Trying preset channel names'.format(channel)) + for channel in ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02']: + if loaded: + continue + try: + strain = TimeSeries.read(source=file_name, channel=channel, start=t1, end=t2, **kwargs) + loaded = True + logging.info('Successfully loaded {}.'.format(channel)) + except RuntimeError: + None + + if loaded: + return strain + else: + logging.warning('No data loaded.') + return None + + +def process_strain_data( + strain, alpha=0.25, filter_freq=1024, **kwargs): + """ + Helper function to obtain an Interferometer instance with appropriate + PSD and data, given an center_time. + + Parameters + ---------- + name: str + Detector name, e.g., 'H1'. + center_time: float + GPS time of the center_time about which to perform the analysis. + Note: the analysis data is from `center_time-T/2` to `center_time+T/2`. + T: float + The total time (in seconds) to analyse. Defaults to 4s. + alpha: float + The tukey window shape parameter passed to `scipy.signal.tukey`. + psd_offset, psd_duration: float + The power spectral density (psd) is estimated using data from + `center_time+psd_offset` to `center_time+psd_offset + psd_duration`. + outdir: str + Directory where the psd files are saved + plot: bool + If true, create an ASD + strain plot + filter_freq: float + Low pass filter frequency + **kwargs: + All keyword arguments are passed to + `gwpy.timeseries.TimeSeries.fetch_open_data()`. + + Returns + ------- + interferometer: `tupak.detector.Interferometer` + An Interferometer instance with a PSD and frequency-domain strain data. + + """ + + sampling_frequency = int(strain.sample_rate.value) + + # Low pass filter + bp = filter_design.lowpass(filter_freq, strain.sample_rate) + strain = strain.filter(bp, filtfilt=True) + strain = strain.crop(*strain.span.contract(1)) + + time_series = strain.times.value + time_duration = time_series[-1] - time_series[0] + + # Apply Tukey window + N = len(time_series) + strain = strain * signal.windows.tukey(N, alpha=alpha) + + frequency_domain_strain = nfft(strain.value, sampling_frequency)[0] + + return frequency_domain_strain + + def set_up_command_line_arguments(): parser = argparse.ArgumentParser( description="Command line interface for tupak scripts")