# Copyright © 2021-2025 HQS Quantum Simulations GmbH. All Rights Reserved.
import numpy as np
[docs]
def fft_causal(f: np.ndarray) -> np.ndarray:
"""Performs the Fourier Transform of a causal signal separating its odd and even part.
This way of transforming is necessary for causal signals or any signal that does not start and
end at the same value. This is because the fft replicates the signal to be able to calculate the
coefficients, leading to sudden jumps if the function itself is not periodic. These jumps
introduce a shift in the FT at all frequencies, being equivalent to a delta within the sampled
time frame. By (anti-)symmetrizing the signal one can circumvent this problem and obtain the
correct Fourier Transform.
Args:
f (np.ndarray): Function to transform.
Returns:
np.ndarray: The complex Fourier Transform.
"""
# total length of the signal
L = len(f)
# -----------------------
# even part of the signal
# -----------------------
# Given a function f(t), its even part fe(t) = 1/2 (f(t) + f(-t))
# Assuming the singal f(t) is a causal function of t in [0, Tmax], its even part has a
# domain [-Tmax, Tmax]. The value at t=0 is not repeated to make the function perfectly even,
# hence the total length is 2L-1.
# Since the fft replicates the signal, if fe is "centered" around t=0, its replicated version
# would not be the correct function to transform, as it would be shifted by Tmax, introducing a
# factor e^(-i Tmax omega). To avoid this, we "shift" back the signal by Tmax, meaning
# fe = f for the first half, and then is f(-t). Equivalently, one could multiply the FT by the
# factor e^(i Tmax omega) to cancel the one introduced by the shift in time.
fe = np.zeros(2 * L - 1, dtype=f.dtype)
fe[:L] = f
fe[L:] = f[:-1][::-1]
fe /= 2
# -----------------------
# odd part of the signal
# -----------------------
# the odd part is defined as fo(t) = 1/2 ( f(t) - f(-t))
# Everything is consistent with what done for the even part, except for the value at t=0, which,
# in this case is set to 0. This is, again, to avoid sudden jumps of the signal around 0. In
# fact, the odd function will be, for t=[-dt, 0, dt], -f(1), 0, f(1), avoiding delta-like jumps.
# NOTE: the same reasoning can be applied to the last point of f(t) if this is not ~0. Although
# this does not make a large difference, in general. Moreover we are mainly interested in
# decaying signals, making this extra step redundant.
fo = np.zeros(2 * L - 1, dtype=f.dtype)
fo[:L] = f
fo[L:] = -f[:-1][::-1]
fo /= 2
fo[0] = 0
# fo[L] = 0
# Since we have effectively doubled the time window in which the signal is defined, the
# resolution in frequency will be finer. This is only an artefact of the construction. We
# therefore only select every second value for the FT. This results in exactly L values for the
# FT at exactly the same frequencies we would have had without the even/odd separation.
# the normalization is left up to the user as it may vary from case to case
return (
np.fft.fftshift(np.fft.fft(fe, norm="ortho"))
+ np.fft.fftshift(np.fft.fft(fo, norm="ortho"))
)[::2]