scipy.ndimage.fourier_gaussian与自定义2D高斯滤波器实现的差异排查

差异根源 #

你通过构造空间域归一化的有限高斯核再做FFT的方式，和scipy.ndimage.fourier_gaussian的实现逻辑本质不同：

• 你的方法先在空间域生成有限大小、求和为1的高斯核，再通过FFT得到传递函数，会受核尺寸截断的影响
• scipy则直接在傅里叶域生成无限大高斯核的传递函数，完全避免了有限核的截断误差，且不需要额外归一化

scipy.ndimage.fourier_gaussian的传递函数定义 #

scipy的实现逻辑可拆解为：

1. 用np.fft.fftfreq(size)生成对应维度的归一化频率坐标（单位：1/像素，范围为[-0.5, 0.5)）
2. 对每个频率分量k，计算一维高斯传递因子：exp(-(2 * np.pi * sigma * k)**2 / 2)
3. 2D传递函数为两个一维因子的外积（利用二维高斯的可分离性）

对应的数学表达式：
$$H(k_x, k_y) = \exp\left(-2\pi^2\sigma2(k_x^2 + k_y^2)\right)$$
其中$k_x, k_y$是归一化频率（$k_x = u/size$，$u$为FFT频率索引的偏移量）

匹配scipy的自定义实现 #

直接在傅里叶域生成传递函数，即可得到和scipy完全一致的结果：

import numpy as np
import scipy

sigma = 2
size = 8
np.random.seed(42)
image = np.random.rand(size, size)
imageTf = np.fft.fft2(image)

# Scipy官方实现
imageBluredTfScipy = scipy.ndimage.fourier_gaussian(imageTf, sigma)

# 自定义匹配实现
def fourier_gaussian_custom(tf, sigma):
    dim_size = tf.shape[0]
    # 生成频率坐标
    freq = np.fft.fftfreq(dim_size)
    # 计算一维高斯传递函数
    gauss_1d = np.exp(-(2 * np.pi * sigma * freq)**2 / 2)
    # 扩展为二维传递函数
    gauss_2d = np.outer(gauss_1d, gauss_1d)
    # 傅里叶域相乘
    return tf * gauss_2d

imageBluredTfCustom = fourier_gaussian_custom(imageTf, sigma)

# 验证差异（误差在机器浮点精度范围内）
print(np.max(np.abs(imageBluredTfCustom - imageBluredTfScipy)))

运行后两者的差异会在$10^{-15}$量级，属于正常浮点误差，结果完全一致。

原实现差异的具体原因 #

你构造的有限高斯核做FFT后，传递函数等价于无限大高斯核的傅里叶变换与矩形窗函数的傅里叶变换的卷积，会引入吉布斯效应；同时空间域的求和归一化仅保证了DC分量（[0,0]点）一致，但高频部分的衰减特性和scipy的无限核传递函数完全不同，因此出现明显差值。

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