Implementing the Fourier Rework Numerically in Python: A Step-by-Step Information

The perform tfquad beneath implements the left rectangle quadrature technique to compute the Fourier remodel of a perform f at a given frequency nu.

import numpy as np

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

    Parameters:
    ----------
    f : callable
        The perform to rework. Should settle for a NumPy array as enter.
    nu : float
        The frequency at which to judge the Fourier remodel.
    n : int
        Variety of quadrature factors.
    T : float
        Width of the time window [-T/2, T/2].

    Returns:
    -------
    tfnu : complicated
        Approximated worth of the Fourier remodel at frequency nu.
    """
    ok = np.arange(n)
    t_k = (ok / n - 0.5) * T
    weights = np.exp(-2j * np.pi * nu * T * ok / n)
    prefactor = (T / n) * np.exp(1j * np.pi * nu * T)

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

We will additionally use SciPy’s quad perform to outline the Fourier remodel of a perform f at a given frequency nu. The perform tf_integral beneath implements this strategy. It makes use of numerical integration to compute the Fourier remodel of f over the interval [-T/2, T/2].

from scipy.combine import quad

def tf_integral(f, nu, T):
    """Compute FT of f at frequency nu over [-T/2, T/2] utilizing scipy quad."""
    real_part = quad(lambda t: np.actual(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 components for the left rectangle quadrature technique, we will now implement it in Python to see the way it performs in follow.

To do that, we take into account a easy instance the place the perform fff is outlined as an indicator perform on the interval [−1, 1]:

For this perform, the Fourier remodel could be computed analytically, which permits us to examine the numerical approximation with the precise outcome.

The next Python script implements each the analytical Fourier remodel and its numerical approximations utilizing:

• the left rectangle quadrature technique, and
• SciPy’s integration perform for reference.

import numpy as np
import matplotlib.pyplot as plt

# ----- Perform Definitions -----

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

def exact_fourier_transform(nu):
    """Analytical FT of the indicator perform 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 utilizing 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, format="constrained")

