Clean up documentation and notation of nfft and infft and add tests

.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

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

@@ -220,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 / sqrt(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
+    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 / sqrt(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 interger 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):