I'm implementing a discrete inverse Fourier transform in Python to approximate the inverse Fourier transform of a Gaussian function.
The input function is sqrt(pi) * e^(-w^2/4)
so the output must be e^(-x^2)
.
While the shape of the resulting function looks correct, the x-axis scaling seems to be off (There might be just some normalization issue). I expect to see a Gaussian function of the form e^(-x^2)
, but my result is much narrower.
This is my implementation:
import matplotlib.pyplot as plt
import numpy as np
from sympy import symbols, exp, pi, lambdify, sqrt
# Defining the Fourier transform of a Gaussian function, sqrt(pi) * exp(-omega ** 2 / 4)
x, omega = symbols('x omega')
f_gaussian_symbolic = exp(-omega ** 2 / 4) * sqrt(pi)
f_gaussian_function = lambdify(omega, f_gaussian_symbolic, 'numpy')
def fourier_inverse(f, n, omega_max):
"""
This function computes the inverse Fourier transform of a function f.
:param f: The function to be transformed
:param n: Number of samples
:param omega_max: The max frequency we want to be sampled
"""
omega_range = np.linspace(-omega_max, omega_max, n)
f_values = f(omega_range)
inverse_f = np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f_values)))
delta_omega = omega_range[1] - omega_range[0]
x_range = np.fft.ifftshift(np.fft.fftfreq(n, d=delta_omega))
inverse_f *= delta_omega * n / (2 * np.pi)
return x_range, inverse_f
plt.figure(figsize=(10, 5))
x_range, inverse_f = fourier_inverse(f_gaussian_function, 10000, 100)
plt.plot(x_range, inverse_f.real)
plt.ylim(-2, 2)
plt.xlim(-4, 4)
plt.show()
The shape of the function looks correct, but it's much narrower than expected. I suspect there's an issue with how I'm calculating or scaling the x_range
in my fourier_inverse
function.
What am I doing wrong in my implementation, and how can I correct the x-axis scaling to get the expected Gaussian function e^(-x^2)?
It looks like you're using frequency in the x-axis when you expect angular frequency. You should modify your x_range
computation like this:
x_range = 2 * np.pi * np.fft.ifftshift(np.fft.fftfreq(n, d=delta_omega))