# Plot utilizing tfquad implementation
axs["tfquad"].plot(nu_vals, np.actual(exact_vals), 'b-', linewidth=2, label=r'$hat{f}$ (precise)')
axs["tfquad"].plot(nu_vals, np.actual(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 utilizing scipy.combine.quad
axs["quad"].plot(nu_vals, np.actual(exact_vals), 'b-', linewidth=2, label=r'$hat{f}$ (quad)')
axs["quad"].plot(nu_vals, np.actual(tf_integral_vals), 'r--', linewidth=1.5, label=r'approximation $hat{S}_n$')
axs["quad"].set_title("TF avec scipy.combine.quad")
axs["quad"].set_xlabel(r'$nu$')
axs["quad"].set_ylabel('Amplitude')
axs["quad"].grid(False)
axs["quad"].set_ylim(-0.5, 2.1)



# --- World legend beneath the plots ---
# Take handles from one subplot solely (assumes labels are constant)
handles, labels = axs["quad"].get_legend_handles_labels()
fig.legend(handles, labels,
           loc='decrease heart', bbox_to_anchor=(0.5, -0.05),
           ncol=3, frameon=False)

plt.suptitle("Comparability of FFT Implementations vs. Actual Fourier Rework", fontsize=14)

plt.present()

2.1.2 Characterization of approximation utilizing the left rectangle quadrature technique

As an example this habits, the determine beneath reveals the perform

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

along with the actual and imaginary components of its Fourier remodel

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

Though f(t) is clean, we will observe robust oscillations in f̂(ν), which reveal the affect of the exponential time period e^-2πiνt within 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, 'ok', 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.actual(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^{-2ipinu t}$, $nu=5/2$")
axs[1].legend()
axs[1].grid(True)

plt.tight_layout()
plt.present()

• The approximation obtained utilizing the left rectangle quadrature technique is periodic in nature.

  We observe that even when the perform f is not periodic, its Fourier remodel approximation f̂ seems to be periodic.

  In reality, the perform Ŝₙ, obtained from the quadrature technique, is periodic with a interval given by:

This periodicity of Ŝₙ implies that it’s not attainable to compute the Fourier remodel for all frequencies ν ∈ ℝ utilizing the quadrature technique when the parameters T and n are mounted.

In reality, it turns into unattainable to compute f̂(ν) precisely when

|ν| ≥ ν_max,

the place v = n/T represents the most frequency that may be resolved due to the periodic nature of Ŝₙ.

Because of this, in follow, to compute the Fourier remodel for bigger frequencies, one should improve both the time window (T) or the variety of factors (n).

Moreover, by evaluating the error in approximating f̂(ν) utilizing the left rectangle quadrature technique, we will present that the approximation is dependable at frequency ν when the next situation holds:

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

In line with Epstein (2005), when utilizing the Quick Fourier Rework (FFT) algorithm, it’s attainable to precisely compute the Fourier remodel of a perform f for all frequencies ν ∈ ℝ, even when (|ν|T)⁄n is near 1, offered that f is piecewise steady and has compact assist.

2.2 Computing the Fourier Rework at Frequency ν Utilizing the FFT Algorithm

On this part, we denote by Ŝₙ(ν) the approximation of the Fourier remodel f̂(ν) of the perform f at some extent ν ∈ [−ν_max/2, ν_max/2], the place ν_max = n/T, that’s,

I now current the Fourier remodel algorithm used to approximate f̂(ν).
I cannot go into the small print of the Quick Fourier Rework (FFT) algorithm on this article.
For a simplified rationalization, I check with Balac (2011), and for a extra technical and complete remedy, to the unique work by Cooley and Tukey (1965).

You will need to perceive that the usage of the FFT algorithm to approximate the Fourier remodel of a perform f is predicated on the outcome established by Epstein (2005). This outcome states that when Ŝₙ(ν) is evaluated on the frequencies v_j = j/T, for j = 0, 1, …, n − 1, it supplies a very good approximation of the continual Fourier remodel f̂(ν).

Furthermore, Ŝₙ is understood to be periodic. This periodicity assigns a symmetric function to the indices j ∈ {0, 1, …, n − 1} and ok ∈ {−n/2, −n/2 + 1, …, −1}.

In reality, the values of the Fourier remodel of f over the interval [−ν_max/2, ν_max/2] could be derived from the values of Ŝₙ on the factors ν_j = j/T, for j = 0, 1, …, n − 1, as follows:

the place we used the relation:

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

Furthermore, when n is an influence of two, the computation turns into considerably quicker (see Balac, 2011). This course of is named the Quick Fourier Rework (FFT).

To summarize, I’ve proven that the Fourier remodel of the perform f could be approximated over the interval [−T⁄2, T⁄2] on the frequencies v_j = j/T, for j = −n/2, …, n/2 − 1, the place n = 2ᵐ for some integer m ≥ 0, by making use of the Quick Fourier Rework (FFT) algorithm as follows:

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

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

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

The Python perform tffft beneath implements these steps to compute the Fourier remodel of a given perform.

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 à assist 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 factors 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 next program illustrates the right way to compute the Fourier remodel of the Gaussian perform f(t) = e^{-10t^2} over the interval [-10, 10], utilizing the Quick Fourier Rework (FFT) algorithm.

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

# Compute the Fourier remodel utilizing 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("Thought-about Perform")
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 Rework utilizing FFT")
axs[1].set_xlim(-15, 15)
axs[1].set_ylim(-0.5, 1)
axs[1].grid(True)

plt.tight_layout()
plt.present()

The tactic I simply introduced makes it attainable to compute and visualize the Fourier remodel of a perform f at discrete factors v_j = j/T, for j = −n/2, …, n/2 − 1, the place n is an influence of two.

These factors lie inside the interval [−n/(2T), n/(2T)].
Nonetheless, this technique doesn’t permit us to judge the Fourier remodel at arbitrary frequencies ν that aren’t of the shape v_j = j/T.

To compute the Fourier remodel of a perform f at a frequency ν that doesn’t match one of many sampled frequencies νⱼ, interpolation strategies can be utilized, resembling linear, polynomial, or spline interpolation.

On this article, we use the strategy proposed by Balac (2011), which depends on Shannon’s interpolation theorem to compute the Fourier remodel of f at any frequency ν.

Utilizing Shannon’s Interpolation Theorem to Compute the Fourier Rework of a Perform f at a Level ν

What does Shannon’s theorem inform us?

It states that for a band-limited perform g — that’s, a perform whose Fourier remodel ĝ is zero outdoors the interval [−B⁄2, B⁄2] — the perform g could be reconstructed from its samples gₖ = g(ok⁄B) for ok ∈ ℤ.

If we let ν_c be the smallest constructive quantity such that ĝ(ν) is zero outdoors the interval [−2πν_c, 2πν_c], then the Shannon interpolation components applies.
For each t ∈ ℝ, and any constructive actual quantity α ≥ 1/(2ν_c), we now have:

the place sinc(x) = sin(x)/x is the sinc perform (cardinal sine).

Balac (2011) reveals that when a perform f has a bounded assist within the interval [−T⁄2, T⁄2], Shannon’s interpolation theorem can be utilized to compute its Fourier remodel f̂(ν) for any ν ∈ ℝ.
That is performed through the use of the values of the discrete Fourier remodel Ŝₙ(νⱼ) for j = −n/2, …, (n/2 − 1).

By setting α = 1/T, Balac derives the next Shannon interpolation components, legitimate for all ν ∈ ℝ:

The next program illustrates the right way to use Shannon’s interpolation theorem to compute the Fourier remodel of a perform f at a given frequency ν.

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

The concept is to estimate the worth of the Fourier remodel from its discrete samples obtained by the FFT.

The perform beneath, shannon, implements this interpolation technique in Python.

def shannon(tf, nu, T):
    """
    Approximates the worth of the Fourier remodel of perform f at frequency 'nu'
    utilizing its discrete values computed from the FFT.

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

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

    return tfnu

The following perform, fourier_at_nu, combines two steps.
It first computes the discrete Fourier remodel of a perform utilizing the FFT (tffft), then applies Shannon interpolation to estimate the Fourier remodel at a particular frequency ν.

This makes it attainable to judge the remodel at any arbitrary level ν, not simply on the FFT sampling frequencies.

We will now take a look at our implementation utilizing a easy instance.
Right here, we outline a perform f(t) = e^-a|t|, whose Fourier remodel is understood analytically.

This lets us examine the precise Fourier remodel with the approximation obtained utilizing our technique.

a = 0.5
f = lambda t: np.exp(-a * np.abs(t))                          # Perform to rework
fhat_exact = lambda nu: (2 * a) / (a**2 + 4 * np.pi**2 * nu**2)  # Actual Fourier remodel

T = 40     # Time window
n = 2**10  # Variety of discretization factors

We consider the Fourier remodel at two frequencies: ν = 3/T and ν = π/T.

For every case, we print each the precise worth and the approximation obtained by our fourier_at_nu perform.

This comparability helps confirm 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"Actual worth at nu={nu}: {exact_value}")
print(f"Approximation at nu={nu}: {np.actual(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"Actual worth at nu={nu}: {exact_value}")
print(f"Approximation at nu={nu}: {np.actual(approx_value)}")

Actual worth at ν = 3/T: 2.118347413776218  
Approximation at ν = 3/T: (2.1185707502943534)

Actual worth at ν = π/T: 2.0262491352594427  
Approximation at ν = π/T: (2.0264201680784835)

It’s value noting that despite the fact that the frequency 3⁄T (whose Fourier remodel worth seems within the FFT output) is near π⁄T, the corresponding values of the Fourier remodel of f at these two frequencies are fairly totally different.

This reveals that Shannon’s interpolation components could be very helpful when the specified Fourier remodel values are in a roundabout way included among the many outcomes produced by the FFT algorithm.

Conclusion

On this article, we explored two strategies to approximate the Fourier remodel of a perform.

The primary was a numerical quadrature technique (the left rectangle rule), and the second was the Quick Fourier Rework (FFT) algorithm.

We confirmed that, in contrast to the built-in fft capabilities in SciPy or NumPy, the FFT could be tailored to compute the Fourier remodel of a steady perform via correct formulation.

Our implementation was primarily based on the work of Balac (2013), who demonstrated the right way to reproduce the FFT computation in Python.
We additionally launched Shannon’s interpolation theorem, which permits us to estimate the Fourier remodel of any perform on ℝ at arbitrary frequencies.

This interpolation could be changed by extra conventional approaches resembling linear or spline interpolation when applicable.

Lastly, you will need to be aware that when the Fourier remodel is required at a single frequency, it’s typically extra environment friendly to compute it straight utilizing a numerical quadrature technique.

Such strategies are properly suited to deal with the oscillatory nature of the Fourier integral and could be extra correct than making use of the FFT after which 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 Complicated Fourier Collection.” Arithmetic of Computation 19 (90): 297–301.

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

Epstein, Charles L. 2005. “How Properly Does the Finite Fourier Rework Approximate the Fourier Rework?” Communications on Pure and Utilized Arithmetic: 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. Assortment 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. Assortment IRIS. Springer Paris, 2000.

A. Iserles and S. Norsett. On quadrature strategies for extremely oscillatory integrals and their implementation. BIT
Numerical Arithmetic, 44 :755–772, 2004.

A. Iserles and S. Norsett. Environment friendly quadrature of extremely oscillatory integrals utilizing derivatives. Proc. R. Soc. A,
461 :1383–1399, 2005.

Picture Credit

All pictures and visualizations on this article had been created by the writer utilizing Python (pandas, matplotlib, seaborn, and plotly) and Excel, until in any other case acknowledged.

Disclaimer

I write to be taught, so errors are the norm, despite the fact that I strive my finest. Please let me know if you happen to discover any. I’m additionally open to any ideas for brand new subjects!