离散逆傅里叶变换实现中X轴缩放异常问题排查与修正

问题 #

我正在用Python实现离散逆傅里叶变换，近似高斯函数的逆傅里叶变换。输入函数为sqrt(pi) * e^(-w^2/4)，预期输出是e^(-x^2)。目前得到的曲线形状正确，但X轴缩放有问题，结果比预期窄很多。以下是实现代码：

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

问题分析与修正 #

问题出在X轴的缩放计算上：

• np.fft.fftfreq(n, d=delta_omega)返回的是循环频率（单位为1/频率采样间隔），但我们需要的是时域坐标x，必须将其乘以2π才能转换为正确的角频率对应的时域值。因为离散傅里叶变换中，时域坐标x和频率采样间隔Δω的关系是x = 2π * k/(nΔω)，其中k是采样点索引，这正是2π * fftfreq(...)的计算结果。

修正后的fourier_inverse函数如下：

def fourier_inverse(f, n, omega_max):
    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轴缩放：乘以2π转换为正确的时域坐标
    x_range = np.fft.ifftshift(np.fft.fftfreq(n, d=delta_omega)) * 2 * np.pi
    inverse_f *= delta_omega * n / (2 * np.pi)
    return x_range, inverse_f

也可以直接计算X轴范围，避免依赖fftfreq，结果一致：

def fourier_inverse(f, n, omega_max):
    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]
    # 直接计算时域采样间隔和范围
    delta_x = 2 * np.pi / (n * delta_omega)
    x_range = np.fft.ifftshift(np.arange(-n//2, n//2) * delta_x)
    inverse_f *= delta_omega * n / (2 * np.pi)
    return x_range, inverse_f

修正后运行代码，得到的曲线就是标准的e^(-x²)高斯函数，X轴缩放完全匹配预期。

内容的提问来源于stack exchange，提问作者Peyman