To get the most out of this paper, you should have a basic understanding of integrals and the Fourier transform and its properties.
If not, we recommend reading Section 1, where these ideas are reviewed.
If you are already familiar with them, you can go directly to Section 2.

In Section 2.1, we implement the Fourier transform using a numerical quadrature method.
In Section 2.2, we implement it using the Fast Fourier Transform (FFT) algorithm and Shannon’s interpolation formula.

this paper came to us during our final year at engineering school.

As part of a course on stochastic process calibration, our professor wanted to test our understanding of the material. To do so, we were asked to select an academic paper, study it in detail, and reproduce its results

We chose the paper by El Kolei (2013), which proposes a parametric method for estimating the parameters of a stochastic volatility model, such as the GARCH model. The approach is based on contrast minimization and deconvolution.

To reproduce the results, we needed to implement an optimization method that involved computing the Fourier transform f̂ of the function f_θ:

To compute the Fourier transform of f_θ, El Kolei (2013) used the left Riemann sum (also known as the rectangle quadrature) method implemented in MATLAB. The paper suggested that using the Fast Fourier Transform (FFT) algorithm could make the computation faster.

We decided to reproduce the results using Python and implemented the Fourier transform to test if it truly improved computation speed.. Implementing the left rectangle quadrature method was simple, and at first, we thought that the scipy.fft and numpy.fft libraries would allow us to compute the Fourier transform of f_θ directly using the FFT algorithm.

However, we quickly discovered that these functions do not compute the Fourier transform of a continuous function. Instead, they compute the Discrete Fourier Transform (DFT) of a finite sequence.

The figure below shows what we observed.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, fftshift

# Define the function f(t) = exp(-pi * t^2)
def f(t):
    return np.exp(-np.pi * t**2)

# Parameters
N = 1024
T = 1.0 / 64
t = np.linspace(-N/2*T, N/2*T, N, endpoint=False)
y = f(t)

# FFT with scipy
yf_scipy = fftshift(fft(y)) * T
xf = fftshift(fftfreq(N, T))
FT_exact = np.exp(-np.pi * xf**2)

# FFT with numpy
yf_numpy = np.fft.fftshift(np.fft.fft(y)) * T
xf_numpy = np.fft.fftshift(np.fft.fftfreq(N, T))

# Plot with subplot_mosaic
fig, axs = plt.subplot_mosaic([["scipy", "numpy"]], figsize=(7, 5), layout="constrained", sharey=True)

# Scipy FFT
axs["scipy"].plot(xf, FT_exact, 'k-', linewidth=1.5, label='Exact FT')
axs["scipy"].plot(xf, np.real(yf_scipy), 'r--', linewidth=1, label='FFT (scipy)')
axs["scipy"].set_xlim(-6, 6)
axs["scipy"].set_ylim(-1, 1)
axs["scipy"].set_title("Scipy FFT")
axs["scipy"].set_xlabel("Frequency")
axs["scipy"].set_ylabel("Amplitude")
axs["scipy"].legend()
axs["scipy"].grid(False)

# NumPy FFT
axs["numpy"].plot(xf_numpy, FT_exact, 'k-', linewidth=1.5, label='Exact FT')
axs["numpy"].plot(xf_numpy, np.real(yf_numpy), 'b--', linewidth=1, label='FFT (numpy)')
axs["numpy"].set_xlim(-6, 6)
axs["numpy"].set_title("NumPy FFT")
axs["numpy"].set_xlabel("Frequency")
axs["numpy"].legend()
axs["numpy"].grid(False)

plt.suptitle("Comparison of FFT Implementations vs. Exact Fourier Transform", fontsize=14)
plt.show()

This experience inspired us to write this article, where we will explain how to compute the Fourier transform of a function in Python using two approaches: the left Riemann sum method and the Fast Fourier Transform (FFT) algorithm.

In the literature, many papers discuss how to approximate the Fourier transform and how to implement it using numerical methods.

