From 1b2022873139700b4583e46fd8f75a2a7bd1b15d Mon Sep 17 00:00:00 2001
From: Colm Talbot <talbotcolm@gmail.com>
Date: Sun, 23 Jul 2023 23:13:01 -0400
Subject: [PATCH 1/2] DEV: apply FD taper to waveform modes to remove Gibbs'
 features

---
 gwmemory/utils.py                 | 23 +++++++++++++++++++++++
 gwmemory/waveforms/approximant.py | 24 ++++++++++++++++++++----
 2 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/gwmemory/utils.py b/gwmemory/utils.py
index 9ba3d57..663395b 100644
--- a/gwmemory/utils.py
+++ b/gwmemory/utils.py
@@ -97,3 +97,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)
diff --git a/gwmemory/waveforms/approximant.py b/gwmemory/waveforms/approximant.py
index e61958c..8a7c16a 100644
--- a/gwmemory/waveforms/approximant.py
+++ b/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

From 919e2192f2ed1f2e01c4875ebd51035178e1010c Mon Sep 17 00:00:00 2001
From: Colm Talbot <talbotcolm@gmail.com>
Date: Sun, 23 Jul 2023 23:36:16 -0400
Subject: [PATCH 2/2] BUGFIX: add missing import

---
 gwmemory/utils.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/gwmemory/utils.py b/gwmemory/utils.py
index 663395b..649327f 100644
--- a/gwmemory/utils.py
+++ b/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