DEV: apply FD taper to waveform modes to remove Gibbs' features

gwmemory/utils.py

@@ -1,6 +1,7 @@
 from typing import Tuple
 
 import numpy as np
+from scipy.special import expit
 
 from .harmonics import lmax_modes, sYlm
 
@@ -97,3 +98,26 @@ def combine_modes(h_lm: dict, inc: float, phase: float) -> np.ndarray:
     total = sum([h_lm[(l, m)] * sYlm(-2, l, m, inc, phase) for l, m in h_lm])
     h_plus_cross = dict(plus=total.real, cross=-total.imag)
     return h_plus_cross
+
+
+def planck_taper(frequencies: np.ndarray, minimum: float, maximum: float):
+    """
+    A low-pass Planck taper for the provided frequencies.
+    This is intended to apply a frequency-domain high-pass filter.
+    This function takes the absolute value of the passed frequencies.
+
+    Parameters
+    ==========
+    frequencies: np.ndarray
+        Array containing the frequencies to calcualte the taper over.
+    minimum: float
+        The start point of the taper (the last frequency where the response
+        is 0)
+    maximum: float
+        The end point of the taper (the first frequency where the response
+        is 1)
+    """
+    normalized = (abs(frequencies) - minimum) / (maximum - minimum)
+    normalized = np.maximum(np.minimum(normalized, 1), 0)
+    with np.errstate(divide='ignore'):
+        return 2 * expit(1 - 1 / normalized)

gwmemory/waveforms/approximant.py

@@ -3,7 +3,7 @@
 import numpy as np
 
 from ..harmonics import sYlm
-from ..utils import MPC, SOLAR_MASS, combine_modes
+from ..utils import MPC, SOLAR_MASS, combine_modes, planck_taper
 from . import MemoryGenerator
 
 
@@ -76,6 +76,8 @@ def __init__(
         self.sampling_frequency = sampling_frequency
         self.duration = duration
         self.l_max = l_max
+        # The width of the taper for frequency-domain approximants [Hz]
+        self.taper_length = 5
 
         self.m1 = self.MTot / (1 + self.q)
         self.m2 = self.m1 * self.q
@@ -166,7 +168,7 @@ def time_domain_oscillatory(
         if self.h_lm is None:
 
             params = dict(
-                f_min=self.minimum_frequency,
+                f_min=self.minimum_frequency - self.taper_length,
                 f_ref=self.reference_frequency,
                 phiRef=0.0,
                 approximant=lalsim.GetApproximantFromString(self.name),
@@ -208,14 +210,28 @@ def time_domain_oscillatory(
                 params["f_max"] = self.sampling_frequency / 2
                 waveform_modes = lalsim.SimInspiralChooseFDModes(**params)
                 times = np.arange(0, duration, 1 / self.sampling_frequency)
+                n_frequencies = int(duration * params["f_max"])
+                frequencies = waveform_modes.fdata.data[:-1]
+                # We apply a frequency-domain window to avoid introducing Gibbs
+                # artefacts.
+                window = planck_taper(
+                    frequencies,
+                    self.minimum_frequency - self.taper_length,
+                    self.minimum_frequency,
+                )
 
                 h_lm = dict()
                 while waveform_modes is not None:
                     mode = (waveform_modes.l, waveform_modes.m)
-                    data = waveform_modes.mode.data.data[:-1]
+                    data = waveform_modes.mode.data.data[:-1] * window
+                    # The LAL frequency series uses [negative, 0, positive]
+                    # but numpy expects [0, positive, negative], we roll the
+                    # LAL waveform to account for this.
+                    # Finally, we roll the time domain waveform to place the
+                    # merger 1s before the end.
                     h_lm[mode] = (
                         np.roll(
-                            np.fft.ifft(np.roll(data, int(duration * params["f_max"]))),
+                            np.fft.ifft(np.roll(data, n_frequencies)),
                             int((duration - 1) * self.sampling_frequency),
                         )
                         * self.sampling_frequency