However, we did not find a source as clear and complete as Balac (2011). This work presents a quadrature-based approach to compute the Fourier transform and explains how to use the Fast Fourier Transform (FFT) algorithm to perform the discrete Fourier transform efficiently.

The FFT algorithm can be applied both to compute the Fourier transform of an integrable function and to calculate the Fourier coefficients of a periodic function.

1. Definition and Properties of the Fourier Transform

We adopt the framework presented by Balac (2011) to define the Fourier transform of a function f and to introduce its properties.

We consider functions that belong to the space of integrable functions, denoted L¹(ℝ, 𝕂). This space includes all functions f: ℝ → 𝕂, where 𝕂 represents either the real numbers (ℝ) or the complex numbers (ℂ).

These functions are integrable in the sense of Lebesgue, meaning that the integral of their absolute value is finite.

So, for a function f to belong to L¹(ℝ, 𝕂), the product f(t) · e^–2iπνt must also be integrable for all ν ∈ ℝ.

In that case, the Fourier transform of f, denoted f̂ (or sometimes 𝓕(f)), is defined for all ν ∈ ℝ by:

We note that the Fourier transform f̂ of a function f is a complex-valued, linear function that depends on the frequency ν.

If f ∈ L¹(ℝ), is real-valued and even, then f̂ is also real-valued and even. Conversely, if f is real-valued and odd, then f̂ is purely imaginary and odd as well.

For some functions, the Fourier transform can be computed analytically. For example, for the function
f : t ∈ ℝ ↦ 𝟙 _{[−a⁄2, a⁄2]}(t),
The Fourier transform is given by:

f̂(ν) = a sinc(a π ν)

where  is the sinc(t) = sin(t)/t or t ∈ ℝ^*, and sinc(0) = 0.

However, for many functions, the Fourier transform cannot be computed analytically. In such cases, we use numerical methods to approximate it. We explore these numerical approaches in the following sections of this article.

2. How to numerically approximate the Fourier Transform?

In his article, Balac (2011) shows that computing the Fourier transform of a function involves approximating it by the following integral over the interval [−T⁄2, T⁄2]:

Where T is a sufficiently large number for the integral to converge.. An approximate value of the integral S̃(ν) can be computed using quadrature methods.

In the next section, we will approximate this integral using the left rectangle quadrature method (also known as the left Riemann sum).

2.1 Quadrature method of left Rectangles

To compute the integral S̃(ν) using the left rectangle quadrature method, we proceed as follows:

1. Discretization of the Interval
   We divide the interval [−T⁄2, T⁄2] into N uniform subintervals of length h_t = T⁄N.
   The discretization points, corresponding to the left endpoints of the rectangles, are defined as:

t_k = −T⁄2 + h_t, for k = 0, 1, …, N−1.

2. Approximation of the Integral
   Using the Chasles relation, we approximate the integral S̃(ν) as follows:

Since tₖ₊₁ − tₖ = hₜ and tₖ = −T⁄2 + k·hₜ = T(k⁄N − ½), the expression becomes:

We refer to this approach as the left rectangle quadrature method because it utilizes the left endpoint tₖ of each subinterval to approximate the value of f(t) within that interval.

3. Final Formula
   The final expression for approximating the Fourier transform is given by:

2.1.1 Implementation of the left rectangle quadrature method in Python.

The function tfquad below implements the left rectangle quadrature method to compute the Fourier transform of a function f at a given frequency nu.

import numpy as np

def tfquad(f, nu, n, T):
    """
    Computes the Fourier transform of a function f at frequency nu
    using left Riemann sum quadrature over the interval [-T/2, T/2].

    Parameters:
    ----------
    f : callable
        The function to transform. Must accept a NumPy array as input.
    nu : float
        The frequency at which to evaluate the Fourier transform.
    n : int
        Number of quadrature points.
    T : float
        Width of the time window [-T/2, T/2].

    Returns:
    -------
    tfnu : complex
        Approximated value of the Fourier transform at frequency nu.
    """
    k = np.arange(n)
    t_k = (k / n - 0.5) * T
    weights = np.exp(-2j * np.pi * nu * T * k / n)
    prefactor = (T / n) * np.exp(1j * np.pi * nu * T)

    return prefactor * np.sum(f(t_k) * weights)

