pythonnumpyfftcontinuous-fourier

How to calculate continuous Fourier inverse numerically with numpy?


I'm trying to implement a function that calculates the inverse Fourier transform of a continuous function F(ω) to obtain f(x).

My input function is continuous (you can assume it is a lambda function) so I have to approximate it discretely using discrete sampling. Here's my current implementation:

import numpy as np

def fourier_inverse(f, n, omega_max):
    """
    :param f: The function to be transformed
    :param n: Number of samples
    :param omega_max: The max frequency we want to be samples
    """
    omega_range = np.linspace(-omega_max, omega_max, n)
    f_values = f(omega_range, sigma, alpha)
    inverse_f = np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f_values)))
    return inverse_f

However, I have two main concerns:

  1. The x-axis of the result ranges from 0 to n, which doesn't seem correct. I don't know how to approach it.
  2. I'm unsure if the overall approach is correct, particularly regarding the sampling and normalization.

Any help or guidance would be greatly appreciated. Thank you!


Solution

  • You can calculate the Fourier transform using:

    np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f(omega_range))))
    
    
    import numpy as np
    
    
    def fourier_inverse(f, n, omega_max):
        omega_range = np.linspace(-omega_max, omega_max, n)
        f_values = np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f(omega_range))))
        delta_omega = 2 * omega_max / n
        x_range = np.fft.ifftshift(np.fft.fftfreq(n, d=delta_omega))
        f_values *= delta_omega * n / (2 * np.pi)
        return x_range, f_values
    
    
    x, f = fourier_inverse(lambda x: np.exp(-x ** 2), 512, 10)
    
    print(x, f)