Merge branch 'add-functionality-to-strain-data' into 'master'

Add functionality to strain data

See merge request Monash/tupak!94

.gitlab-ci.yml

@@ -31,6 +31,7 @@ exitcode-jessie:
     - coverage erase
     - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/conversion_tests.py
     - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/detector_tests.py
+    - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/utils_tests.py
     - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/prior_tests.py
     - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/sampler_tests.py
     - coverage run --source /usr/local/lib/python2.7/dist-packages/tupak/ -a test/waveform_generator_tests.py

requirements.txt

@@ -3,7 +3,7 @@ dynesty
 corner
 numpy>=1.9
 matplotlib>=2.0
-scipy
+scipy>=0.16
 gwpy
 pandas
 astropy

test/detector_tests.py

@@ -358,6 +358,23 @@ class TestInterferometerStrainData(unittest.TestCase):
             sampling_frequency=sampling_frequency, start_time=10)
         self.assertTrue(self.ifosd.start_time == 10)
 
+    def test_time_array_frequency_array_consistency(self):
+        duration = 1
+        sampling_frequency = 10
+        time_array = tupak.core.utils.create_time_series(
+            sampling_frequency=sampling_frequency, duration=duration)
+        time_domain_strain = np.random.normal(0, 1, len(time_array))
+        self.ifosd.roll_off = 0
+        self.ifosd.set_from_time_domain_strain(
+            time_domain_strain=time_domain_strain, duration=duration,
+            sampling_frequency=sampling_frequency)
+
+        frequency_domain_strain, freqs = tupak.core.utils.nfft(
+            time_domain_strain, sampling_frequency)
+
+        self.assertTrue(np.all(
+            frequency_domain_strain == self.ifosd.frequency_domain_strain))
+
     #def test_sampling_duration_init(self):
     #    self.assertIsNone(self.ifo.duration)
 

test/utils_tests.py