We can also use SciPy’s quad function to define the Fourier transform of a function f at a given frequency nu. The function tf_integral below implements this approach. It uses numerical integration to compute the Fourier transform of f over the interval [-T/2, T/2].

from scipy.integrate import quad

def tf_integral(f, nu, T):
    """Compute FT of f at frequency nu over [-T/2, T/2] using scipy quad."""
    real_part = quad(lambda t: np.real(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
    imag_part = quad(lambda t: np.imag(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
    return real_part + 1j * imag_part

After deriving the formula for the left rectangle quadrature method, we can now implement it in Python to see how it performs in practice.

To do this, we consider a simple example where the function fff is defined as an indicator function on the interval [−1, 1]:

For this function, the Fourier transform can be computed analytically, which allows us to compare the numerical approximation with the exact result.

The following Python script implements both the analytical Fourier transform and its numerical approximations using:

• the left rectangle quadrature method, and
• SciPy’s integration function for reference.

import numpy as np
import matplotlib.pyplot as plt

# ----- Function Definitions -----

def f(t):
    """Indicator function on [-1, 1]."""
    return np.where(np.abs(t) <= 1, 1.0, 0.0)

def exact_fourier_transform(nu):
    """Analytical FT of the indicator function over [-1, 1]."""
    # f̂(ν) = ∫_{-1}^{1} e^{-2πiνt} dt = 2 * sinc(2ν)
    return 2 * np.sinc(2 * nu)


# ----- Computation -----

T = 2.0
n = 32
nu_vals = np.linspace(-6, 6, 500)
exact_vals = exact_fourier_transform(nu_vals)
tfquad_vals = np.array([tfquad(f, nu, n, T) for nu in nu_vals])

# Compute the approximation using scipy integral
tf_integral_vals = np.array([tf_integral(f, nu, T) for nu in nu_vals])

# ----- Plotting -----
fig, axs = plt.subplot_mosaic([["tfquad", "quad"]], figsize=(7.24, 4.07), dpi=100, layout="constrained")

# Plot using tfquad implementation
axs["tfquad"].plot(nu_vals, np.real(exact_vals), 'b-', linewidth=2, label=r'$\hat{f}$ (exact)')
axs["tfquad"].plot(nu_vals, np.real(tfquad_vals), 'r--', linewidth=1.5, label=r'approximation $\hat{S}_n$')
axs["tfquad"].set_title("TF avec tfquad (rectangles)")
axs["tfquad"].set_xlabel(r'$\nu$')
axs["tfquad"].grid(False)
axs["tfquad"].set_ylim(-0.5, 2.1)

# Plot using scipy.integrate.quad
axs["quad"].plot(nu_vals, np.real(exact_vals), 'b-', linewidth=2, label=r'$\hat{f}$ (quad)')
axs["quad"].plot(nu_vals, np.real(tf_integral_vals), 'r--', linewidth=1.5, label=r'approximation $\hat{S}_n$')
axs["quad"].set_title("TF avec scipy.integrate.quad")
axs["quad"].set_xlabel(r'$\nu$')
axs["quad"].set_ylabel('Amplitude')
axs["quad"].grid(False)
axs["quad"].set_ylim(-0.5, 2.1)



# --- Global legend below the plots ---
# Take handles from one subplot only (assumes labels are consistent)
handles, labels = axs["quad"].get_legend_handles_labels()
fig.legend(handles, labels,
           loc='lower center', bbox_to_anchor=(0.5, -0.05),
           ncol=3, frameon=False)

plt.suptitle("Comparison of FFT Implementations vs. Exact Fourier Transform", fontsize=14)

plt.show()

2.1.2 Characterization of approximation using the left rectangle quadrature method

• The approximation of the Fourier transform  using the left rectangle quadrature method exhibits an oscillatory behavior by nature.

  As noted by Balac (2011), the Fourier transform f̂ of a function f is inherently oscillatory. This behavior comes from the complex exponential term e⁻²πiνt in the integral definition of the transform.

To illustrate this behavior, the figure below shows the function

f : t ∈ ℝ ↦ e^{-t^2} ∈ ℝ

together with the real and imaginary parts of its Fourier transform

f̂ : ν ∈ ℝ ↦ f̂(ν) ∈ ℂ, evaluated at ν = 5⁄2.

Although f(t) is smooth, we can observe strong oscillations in f̂(ν), which reveal the influence of the exponential term e^-2πiνt in the Fourier integral.

import numpy as np
import matplotlib.pyplot as plt

nu = 5 / 2
t1 = np.linspace(-8, 8, 1000)
t2 = np.linspace(-4, 4, 1000)

f = lambda t: np.exp(-t**2)
phi = lambda t: f(t) * np.exp(-2j * np.pi * nu * t)

f_vals = f(t1)
phi_vals = phi(t2)

# Plot
fig, axs = plt.subplots(1, 2, figsize=(10, 4))

axs[0].plot(t1, f_vals, 'k', linewidth=2)
axs[0].set_xlim(-8, 8)
axs[0].set_ylim(0, 1)
axs[0].set_title(r"$f(t) = e^{-t^2}$")
axs[0].grid(True)

axs[1].plot(t2, np.real(phi_vals), 'b', label=r"$\Re(\phi)$", linewidth=2)
axs[1].plot(t2, np.imag(phi_vals), 'r', label=r"$\Im(\phi)$", linewidth=2)
axs[1].set_xlim(-4, 4)
axs[1].set_ylim(-1, 1)
axs[1].set_title(r"$\phi(t) = f(t)e^{-2i\pi\nu t}$, $\nu=5/2$")
axs[1].legend()
axs[1].grid(True)

plt.tight_layout()
plt.show()

• The approximation obtained using the left rectangle quadrature method is periodic in nature.

  We observe that even if the function f is not periodic, its Fourier transform approximation f̂ appears to be periodic.

  In fact, the function Ŝₙ, obtained from the quadrature method, is periodic with a period given by:

This periodicity of Ŝₙ implies that it is not possible to compute the Fourier transform for all frequencies ν ∈ ℝ using the quadrature method when the parameters T and n are fixed.

In fact, it becomes impossible to compute f̂(ν) accurately when

|ν| ≥ ν_max,

where v = n/T represents the maximum frequency that can be resolved because of the periodic nature of Ŝₙ.

As a result, in practice, to compute the Fourier transform for larger frequencies, one must increase either the time window (T) or the number of points (n).

Furthermore, by evaluating the error in approximating f̂(ν) using the left rectangle quadrature method, we can show that the approximation is reliable at frequency ν when the following condition holds:

|ν| ≪ n/T or equivalently, (|ν|T)/n ≪ 1.

According to Epstein (2005), when using the Fast Fourier Transform (FFT) algorithm, it is possible to accurately compute the Fourier transform of a function f for all frequencies ν ∈ ℝ, even when (|ν|T)⁄n is close to 1, provided that f is piecewise continuous and has compact support.

2.2 Computing the Fourier Transform at Frequency ν Using the FFT Algorithm

In this section, we denote by Ŝₙ(ν) the approximation of the Fourier transform f̂(ν) of the function f at a point ν ∈ [−ν_max/2, ν_max/2], where ν_max = n/T, that is,

I now present the Fourier transform algorithm used to approximate f̂(ν).
I will not go into the details of the Fast Fourier Transform (FFT) algorithm in this article.
For a simplified explanation, I refer to Balac (2011), and for a more technical and comprehensive treatment, to the original work by Cooley and Tukey (1965).

It is important to understand that the use of the FFT algorithm to approximate the Fourier transform of a function f is based on the result established by Epstein (2005). This result states that when Ŝₙ(ν) is evaluated at the frequencies v_j = j/T, for j = 0, 1, …, n − 1, it provides a good approximation of the continuous Fourier transform f̂(ν).

Moreover, Ŝₙ is known to be periodic. This periodicity assigns a symmetric role to the indices j ∈ {0, 1, …, n − 1} and k ∈ {−n/2, −n/2 + 1, …, −1}.

In fact, the values of the Fourier transform of f over the interval [−ν_max/2, ν_max/2] can be derived from the values of Ŝₙ at the points ν_j = j/T, for j = 0, 1, …, n − 1, as follows:

where we used the relation:

This relationship shows that the Fourier transform Ŝₙ(j/T) can be computed for j = −n⁄2, …, n⁄2 − 1.

Moreover, when n is a power of two, the computation becomes significantly faster (see Balac, 2011). This process is known as the Fast Fourier Transform (FFT).

To summarize, I have shown that the Fourier transform of the function f can be approximated over the interval [−T⁄2, T⁄2] at the frequencies v_j = j/T, for j = −n/2, …, n/2 − 1, where n = 2ᵐ for some integer m ≥ 0, by applying the Fast Fourier Transform (FFT) algorithm as follows:

• Step 1: Construct the finite sequence F of values
  f((2k − n)T/(2n)), for k = 0, 1, …, n − 1.
• Step 2: Compute the discrete Fourier transform (DFT) of the sequence F using the FFT algorithm, given by:

• Step 3: Reindex and symmetrize the values to span j = −n/2, …, −1.
• Step 4: Multiply each value in the array by (T/n)(−1)^j-1, where j ∈ {1, …, n}.

This process yields an array representing the values of the Fourier transform f̂(νⱼ), with ν_j = j/T for j = −n⁄2, …, n⁄2 − 1.

The Python function tffft below implements these steps to compute the Fourier transform of a given function.

import numpy as np
from scipy.fft import fft, fftshift

def tffft(f, T, n):
    """
    Calcule la transformée de Fourier approchée d'une fonction f à support dans [-T/2, T/2],
    en utilisant l’algorithme FFT.

    Paramètres
    ----------
    f : callable
        Fonction à transformer (doit être vectorisable avec numpy).
    T : float
        Largeur de la fenêtre temporelle (intervalle [-T/2, T/2]).
    n : int
        Nombre de points de discrétisation (doit être une puissance de 2 pour FFT efficace).

    Retours
    -------
    tf : np.ndarray
        Valeurs approximées de la transformée de Fourier aux fréquences discrètes.
    freq_nu : np.ndarray
        Fréquences discrètes correspondantes (de -n/(2T) à (n/2 - 1)/T).
    """
    h = T / n
    t = -0.5 * T + np.arange(n) * h  # noeuds temporels
    F = f(t)                         # échantillonnage de f
    tf = h * (-1) ** np.arange(n) * fftshift(fft(F))  # TF approximée
    freq_nu = -n / (2 * T) + np.arange(n) / T              # fréquences ν_j = j/T

    return tf, freq_nu, t

The following program illustrates how to compute the Fourier transform of the Gaussian function f(t) = e^{-10t^2} over the interval [-10, 10], using the Fast Fourier Transform (FFT) algorithm.

# Parameters
a = 10
f = lambda t: np.exp(-a * t**2)
T = 10
n = 2**8  # 256

# Compute the Fourier transform using FFT
tf, nu, t = tffft(f, T, n)

# Plotting
fig, axs = plt.subplots(1, 2, figsize=(7.24, 4.07), dpi=100)

axs[0].plot(t, f(t), '-g', linewidth=3)
axs[0].set_xlabel("time")
axs[0].set_title("Considered Function")
axs[0].set_xlim(-6, 6)
axs[0].set_ylim(-0.5, 1.1)
axs[0].grid(True)

axs[1].plot(nu, np.abs(tf), '-b', linewidth=3)
axs[1].set_xlabel("frequency")
axs[1].set_title("Fourier Transform using FFT")
axs[1].set_xlim(-15, 15)
axs[1].set_ylim(-0.5, 1)
axs[1].grid(True)

plt.tight_layout()
plt.show()

The method I just presented makes it possible to compute and visualize the Fourier transform of a function f at discrete points v_j = j/T, for j = −n/2, …, n/2 − 1, where n is a power of two.

These points lie within the interval [−n/(2T), n/(2T)].
However, this method does not allow us to evaluate the Fourier transform at arbitrary frequencies ν that are not of the form v_j = j/T.

To compute the Fourier transform of a function f at a frequency ν that does not match one of the sampled frequencies νⱼ, interpolation methods can be used, such as linear, polynomial, or spline interpolation.

In this article, we use the approach proposed by Balac (2011), which relies on Shannon’s interpolation theorem to compute the Fourier transform of f at any frequency ν.

Using Shannon’s Interpolation Theorem to Compute the Fourier Transform of a Function f at a Point ν

What does Shannon’s theorem tell us?

It states that for a band-limited function g — that is, a function whose Fourier transform ĝ is zero outside the interval [−B⁄2, B⁄2] — the function g can be reconstructed from its samples gₖ = g(k⁄B) for k ∈ ℤ.

If we let ν_c be the smallest positive number such that ĝ(ν) is zero outside the interval [−2πν_c, 2πν_c], then the Shannon interpolation formula applies.
For every t ∈ ℝ, and any positive real number α ≥ 1/(2ν_c), we have:

where sinc(x) = sin(x)/x is the sinc function (cardinal sine).

Balac (2011) shows that when a function f has a bounded support in the interval [−T⁄2, T⁄2], Shannon’s interpolation theorem can be used to compute its Fourier transform f̂(ν) for any ν ∈ ℝ.
This is done by using the values of the discrete Fourier transform Ŝₙ(νⱼ) for j = −n/2, …, (n/2 − 1).

By setting α = 1/T, Balac derives the following Shannon interpolation formula, valid for all ν ∈ ℝ:

The following program illustrates how to use Shannon’s interpolation theorem to compute the Fourier transform of a function f at a given frequency ν.

To compute the Fourier transform at any frequency ν, we can use Shannon’s interpolation theorem.

The idea is to estimate the value of the Fourier transform from its discrete samples obtained by the FFT.

The function below, shannon, implements this interpolation method in Python.

def shannon(tf, nu, T):
    """
    Approximates the value of the Fourier transform of function f at frequency 'nu'
    using its discrete values computed from the FFT.

    Parameters:
    - tf : numpy array, discrete Fourier transform values (centered with fftshift) at frequencies j/T for j = -n/2, ..., n/2 - 1
    - nu : float, frequency at which to approximate the Fourier transform
    - T  : float, time window width used for the FFT

    Returns:
    - tfnu : approximation of the Fourier transform at frequency 'nu'
    """
    n = len(tf)
    tfnu = 0.0
    for j in range(n):
        k = j - n // 2  # correspond à l'indice j dans {-n/2, ..., n/2 - 1}
        tfnu += tf[j] * np.sinc(T * nu - k)  # np.sinc(x) = sin(pi x)/(pi x) en numpy

    return tfnu

The next function, fourier_at_nu, combines two steps.
It first computes the discrete Fourier transform of a function using the FFT (tffft), then applies Shannon interpolation to estimate the Fourier transform at a specific frequency ν.

This makes it possible to evaluate the transform at any arbitrary point ν, not just at the FFT sampling frequencies.

We can now test our implementation using a simple example.
Here, we define a function f(t) = e^-a|t|, whose Fourier transform is known analytically.

This lets us compare the exact Fourier transform with the approximation obtained using our method.

a = 0.5
f = lambda t: np.exp(-a * np.abs(t))                          # Function to transform
fhat_exact = lambda nu: (2 * a) / (a**2 + 4 * np.pi**2 * nu**2)  # Exact Fourier transform

T = 40     # Time window
n = 2**10  # Number of discretization points

We evaluate the Fourier transform at two frequencies: ν = 3/T and ν = π/T.

For each case, we print both the exact value and the approximation obtained by our fourier_at_nu function.

This comparison helps verify the accuracy of Shannon interpolation.

# Compute for nu = 3/T

nu = 3 / T

exact_value = fhat_exact(nu)
approx_value = fourier_at_nu(f, T, n, nu)
print(f"Exact value at nu={nu}: {exact_value}")
print(f"Approximation at nu={nu}: {np.real(approx_value)}")

# Compute for nu = pi/T

nu = np.pi / T

exact_value = fhat_exact(nu)
approx_value = fourier_at_nu(f, T, n, nu)
print(f"Exact value at nu={nu}: {exact_value}")
print(f"Approximation at nu={nu}: {np.real(approx_value)}")

Exact value at ν = 3/T: 2.118347413776218  
Approximation at ν = 3/T: (2.1185707502943534)

Exact value at ν = π/T: 2.0262491352594427  
Approximation at ν = π/T: (2.0264201680784835)

It is worth noting that even though the frequency 3⁄T (whose Fourier transform value appears in the FFT output) is close to π⁄T, the corresponding values of the Fourier transform of f at these two frequencies are quite different.

This shows that Shannon’s interpolation formula can be very useful when the desired Fourier transform values are not directly included among the results produced by the FFT algorithm.

Conclusion

In this article, we explored two methods to approximate the Fourier transform of a function.

The first was a numerical quadrature method (the left rectangle rule), and the second was the Fast Fourier Transform (FFT) algorithm.

We showed that, unlike the built-in fft functions in SciPy or NumPy, the FFT can be adapted to compute the Fourier transform of a continuous function through proper formulation.

Our implementation was based on the work of Balac (2013), who demonstrated how to reproduce the FFT computation in Python.
We also introduced Shannon’s interpolation theorem, which allows us to estimate the Fourier transform of any function on ℝ at arbitrary frequencies.

This interpolation can be replaced by more traditional approaches such as linear or spline interpolation when appropriate.

Finally, it is important to note that when the Fourier transform is required at a single frequency, it is often more efficient to compute it directly using a numerical quadrature method.

Such methods are well suited to handle the oscillatory nature of the Fourier integral and can be more accurate than applying the FFT and then interpolating.

References

Balac, Stéphane. 2011. “La Transformée de Fourier Vue Sous l’angle Du Calcul Numérique.”

Cooley, James W, and John W Tukey. 1965. “An Algorithm for the Machine Calculation of Complex Fourier Series.” Mathematics of Computation 19 (90): 297–301.

El Kolei, Salima. 2013. “Parametric Estimation of Hidden Stochastic Model by Contrast Minimization and Deconvolution: Application to the Stochastic Volatility Model.” Metrika 76 (8): 1031–81.

Epstein, Charles L. 2005. “How Well Does the Finite Fourier Transform Approximate the Fourier Transform?” Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 58 (10): 1421–35.

M. Crouzeix and A.L. Mignot. Analyse numérique des équations différentielles. Collection mathématiques appliquées
pour la maîtrise. Masson, 1992.

A.M. Quarteroni, J.F. Gerbeau, R. Sacco, and F. Saleri. Méthodes numériques pour le calcul scientifique : Programmes
en MATLAB. Collection IRIS. Springer Paris, 2000.

A. Iserles and S. Norsett. On quadrature methods for highly oscillatory integrals and their implementation. BIT
Numerical Mathematics, 44 :755–772, 2004.

A. Iserles and S. Norsett. Efficient quadrature of highly oscillatory integrals using derivatives. Proc. R. Soc. A,
461 :1383–1399, 2005.

Image Credits

All images and visualizations in this article were created by the author using Python (pandas, matplotlib, seaborn, and plotly) and Excel, unless otherwise stated.

Disclaimer

I write to learn, so mistakes are the norm, even though I try my best. Please let me know if you notice any. I am also open to any suggestions for new topics!

---