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()

I expect the plot to be:

Buy my output is this:

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)?