+from __future__ import absolute_import, division
+
+import tupak
+import unittest
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+class TestFFT(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def tearDown(self):
+        pass
+
+    def test_nfft_frequencies(self):
+        f = 2.1
+        sampling_frequency = 10
+        times = np.arange(0, 100, 1/sampling_frequency)
+        tds = np.sin(2*np.pi*times * f + 0.4)
+        fds, freqs = tupak.core.utils.nfft(tds, sampling_frequency)
+        self.assertTrue(np.abs((f-freqs[np.argmax(np.abs(fds))])/f < 1e-15))
+
+    def test_nfft_infft(self):
+        sampling_frequency = 10
+        tds = np.random.normal(0, 1, 10)
+        fds, _ = tupak.core.utils.nfft(tds, sampling_frequency)
+        tds2 = tupak.core.utils.infft(fds, sampling_frequency)
+        self.assertTrue(np.all(np.abs((tds - tds2) / tds) < 1e-12))
+
+
+if __name__ == '__main__':
+    unittest.main()

tupak/core/utils.py

@@ -33,6 +33,25 @@ def get_sampling_frequency(time_series):
         return 1. / (time_series[1] - time_series[0])
 
 
+def get_sampling_frequency_and_duration_from_time_array(time_array):
+    """
+    Calculate sampling frequency and duration from a time array
+
+    Returns
+    -------
+    sampling_frequency, duration:
+
+    Raises
+    -------
+    ValueError: If the time_array is not evenly sampled.
+
+    """
+
+    sampling_frequency = get_sampling_frequency(time_array)
+    duration = time_array[-1] - time_array[0]
+    return sampling_frequency, duration
+
+
 def get_sampling_frequency_and_duration_from_frequency_array(frequency_array):
     """
     Calculate sampling frequency and duration from a frequency array
@@ -201,59 +220,68 @@ def create_white_noise(sampling_frequency, duration):
     return white_noise, frequencies
 
 
-def nfft(ht, Fs):
-    """ Performs an FFT while keeping track of the frequency bins assumes input time series is real
-        (positive frequencies only)
+def nfft(time_domain_strain, sampling_frequency):
+    """ Perform an FFT while keeping track of the frequency bins. Assumes input
+        time series is real (positive frequencies only).
 
     Parameters
-    -------
-    ht: array_like
-        Time series
-    Fs: float
-        Sampling frequency
+    ----------
+    time_domain_strain: array_like
+        Time series of strain data.
+    sampling_frequency: float
+        Sampling frequency of the data.
 
     Returns
     -------
-    array_like: Single-sided FFT of ft normalised to units of strain / sqrt(Hz) (hf)
-    array_like: Frequencies associated with hf
+    frequency_domain_strain, frequency_array: (array, array)
+        Single-sided FFT of time domain strain normalised to units of
+        strain / Hz, and the associated frequency_array.
+
     """
-    # add one zero padding if time series does not have even number of sampling times
-    if np.mod(len(ht), 2) == 1:
-        ht = np.append(ht, 0)
-    LL = len(ht)
+
+    if np.ndim(sampling_frequency) != 0:
+        raise ValueError("Sampling frequency must be interger or float")
+
+    # add one zero padding if time series doesn't have even number of samples
+    if np.mod(len(time_domain_strain), 2) == 1:
+        time_domain_strain = np.append(time_domain_strain, 0)
+    LL = len(time_domain_strain)
     # frequency range
-    ff = Fs / 2 * np.linspace(0, 1, int(LL/2+1))
+    frequency_array = sampling_frequency / 2 * np.linspace(0, 1, int(LL/2+1))
 
     # calculate FFT
     # rfft computes the fft for real inputs
-    hf = np.fft.rfft(ht)
+    frequency_domain_strain = np.fft.rfft(time_domain_strain)
 
-    # normalise to units of strain / sqrt(Hz)
-    hf = hf / Fs
+    # normalise to units of strain / Hz
+    norm_frequency_domain_strain = frequency_domain_strain / sampling_frequency
 
-    return hf, ff
+    return norm_frequency_domain_strain, frequency_array
 
 
-def infft(hf, Fs):
-    """ Inverse FFT for use in conjunction with nfft (eric.thrane@ligo.org)
+def infft(frequency_domain_strain, sampling_frequency):
+    """ Inverse FFT for use in conjunction with nfft.
 
     Parameters
-    -------
-    hf: array_like
-        single-side FFT calculated by fft_eht
-    Fs: float
-        sampling frequency
+    ----------
+    frequency_domain_strain: array_like
+        Single-sided, normalised FFT of the time-domain strain data (in units
+        of strain / Hz).
+    sampling_frequency: float
+        Sampling frequency of the data.
 
     Returns
     -------
-    array_like: time series
+    time_domain_strain: array
+        An array of the time domain strain
     """
-    # use irfft to work with positive frequencies only
-    h = np.fft.irfft(hf)
-    # undo LAL/Lasky normalisation
-    h = h*Fs
 
-    return h
+    if np.ndim(sampling_frequency) != 0:
+        raise ValueError("Sampling frequency must be integer or float")
+
+    time_domain_strain_norm = np.fft.irfft(frequency_domain_strain)
+    time_domain_strain = time_domain_strain_norm * sampling_frequency
+    return time_domain_strain
 
 
 def setup_logger(outdir=None, label=None, log_level=None, print_version=False):

tupak/gw/utils.py

@@ -291,7 +291,7 @@ def get_event_time(event):
 
 
 def get_open_strain_data(
-        name, t1, t2, outdir, cache=False, raw_data_file=None, **kwargs):
+        name, t1, t2, outdir, cache=False, buffer_time=0, **kwargs):
     """ A function which accesses the open strain data
 
     This uses `gwpy` to download the open data and then saves a cached copy for
@@ -307,9 +307,10 @@ def get_open_strain_data(
         The output directory to place data in
     cache: bool
         If true, cache the data
+    buffer_time: float
+        Time to add to the begining and end of the segment.
     **kwargs:
         Passed to `gwpy.timeseries.TimeSeries.fetch_open_data`
-    raw_data_file
 
     Returns
     -----------
@@ -317,19 +318,18 @@ def get_open_strain_data(
 
     """
     filename = '{}/{}_{}_{}.txt'.format(outdir, name, t1, t2)
-    if raw_data_file is not None:
-        logging.info('Using raw_data_file {}'.format(raw_data_file))
-        strain = TimeSeries.read(raw_data_file)
-        if (t1 > strain.times[0].value) and (t2 < strain.times[-1].value):
-            logging.info('Using supplied raw data file')
-            strain = strain.crop(t1, t2)
-        else:
-            raise ValueError('Supplied file does not contain requested data')
-    elif os.path.isfile(filename) and cache:
+
+    if buffer_time < 0:
+        raise ValueError("buffer_time < 0")
+    t1 = t1 - buffer_time
+    t2 = t2 + buffer_time
+
+    if os.path.isfile(filename) and cache:
         logging.info('Using cached data from {}'.format(filename))
         strain = TimeSeries.read(filename)
     else:
-        logging.info('Fetching open data ...')
+        logging.info('Fetching open data from {} to {} with buffer time {}'
+                     .format(t1, t2, buffer_time))
         strain = TimeSeries.fetch_open_data(name, t1, t2, **kwargs)
         logging.info('Saving data to {}'.format(filename))
         strain.write(filename)