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
= 1024
N = 1.0 / 64
T = np.linspace(-N/2*T, N/2*T, N, endpoint=False)
t = f(t)
y
# FFT with scipy
= fftshift(fft(y)) * T
yf_scipy = fftshift(fftfreq(N, T))
xf = np.exp(-np.pi * xf**2)
FT_exact
# FFT with numpy
= np.fft.fftshift(np.fft.fft(y)) * T
yf_numpy = np.fft.fftshift(np.fft.fftfreq(N, T))
xf_numpy
# Plot with subplot_mosaic
= plt.subplot_mosaic([["scipy", "numpy"]], figsize=(7, 5), layout="constrained", sharey=True)
fig, axs
# Scipy FFT
"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)
axs[
# NumPy FFT
"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)
axs[
"Comparison of FFT Implementations vs. Exact Fourier Transform", fontsize=14)
plt.suptitle( plt.show()
Fourrier Transformation
Introduction
The primary goal of the Fourier transformation is to analyze signals by decomposing them into their constituent frequencies. It was Joseph Fourier who first recognized the significance of spectral decomposition of a signal. Indeed, he demonstrated that any periodic signal can be decomposed into a finite sum of sinusoidal signals with constant frequencies and amplitudes. The finite set of these frequencies constitutes the spectrum of the signal.
Beyond spectral analysis, the Fourier transformation is employed in solving partial differential equations, evaluating integrals, and computing sums of series. It has applications in various fields, including physics, engineering, and signal processing. In recent years, it has found new uses in finance, such as estimating financial asset volatility, which We will cover in future posts.
In this article, we explore the numerical computation of the Fourier transform for both real and complex functions. Drawing on the work of Balac (2011), who proposed a quadrature-based approach, we illustrate how this method can be applied. We also demonstrate how the Fast Fourier Transform (FFT) algorithm (for more details see Cooley and Tukey 1965) enables efficient computation of the discrete Fourier transform, making it well-suited for calculating the Fourier transform of integrable functions and the Fourier coefficients of periodic functions.
This article was motivated by the realization that the numerical methods implemented in Python—such as [scipy.fft
and numpy.fft
] do not directly compute the Fourier transform of a function. This observation had already been made by Balac (2011) in the context of MATLAB.
Definition and Properties of the Fourier Transform
We must first define the space in which the function we want to compute the Fourier transform is defined and the caracteristics they must satisfy so that the Fourier transform exists. These functions are defined in the space of integrable functions, denoted by \(L^1(\mathbb{R,K})\) which are functions \(f:\mathbb{R} \to \mathbb{K}\) [where \(\mathbb{K}\) is either \(\mathbb{R}\) or \(\mathbb{C}\)]. These functions are integrable in the sense of Lebesgue, meaning that the integral of their absolute value is finite: \[ \int_{-\infty}^{+\infty} |f(t)| dt < +\infty. \]
So for a function \(f\) to be in \(L^1(\mathbb{R,K})\), the product \(f(t) \cdot e^{-2i\pi\nu t}\) is integrable for all \(\nu \in \mathbb{R}\). And for all \(\nu \in \mathbb{R}\), the function \(\hat{f}\) [notée egalement \(F(f)\)] defined by \[ \hat{f}(\nu) = \int_{-\infty}^{+\infty} f(x) e^{-2i\pi\nu t} dt \]
is well-defined and is called the Fourier transform of \(f\).
We deduce in this formula that a Fourier transform is a complex-valued function, and it is a linear operator. We can also deduce others properties of the Fourier transform in particular for a function \(f\) in \(L^1(\mathbb{R,R})\), une fonction à valeurs réelles :
Translation property: If \(f\) is in \(L^1(\mathbb{R,K})\), \(t_0 \in \mathbb{R}\), and \(h : t \in \mathbb{R} \mapsto f(t - t_0)\), we have \(h \in L^1(\mathbb{R,K})\) and \(\hat{h}(\nu) = e^{-2i\pi\nu t_0} \hat{f}(\nu)\).
Homothety property: If \(f\) is in \(L^1(\mathbb{R,K})\) and \(a \in \mathbb{R}\), then the function \(g : t \in \mathbb{R} \mapsto f(at)\) is in \(L^1(\mathbb{R,K})\) and \(\hat{g}(\nu) = \frac{1}{|a|} \hat{f}\left(\frac{\nu}{a}\right)\).
Modulation property: If \(f\) is in \(L^1(\mathbb{R,K})\) and \(\nu_0 \in \mathbb{R}\), and \(h : t \in \mathbb{R} \mapsto f(t) e^{2i\pi\nu_0 t}\), we have \(h \in L^1(\mathbb{R,K})\) and \(\hat{h}(\nu) = \hat{f}(\nu - \nu_0)\).
If \(f \in L^1(\mathbb{R,R})\) and est pair, then \(\hat{f}\) is also real-valued and the real part of \(\hat{f}\) is also an even function.
If \(f\) is in \(L^1(\mathbb{R,R})\) and is odd, then \(\hat{f}\) is est une fonction imaginaire pure and the imaginary part of \(\hat{f}\) is also an odd function.
For some functions, the Fourier transform can be computed analytically. For example, for the function \(f : t \in \mathbb{R} \mapsto \mathbb{1}_{[-\frac{a}{2}, \frac{a}{2}]}(t)\), the Fourier transform is given by: \[ \hat{f} : \nu \in \mathbb{R} \mapsto a sinc(a \pi \nu) \] where \(sinc(t) = \frac{\sin(t)}{t}\) is the sinc function.
However, for many functions, the Fourier transform cannot be computed analytically. In such cases, we can use numerical methods to approximate the Fourier transform. We will explore these numerical dans la suite of this article.
How to approximate the Fourier Transform?
The Fourier transform of a function \(f\) is defined as an integral over the entire real line. However, for the functions that are integral in the sense of lebesgue and that have a practical applications tend to 0 as \(|t| \to +\infty\). And we can approximate the Fourier transform by integrating over a finite interval \([-T, T]\). If the lenght of the interval is large enough, or if the function decays quickly when t tends to infinity, this approximation will be accurate.
\[ \hat{f}(\nu) \approx \int_{-T}^{T} f(t) e^{-2i\pi\nu t} dt \]
In his article, Balac (2011) va plus loin en montrant qu’approximating the Fourier transform met en oeuvre trois outils mathématiques :
- Les séries de Fourier dans le cadre d’un signal périodique.
- La transformée de Fourier dans le cadre d’un signal non périodique.
- La transformée de Fourier discrète dans le cadre d’un signal discret.
And for each of these tools, computing the Fourier transform revient à calculer l’intgrale below: \[ \hat{f}(\nu) \approx \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) e^{-2i\pi\nu t} dt \] I recommends reading his article for more details on these three tools. Then we will focus on the numerical computation of the Fourier transform using quadrature methods, which is a numerical integration technique.
Numerical Computation of the Fourier Transform
We show that to compute the Fourier transform of a function \(f\) consists to approximating it by the integral below in the interval \([-\frac{T}{2}, \frac{T}{2}]\): \[ \underbrace{ \int_{-\infty}^{+\infty} f(t)\, e^{-2\pi i \nu t} \, dt }_{=\,\hat{f}(\nu)} \approx \underbrace{ \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\, e^{-2\pi i \nu t} \, dt }_{=\,\tilde{S}(\nu)} \]
where T is a large enough number such that the integral converges. Une valeur approchée of the integral \(\tilde{S}(\nu) = \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\, e^{-2\pi i \nu t} \, dt\) can be computed using quadrature methods. In the next section, we will approximate the integral using the quadrature method of rectangles à gauche.
Quadrature Method of Rectangles à Gauche
For computing the integral \(\tilde{S}(\nu) = \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\, e^{-2\pi i \nu t} \, dt\), using the quadrature method of rectangles à gauche, we must respect the following steps: 1. Discretization of the Interval: We discretize the interval \([-\frac{T}{2}, \frac{T}{2}]\) into \(N\) uniform subintervals of length \(h_t = \frac{T}{N}\). The points in the interval are given by:
\[ t_k = -\frac{T}{2} + k \cdot h_t, \quad k = 0, 1, \ldots, N-1. \]
- Approximation of the Integral: Using the Chasles relation, we can approximate the integral \(\tilde{S}(\nu)\) as follows:
\[ \tilde{S}(\nu) = \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\, e^{-2\pi i \nu t} \, dt = \sum_{k=0}^{N-1} \int_{t_k}^{t_{k+1}} f(t)\, e^{-2\pi i \nu t} \, dt. \]
By taking into account that we have \(t_{k+1} - t_k = h_t\), and \(t_k = -\frac{T}{2} + k \cdot h_t = T(\frac{k}{N} - \frac{1}{2})\), we can rewrite the integral as: \[ \tilde{S}(\nu) = \sum_{k=0}^{N-1} f(t_k) e^{-2\pi i \nu t_k} h_t. \] 3. Final Formula: The final formula for the approximation of the Fourier transform is given by:
\[ \boxed{ \forall \nu \in \mathbb{R} \quad \underbrace{ \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t)\, e^{-2\pi i \nu t} \, dt }_{=\,\tilde{S}(\nu)} \approx \underbrace{ \frac{T}{n} e^{i \pi \nu T} \sum_{k=0}^{n-1} f_k\, e^{-2 i \pi \nu T k / n} }_{=\,\tilde{S}_n(\nu)} \quad \text{où } f_k = f\left( \frac{2k - n}{2n} T \right). } \]
Mis en Oeuvre in Python
The function tfquad
below implements the quadrature method of rectangles à gauche 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.
"""
= np.arange(n)
k = (k / n - 0.5) * T
t_k = np.exp(-2j * np.pi * nu * T * k / n)
weights = (T / n) * np.exp(1j * np.pi * nu * T)
prefactor
return prefactor * np.sum(f(t_k) * weights)
We can also use the function scipy’s quad
functionnfor defining 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 a function 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."""
= quad(lambda t: np.real(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
real_part = quad(lambda t: np.imag(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
imag_part return real_part + 1j * imag_part
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 -----
= 2.0
T = 32
n = np.linspace(-6, 6, 500)
nu_vals = exact_fourier_transform(nu_vals)
exact_vals = np.array([tfquad(f, nu, n, T) for nu in nu_vals])
tfquad_vals
# Compute the approximation using scipy integral
= np.array([tf_integral(f, nu, T) for nu in nu_vals]) tf_integral_vals
= plt.subplot_mosaic([["tfquad", "quad"]], figsize=(7, 5), layout="constrained")
fig, axs
# Plot using tfquad implementation
"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)
axs[
# Plot using scipy.integrate.quad
"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)
axs[
# --- Global legend below the plots ---
# Take handles from one subplot only (assumes labels are consistent)
= axs["quad"].get_legend_handles_labels()
handles, labels
fig.legend(handles, labels,='lower center', bbox_to_anchor=(0.5, -0.05),
loc=3, frameon=False)
ncol
"Comparison of FFT Implementations vs. Exact Fourier Transform", fontsize=14)
plt.suptitle(=0.2) # Make room for the legend
plt.subplots_adjust(bottom plt.show()
/var/folders/v8/l5q0bw4s2ln17s59y7cc86rm0000gn/T/ipykernel_70773/938168490.py:30: UserWarning:
This figure was using a layout engine that is incompatible with subplots_adjust and/or tight_layout; not calling subplots_adjust.
We are able now to compute the Fourier transform of a function using the quadrature method of rectangles. Let’s characterize this approximation.
Characterization of the Approximation
Le caractère oscillant de l’approximation
We observe that the transform \(\hat{f}\) of the function \(f\) is oscillatory. This oscillatory nature is due to the complex exponential term \(e^{-2\pi i \nu t}\) in the integral.
To illustrate this, the figure below represents the function \(f : t \in \mathbb{R} \mapsto e^{-t^2} \in \mathbb{R}\), ainsi que les parties réelle et imaginaire de sa transformée de Fourier \(\hat{f} : \nu \in \mathbb{R} \mapsto \hat{f}(\nu) \in \mathbb{C}\), for \(\nu = \frac{5}{2}\). Bien que la fonction \(f\) soit lisse, on observe des fortes oscillations pour la fonction \(\hat{f}\),.
import numpy as np
import matplotlib.pyplot as plt
= 5 / 2
nu = np.linspace(-8, 8, 1000)
t1 = np.linspace(-4, 4, 1000)
t2
= lambda t: np.exp(-t**2)
f = lambda t: f(t) * np.exp(-2j * np.pi * nu * t)
phi
= f(t1)
f_vals = phi(t2)
phi_vals
# Plot
= plt.subplots(1, 2, figsize=(10, 4))
fig, 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)
axs[
plt.tight_layout() plt.show()
Ces fortes variations peuvent poser des problèmes pour l’approximation numérique de la transformée de Fourier en utilisant des méthodes de quadrature même lorsqu’on prend un nombre de points \(n\) assez grand. L’utilisation de l’agorithme de la transformée de Fourier rapide (FFT) est une solution pour surmonter ce problème.
The approximation using the quadrature method of rectangles à gauche is periodic.
On remarque que même si la fonction \(f\) n’est pas périodique, la transformée de Fourier \(\hat{f}\) est périodique. En effet, la fonction \(\hat{S}_n\) est périodique de période \(T = \frac{n}{T}\) : \[ \forall \nu \in \mathbb{R} \quad \widehat{S}_n\left(\nu + \frac{n}{T} \right) = \frac{T}{n} e^{i\pi \nu T} \sum_{k=0}^{n-1} f_k \, e^{-2i\pi\left(\nu + \frac{n}{T} \right)T \frac{k}{n}} = \frac{T}{n} e^{i\pi \nu T} \sum_{k=0}^{n-1} f_k \, e^{-2i\pi \nu T \frac{k}{n}} \underbrace{e^{-2i\pi k}}_{=1} \]
\[ \phantom{\forall \nu \in \mathbb{R} \quad \widehat{S}_n\left(\nu + \frac{n}{T} \right)} = \frac{T}{n} e^{i\pi \nu T} \sum_{k=0}^{n-1} f_k \, e^{-2i\pi \nu T \frac{k}{n}} = \widehat{S}_n(\nu). \]
This periodicity of \(\hat{S}_n\) has pour consequence that it is not possible to compute the Fourier transform by the quadrature method for all frequencies \(\nu \in \mathbb{R}\) when the parameter \(T\) and \(n\) are fixed. In fact, it is impossible to compute \(\hat{f}(\nu)\) when \(\nu \geq {\nu}_{max}\) or \(\nu \leq -{\nu}_{max}\) where \({\nu}_{max} = \frac{n}{T}\) is the period of the Fourier transform \(\hat{S}_n\). So in practice, for computing the Fourier transform of a function \(f\) for \(\nu\) very large, we need to increase the parameter \(T\) or \(n\). In fact, by evaluating the error of the approximation of the Fourier transform \(\hat{f}\) by the quadrature method of rectangles à gauche, we can see that we have a good approximation au point \(\nu\) when this relation holds:
\[ \frac{\nu T}{n} \ll 1. \]
Epstein (2005) shows that when using the algorithm of the Fast Fourier Transform (FFT), we can correctly compute the Fourier transform of a function \(f\) for all frequencies \(\nu \in \mathbb{R}\), even when \(\frac{\nu T}{n}\) is proche to 1 for toute fonction \(f\) continue par morceaux et à support borné.
Utilisation the l’algorithme de la Transformée de Fourier Rapide (FFT) pour calculer la transformée de Fourier d’une fonction f en un point \(\nu\).
Dans cette partie, nous considérons toujours que \(\hat{S}_n(\nu)\) est l’approximation de la transformée de Fourier \(\hat{f}(\nu)\) de la fonction \(f\) en un point \(\nu \in [-\frac{\nu_{max}}{2}, \frac{\nu_{max}}{2}]\) où \(\nu_{max} = \frac{n}{T}\) ie \[ \hat{f}(\nu) \approx \hat{S}_n(\nu) = \frac{T}{n} e^{i\pi \nu T} \sum_{k=0}^{n-1} f_k\, e^{-2 i \pi \nu T k / n}. \]
Nous allons présenter l’algorithme de la transformée qui permet d’approximer \(\hat{f}(\nu)\). Nous n’allons pas détailler l’algorithme de la Transformée de Fourier Rapide (FFT) dans cet article. Balac (2011) provides a simplified explanation of the FFT algorithm. Pour plus de détails, je vous recommande de lire l’article de Cooley and Tukey (1965).
Ce qu’il faut retenir est que l’utilisation de l’algorithme de la Transformée de Fourier Rapide (FFT) pour calculer la transformée de Fourier d’une fonction \(f\) en un point \(\nu\) se fait en deux étapes :
La première étape utilise le fait que Epstein (2005) shows that when \(\hat{S_n}(\nu)\) is evaluated at \(\nu = \frac{j}{T}\) for \(j = 0, 1, \ldots, n-1\), we have the good approximation of the Fourier transform \(\hat{f}(\nu)\).
De plus nous savons que \(\hat{S_n}\) est périodique. Et Cette périodicité fait jouer un rôle symétrique aux indices \(j \in \{0, 1, \ldots, n-1\}\) et \(k \in \{-n/2, -n/2 + 1, \ldots, -1\}\). En effet, les valeurs de la transformée de Fourier de f sur l’intervalle \([-\frac{\nu_{max}}{2}, \frac{\nu_{max}}{2}]\) peuvent être déduites des valeurs de \(\hat{S_n}\) aux points \(\nu_j = \frac{j}{T}\) pour \(j = 0, 1, \ldots, n-1\) comme suit :
\[ \widehat{S}_n(\nu'_j) = \frac{T}{n} (-1)^j \sum_{k=0}^{n-1} f_k\, e^{-2i\pi j \frac{k}{n}} = \begin{cases} \widehat{S}_n(\nu_j) & \text{si } j \in \left\{0, \dots, \frac{n}{2} - 1 \right\} \\ \widehat{S}_n(\nu_{j-n}) & \text{si } j \in \left\{ \frac{n}{2}, \dots, n-1 \right\} \end{cases} \]
\[ \text{où on a exploité la relation} \quad e^{-2i\pi j \frac{k}{n}} = e^{-2i\pi (j-n) \frac{k}{n}} \times \underbrace{e^{-2i\pi k}}_{=1} = e^{-2i\pi (j-n) \frac{k}{n}} \quad \text{pour } j \in \left\{ \frac{n}{2}, \dots, n-1 \right\}. \]
Cette relation nous montre que l’on peut calculer la transformée de Fourier \(\hat{S_n}(\frac{j}{T})\) pour \(j = -\frac{n}{2}, \ldots, \frac{n}{2} - 1\). De plus, lorsque \(n\) est une puissance de 2, calcul est plus rapide. Ce procédé est appelé la Transformée de Fourier Rapide (FFT).
Si nous récapitulons, nous avons montré que nous pouvons approximer la transformée de Fourier de la fonction \(f\) dans l’intervalle \([-\frac{T}{2}, \frac{T}{2}]\) en les fréquences \(\nu_j = \frac{j}{T}\) pour \(j = -\frac{n}{2}, \ldots, \frac{n}{2} - 1\) avec \(n = 2^m\) pour un entier \(m \geq 0\) en utilisant l’algorithme de la Transformée de Fourier Rapide (FFT) en procédant comme suit :
Créer la suite finie F de valeurs \(f(\frac{2k - n}{2n} T)\) pour \(k = 0, 1, \ldots, n-1\).
Calculer la transformée de Fourier discrète \(\hat{F}\) de la suite F en utilisant l’algorithme de la Transformée de Fourier Rapide (FFT) qui est donnée par \(\sum_{k=0}^{n-1} f_k e^{-2i\pi \frac{jk}{n}}\) pour \(j = 0, 1, \ldots, n-1\).
Symmétriser les valeurs de la seconde partie de \(-\frac{n}{2}, \ldots, -1\);
Multiplier chacune des valeurs du tableau par \(\frac{T}{n} (-1)^{j-1}\) où \(j \in \{1, \ldots, n\}\).
On dispose ainsi d’un tableau correspondant aux valeurs de la transformée de Fourier \(\hat{f}(\nu_j)\) pour \(\nu_j = \frac{j}{T}\) pour \(j = -\frac{n}{2}, \ldots, \frac{n}{2} - 1\).
Donc dans la première étape, consiste à discretiser l’intervalle \([-\frac{T}{2}, \frac{T}{2}]\) en \(n\) points \(\nu_j = \frac{j}{T}\) pour \(j = 0, 1, \ldots, n-1\). En effet on a :
La fonction matlab tffft suivante calcule la transformée de Fourier en python d’une fonction donnée en mettant en œuvre les différentes étapes qui ont été détaillées précédemment.
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.
nu : np.ndarray
Fréquences discrètes correspondantes (de -n/(2T) à (n/2 - 1)/T).
"""
= T / n
h = -0.5 * T + np.arange(n) * h # noeuds temporels
t = f(t) # échantillonnage de f
F = h * (-1) ** np.arange(n) * fftshift(fft(F)) # TF approximée
tf = -n / (2 * T) + np.arange(n) / T # fréquences ν_j = j/T
nu
return tf, nu, t
On illustre dans le programme ci-dessous le calcul de la transformée de Fourier d’une fonction donnée (une fonction de Gauss).
# Paramètres
= 10
a = lambda t: np.exp(-a * t**2)
f = 10
T = 2**8 # 256
n
# Calcul de la transformée de Fourier
= tffft(f, T, n)
tf, nu, t
# Représentation graphique
= plt.subplots(1, 2, figsize=(10, 4))
fig, axs
0].plot(t, f(t), '-g', linewidth=3)
axs[0].set_xlabel("temps")
axs[0].set_title("Fonction considérée")
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("fréquence")
axs[1].set_title("Transformée de Fourier")
axs[1].set_xlim(-15, 15)
axs[1].set_ylim(-0.5, 1)
axs[1].grid(True)
axs[
plt.tight_layout() plt.show()
Dans ce cas, lorsque n est une puissance de 2 \(\hat{S_n}\) est approximée rapidement. Et c’est l’algorithme qui permet d’estimer \(\hat{S_n}\) au point \(\nu_j = \frac{j}{T}\) pour \(j = 0, 1, \ldots, n-1\) qui est appelé la Transformée de Fourier Rapide (FFT).