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

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

import numpy as np

def tfquad(f, nu, n, T):
    """
    Computes the Fourier rework 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 remodel. Should settle for a NumPy array as enter.
    nu : float
        The frequency at which to guage the Fourier rework.
    n : int
        Variety of quadrature factors.
    T : float
        Width of the time window [-T/2, T/2].

    Returns:
    -------
    tfnu : advanced
        Approximated worth of the Fourier rework 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 rework of a perform f at a given frequency nu. The perform tf_integral under implements this strategy. It makes use of numerical integration to compute the Fourier rework of f over the interval [-T/2, T/2].

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

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

along with the actual and imaginary components of its Fourier rework

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

Though f(t) is clean, we are able to 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 rework approximation f̂ seems to be periodic.

  Actually, 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 rework for all frequencies ν ∈ ℝ utilizing the quadrature technique when the parameters T and n are mounted.

Actually, it turns into not possible 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 Ŝₙ.

In consequence, in observe, to compute the Fourier rework for bigger frequencies, one should enhance 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 are able to present that the approximation is dependable at frequency ν when the next situation holds:

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

In accordance with Epstein (2005), when utilizing the Quick Fourier Rework (FFT) algorithm, it’s attainable to precisely compute the Fourier rework 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 rework f̂(ν) of the perform f at some extent ν ∈ [−ν_max/2, ν_max/2], the place ν_max = n/T, that’s,

I now current the Fourier rework algorithm used to approximate f̂(ν).
I can’t go into the small print of the Quick Fourier Rework (FFT) algorithm on this article.
For a simplified clarification, I discuss with Balac (2011), and for a extra technical and complete remedy, to the authentic work by Cooley and Tukey (1965).

You will need to perceive that using the FFT algorithm to approximate the Fourier rework of a perform f relies on the end result established by Epstein (2005). This end result states that when Ŝₙ(ν) is evaluated on the frequencies v_j = j/T, for j = 0, 1, …, n − 1, it supplies an excellent approximation of the continual Fourier rework f̂(ν).

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

Actually, the values of the Fourier rework of f over the interval [−ν_max/2, ν_max/2] may 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 rework Ŝₙ(j/T) may 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 rework of the perform f may 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 rework (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 rework f̂(νⱼ), with ν_j = j/T for j = −n⁄2, …, n⁄2 − 1.

The Python perform tffft under implements these steps to compute the Fourier rework 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 rework 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 rework 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 of Operate")
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 rework 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 throughout the interval [−n/(2T), n/(2T)].
Nevertheless, this technique doesn’t permit us to guage the Fourier rework at arbitrary frequencies ν that aren’t of the shape v_j = j/T.

To compute the Fourier rework of a perform f at a frequency ν that doesn’t match one of many sampled frequencies νⱼ, interpolation strategies can be utilized, akin to 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 rework of f at any frequency ν.

Utilizing Shannon’s Interpolation Theorem to Compute the Fourier Rework of a Operate 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 rework ĝ is zero exterior the interval [−B⁄2, B⁄2] — the perform g may be reconstructed from its samples gₖ = g(ok⁄B) for ok ∈ ℤ.

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

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 rework f̂(ν) for any ν ∈ ℝ.
That is finished through the use of the values of the discrete Fourier rework Ŝₙ(νⱼ) for j = −n/2, …, (n/2 − 1).

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

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

To compute the Fourier rework at any frequency ν, we are able to use Shannon’s interpolation theorem.

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

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

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

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

    Returns:
    - tfnu : approximation of the Fourier rework 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 subsequent perform, fourier_at_nu, combines two steps.
It first computes the discrete Fourier rework of a perform utilizing the FFT (tffft), then applies Shannon interpolation to estimate the Fourier rework at a selected frequency ν.

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

We will now check our implementation utilizing a easy instance.
Right here, we outline a perform f(t) = e^-a|t|, whose Fourier rework is thought analytically.

This lets us evaluate the precise Fourier rework with the approximation obtained utilizing our technique.

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

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

We consider the Fourier rework 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 though the frequency 3⁄T (whose Fourier rework worth seems within the FFT output) is near π⁄T, the corresponding values of the Fourier rework of f at these two frequencies are fairly completely different.

This reveals that Shannon’s interpolation method may be very helpful when the specified Fourier rework 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 rework 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, not like the built-in fft features in SciPy or NumPy, the FFT may be tailored to compute the Fourier rework of a steady perform by means of correct formulation.

Our implementation was based mostly 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 rework of any perform on ℝ at arbitrary frequencies.

This interpolation may be changed by extra conventional approaches akin to linear or spline interpolation when applicable.

Lastly, it is very important word that when the Fourier rework is required at a single frequency, it’s typically extra environment friendly to compute it straight utilizing a numerical quadrature technique.

Such strategies are nicely suited to deal with the oscillatory nature of the Fourier integral and may 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 Advanced Fourier Sequence.” 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 photographs and visualizations on this article have been created by the creator utilizing Python (pandas, matplotlib, seaborn, and plotly) and Excel, until in any other case said.

Disclaimer

I write to study, so errors are the norm, though I strive my finest. Please let me know if you happen to discover any. I’m additionally open to any options for brand new subjects!